Efficient simulation of ultrafast quantum nonlinear optics with matrix product states
Ryotatsu Yanagimoto, Edwin Ng, Logan G. Wright, Tatsuhiro Onodera, Hideo Mabuchi
EEfficient simulation of ultrafast quantum nonlinear optics with matrix product states
Ryotatsu Yanagimoto, ∗ Edwin Ng, Logan G. Wright,
2, 3
Tatsuhiro Onodera,
2, 3 and Hideo Mabuchi E. L. Ginzton Laboratory, Stanford University, Stanford, California 94305, USA NTT Physics and Informatics Laboratories, NTT Research, Inc.,1950 University Ave., East Palo Alto, California 94303, USA School of Applied and Engineering Physics, Cornell University, Ithaca, New York 14853, USA (Dated: February 12, 2021)Ultra-short pulses propagating in nonlinear nanophotonic waveguides can simultaneously leverageboth temporal and spatial field confinement, promising a route towards single-photon nonlinearitiesin an all-photonic platform. In this multimode quantum regime, however, faithful numerical sim-ulations of pulse dynamics na¨ıvely require a representation of the state in an exponentially largeHilbert space. Here, we employ a time-domain, matrix product state (MPS) representation to en-able efficient simulations by exploiting the entanglement structure of the system. In order to extractphysical insight from these simulations, we develop an algorithm to unravel the MPS quantum stateinto constituent temporal supermodes, enabling, e.g., access to the phase-space portraits of arbitrarypulse waveforms. As a demonstration, we perform exact numerical simulations of a Kerr soliton inthe quantum regime. We observe the development of non-classical Wigner-function negativity in thesolitonic mode as well as quantum corrections to the semiclassical dynamics of the pulse. A similaranalysis of χ (2) simultons reveals a unique entanglement structure between the fundamental andsecond harmonic. Our approach is also readily compatible with quantum trajectory theory, allow-ing full quantum treatment of propagation loss and decoherence. We expect this work to establishthe MPS technique as part of a unified engineering framework for the emerging field of broadbandquantum photonics. I. INTRODUCTION
The ability to manipulate photon-photon interactionsat the quantum level holds the key to lifting conventionalclassical limits in a wide range of photonic technolo-gies and applications [1–6]. Recent efforts in the fieldof nonlinear nanophotonics have resulted in the devel-opment of ultra-low-loss and highly efficient platformsfor nonlinear optics [7, 8], with experimental numberscoming remarkably close to bridging the long-standinggap between classical optics and the “strong interactionregime” of quantum optics with single-photon-level non-linearity [9–13]. In particular, advances in dispersion en-gineering on these platforms enable ultra-short-pulse op-eration [14, 15], where the available peak power furtherleverages the material nonlinearities by orders of magni-tude, bringing the possibility of engineering highly non-classical states of light into the foreseeable future [16].In the presence of such strong nonlinearities, the quan-tum nature of individual photons plays a critical role inthe physical behavior of these systems [17], i.e., classi-cal mean-field theories and semiclassical approximationsare no longer valid in predicting the results of exper-iments [18, 19]. At the same time, controllably har-nessing these exotic quantum phenomena for photonicquantum engineering requires the development of uni-fied simulation and modeling methodologies faithful tothe hardware level, posing a significant and imminenttheoretical and modeling problem. In general, model-ing non-classical photon dynamics in a broadband system ∗ Email: [email protected] like a strongly nonlinear propagating pulse is a nontrivialtask due to the immense dimension of the Hilbert spacethat the quantum optical states in general occupy. Forinstance, when we discretize an optical pulse into 100bins, as is conventionally done in classical pulse propa-gation, but with, say, <
10 interacting photons in eachbin, we would na¨ıvely need to compute the evolution of a ∼ -dimensional quantum state vector, which clearlyexceeds any available computational resource. Thus, tofully exploit the technological advantages offered by theseemerging quantum-optical devices, it is essential to ap-ply sophisticated model reduction techniques to the na¨ıvefull quantum model to obtain computationally tractablemodels that nevertheless retain quantum features ex-pected to be important given the hardware specifications.For modeling the propagation of a quantum pulse,some key features that we can utilize to reduce thecomplexity of the na¨ıve quantum state are (1) the one-dimensionality of the optical field, and (2) the locality ofthe interactions and dynamics. Quantum mechanically,a nonlinear waveguide can be seen as a system of bosons(i.e., photons) confined to a one-dimensional space evolv-ing under a Hamiltonian that mediates only local inter-actions, and intuitively, the classical envelope of a pulse(i.e., the g (1) correlation function) simply describes howthe bosons are spatially distributed [20]. Notably, it isknown in the context of many-body physics that sucha one-dimensional system of locally interacting particlesoften exhibits a limited amount of entanglement [21], in-dicating that a large portion of the na¨ıvely large Hilbertspace (i.e., the parts representing highly nonlocal entan-gled states) are not relevant. One technique extensivelyused in the field of many-body physics to exploit this a r X i v : . [ qu a n t - ph ] F e b feature is to cast the quantum state into the form of amatrix product state (MPS) [22, 23], which yields an ef-ficient state representation for numerical studies [24, 25].MPS has recently been extended to analyze opti-cal systems, such as in the treatment of pulse prop-agation through a mesoscopic atomic cloud [26] andatom-coupled chiral waveguides [27]. While these stud-ies have uncovered rich many-body dynamics, the focushas predominantly been on the “particle” aspect of thephysics, such as analyzing their mean-field distributionsand n -particle correlation functions. On the other hand,relatively little attention has been paid to the phase-space properties of the quantum state [28] which char-acterize the “oscillator” aspect of the physics relatedto the continuous-variable (CV) nature of the opticalfield [29, 30]; in many cases, it is this latter picture thatconceptually connects more directly to the perspectivesemployed both in classical wave optics and CV quan-tum information. In the CV approach, it is natural toask about, for example, the reduced density matrix [31]of a given pulse supermode, which is indispensable forphotonic quantum state engineering in the time [32–35] and/or spectral [36, 37] domains, but which is notstraightforward to infer from the n -particle correlationfunctions. Other CV quantities such as the Wigner func-tion and various entanglement measures computed fromthe reduced density matrix are routinely used to quan-tify photonic states as resources for quantum informationprocessing [38, 39].In this article, we first review the application of MPSand MPS-based time evolution methods to the problemof quantum pulse propagation. In order to appreciatethe optical phase-space dynamics of this novel nonlinearquantum regime, we develop a “demultiplexing” schemewhich can be applied to an MPS representation of a quan-tum pulse to extract the reduced density matrix in abasis of pulse supermodes [33]. In our scheme, a care-fully chosen sequence of one- and two-mode (i.e., local)linear operations (phaseshifters and beamsplitters) areused to overcome the problem of manipulating nonlo-cal supermodes—which are not naturally accessible inthe MPS framework—by mapping them onto local bins.As a demonstration, we analyze the quantum propaga-tion of a pulse initialized in a coherent-state Kerr soli-ton of a 1D χ (3) -nonlinear waveguide [40], using theMPS-based approach of time-evolving block decimation(TEBD) [21, 41, 42] which is compatible with standardquantum-optical treatments of linear loss [43]. We showthat the Wigner function of the state contained in thecanonical sech-pulse supermode can exhibit a consider-able amount of negativity with an enhancement corre-lated with the peak intensity. We highlight features of thepulse dynamics that indicate departures from establishedmodels, such as broadband corrections to the conven-tional time-dependent Hartree-Fock approximation forfundamental Kerr solitons [44, 45], and, for higher-ordersolitons [46], stark qualitative deviations from classicalbreather dynamics. Then, we extend our analysis to the quantum propagation of a pulse initialized as a coherent-state simulton [47], a quadratic soliton of a 1D χ (2) -nonlinear waveguide composed of two co-propagatingpulses at the fundamental harmonic (FH) and secondharmonic (SH). We compute a two-mode reduced densitymatrix for the FH and SH pulse supermodes and revealtheir entanglement structure. As a result, we find thatnon-classicality in the quantum simulton primarily accu-mulates in a hybrid supermode consisting of a specificlinear combination of the FH and SH supermodes. Thenumerical techniques highlighted in these examples canbe generalized to address various questions about quan-tum pulse propagation expected to arise with the adventof strongly-interacting broadband quantum photonics. II. MATRIX PRODUCT STATES FORQUANTUM OPTICAL PULSE PROPAGATION
In this section, we introduce basic concepts of matrixproduct states (MPS), together with methods to imple-ment time evolutions. While we base our discussionson quantum pulse propagation on a 1D χ (3) -nonlinearwaveguide to keep discussions concrete, the basic con-cepts introduced in this section is general and can beapplied to a broader class of systems. We extend ouranalysis to χ (2) nonlinear waveguides in Sec. V.The Hamiltonian of a 1D χ (3) -nonlinear waveguide inthe moving frame takes a form [20]ˆ H = − (cid:90) d z (cid:16) ˆ φ † z ∂ z ˆ φ z + ˆ φ † z ˆ φ z (cid:17) , (1)where photon field annihilation operators fulfill commu-tation relationships [ ˆ φ z , ˆ φ † z (cid:48) ] = δ ( z − z (cid:48) ). The mean-field(i.e., c-number) equation under (1) takes a well-knownform of nonlinear Schr¨odinger equation (NLSE) [40]i ∂ t φ z = − ∂ z φ z − | φ z | φ z , (2)where time t and space z have been normalized. In thecontext of many-body physics, (1) is referred to as Lieb-Liniger Hamiltonian, describing bosons (e.g., photons)with point-like interactions [24, 48].While the continuum coordinate z is convenient for an-alytic studies, it is often easier to work in a discretizedcoordinate for numerical evaluations. More importantly,discretization of the field allows us to encode systemstates on an MPS. We consider a finite space interval − L/ ≤ z ≤ L/ N spatial binswith size ∆ z = L/N . As shown in Fig. 1 (a), we assigna photon annihilation operator ˆ a m to the m th the spa-tial bin for m ∈ { , , . . . , N } . When N is large enough,(1) can be approximated by a discretized Bose-Hubbard z a a a a a SPM
Dispersion f f f f f f(z) (a) L Δ z f a a a a a a a SPM SPM SPM SPM SPM SPM
Dispersion Dispersion DispersionDispersion Dispersion
SPM SPM SPM SPM SPM SPM
Dispersion Dispersion DispersionDispersion Dispersion T i m e (b) FIG. 1. (a) Finite space interval − L/ ≤ z ≤ L/ n spatial bins with size ∆ z = L/N . Pulse en-velope function f ( z ) (red lines) is discretized to form a vec-tor f = ( f , . . . , f N ) (cid:124) (orange boxes). Self-phase modulation(SPM) acts locally on discretized photon field (blue arrows),while dispersion mediates interactions between neighboringbins (brown arrows). (b) An example of how time evolutionunder the Hamiltonian (3) can be decomposed into one-modeSPM operations (blue boxes) and two-mode operations cor-responding to dispersions (brown boxes). Spatial bins are la-beled by corresponding annihilation operators for the purposeof illustration. Hamiltonian of the form [49, 50]ˆ H = (cid:88) m (cid:20) − z (cid:0) ˆ a † m ˆ a m +1 + H . c . (cid:1) + 1∆ z ˆ a † m ˆ a m − z ˆ a † m ˆ a m (cid:21) . (3)Similarly, as shown in the figure, a pulse waveform inthe continuous coordinate f ( z ) is discretized as an N -dimensional vector f = ( f , . . . , f N ) (cid:124) .A generic system state can be written as | Ψ (cid:105) = (cid:88) i c i | i (cid:105) , (4)where i = ( i , i , . . . , i N ) (cid:124) and | i (cid:105) = | i (cid:105) ⊗ | i (cid:105) · · · | i N (cid:105) .Each local Hilbert space is spanned by Fock states | i m (cid:105) ( i m ≥ c i re-quire computational resources that grow exponentiallywith respect to N . On the other hand, depending on thespecific properties of the Hamiltonian, such as locality, much of the states in the entire Hilbert space might notbe populated and thus could be excluded. To this end,matrix product states (MPS) allows an effective repre-sentation of quantum states given the amount of entan-glement is limited and local.Using an MPS representation with bond dimension χ ,(4) can be approximated as [21, 41] c i = (cid:88) α ,...,α n − Γ [1] i α λ [1] α Γ [2] i α α λ [2] α · · · λ [ N − α N − Γ [ N ] i N α N − , (5)where Γ [ m ] is a rank-3 tensor, λ [ m ] is a vector, and each α m runs through 1 to χ . Intuitively, (5) is a decompo-sition of a rank- N tensor c i into a product of low-ranktensors. For instance, a coherent pulse with a normalizedenvelope f with (cid:80) m | f m | = 1 takes a formΓ [ m ] i m = e −| f ∗ m α | / ( f ∗ m α ) i m √ i m ! , λ [ m ]1 = 1 , (6)where all the other tensor elements are zero. Physically,(6) is obtained by displacing a supermode ˆ A = (cid:80) m f m ˆ a m by α from the vacuum. Notice that parameters neededfor (5) has a favorable polynomial scaling of O ( χ N ).The bond dimension χ is related to the maximum amountof entanglement that (5) can support, and larger χ isneeded to describe | Ψ (cid:105) with longer-range entanglement.Notably, it is heuristically known that entanglement in1D quantum many-body systems is often limited [21],making MPS an ideal representation to study quantumpulse propagation.Similarly, a generic operatorˆ O = (cid:88) ii (cid:48) O ii (cid:48) (cid:12)(cid:12) i (cid:11)(cid:10) i (cid:48) (cid:12)(cid:12) (7)can also be expressed in the form of a matrix product as O ii (cid:48) = (cid:88) α ,...,α N − O [1] i i (cid:48) α O [2] i i (cid:48) α α · · · O [ N ] i N i (cid:48) N α N − (8)where O [ m ] is a rank-4 tensor. Operators expressed inthe form of (8) are referred to as matrix product oper-ators (MPO), and their expectation values with respectto MPS are computed via tensor contractions [22]. Justin the same way as χ for an MPS is determined by howentangled the state is, bond dimension of an MPO is de-termined by how non-local the operator ˆ O is.In his seminal paper [41], Vidal introduced an algo-rithm, which is often referred to as time-evolving blockdecimation (TEBD), to efficiently simulate the time evo-lution of an MPS. Multiple open-source packages havebeen developed for TEBD [51, 52] to this date. TEBDutilizes the fact that updating an MPS for local one-mode or two-mode unitary operations can be done effi-ciently. While this process involves a truncation of minorsingular-value components, truncation error can in prin-ciple be arbitrary small by taking large enough bond di-mension χ . For instance, as shown in Fig. 1(b), dynamicsunder (3) is simulated by Trotter decomposing a short-time evolution ˆ U = e − i ˆ Hδt into a sequence of one-modeand two-mode operationsˆ U ≈ ( ˆ D ˆ D ˆ D · · · )( ˆ D ˆ D ˆ D · · · )( ˆ S ˆ S ˆ S · · · ) , (9)where ˆ S m = exp (cid:0) i δt ˆ a † m ˆ a m / z (cid:1) implements Kerr self-phase modulation (SPM) on the m th bin, andˆ D m = exp (cid:20) i δt z (ˆ a † m ˆ a m +1 + ˆ a m ˆ a † m +1 − ˆ a † m ˆ a m − ˆ a † m +1 ˆ a m +1 ) (cid:105) (10)represents hopping interactions between m th and ( m +1)th bins due to quadratic energy dispersions. Note thatoperations in each parentheses of (9) commute, and thus,they can be implemented in parallel numerically. Addi-tionally, higher-order decomposition methods can reducethe Trotter discretization errors [53]. For comprehensiveoverview on time-evolution methods for MPS, we leadreaders to Ref. [42].It is worth mentioning that we can readily includedissipations to the simulation, e.g., by the Monte-Carlowavefunction method (MCWF) [43]. The capability ofincluding loss is particularly important for optical simu-lations, not only because realistic optical systems oftenhave a non-negligible amount of loss, but also becausedissipation plays critical roles in a host of emergent phe-nomena, such as dissipative Kerr solitons [54] and thephysics of PT -symmetric systems [55, 56]. III. DEMULTIPLEXING SUPERMODES FROMAN MPS
After a numerical simulation of a quantum pulse prop-agation using MPS, we wish to calculate various physicalproperties of the resultant quantum state | Ψ( t ) (cid:105) . For in-stance, the two photon correlation function g (2) ( (cid:96), m ) = (cid:104) ˆ a † (cid:96) ˆ a † m ˆ a m ˆ a (cid:96) (cid:105) / (cid:104) ˆ a † (cid:96) ˆ a (cid:96) (cid:105)(cid:104) ˆ a † m ˆ a m (cid:105) can be calculated by express-ing ˆ a † (cid:96) ˆ a † m ˆ a m ˆ a (cid:96) and ˆ a † j ˆ a j ( j = (cid:96), m ) as MPOs and comput-ing their expectation values. This procedure can be ex-tended to compute n -particle correlation function as well.Generically, physical quantities associated to MPOs withlarger bond dimension are more expensive to compute.In the context of quantum engineering and informa-tion, it is of great importance to know quantum statespopulating certain supermodes of interest. To be moreprecise, let us consider s (cid:28) N supermodes, where the r th pulse supermode [33] is defined asˆ A ( r ) = N (cid:88) m =1 f ( r ) m ˆ a m (11)for an arbitrary set of orthonormal vectors { f (1) , . . . , f ( s ) } . We denote the Hilbert space forthese supermodes as S = S ⊗ S ⊗ . . . S s , where S r is the space for the r th supermode, and the Hilbert spacefor the rest of the system is denoted as E . Here, we areinterested in a reduced density matrixˆ ρ S = Tr E | Ψ (cid:105)(cid:104) Ψ | = (cid:88) jj (cid:48) ρ jj (cid:48) | j (cid:105) S (cid:10) j (cid:48) (cid:12)(cid:12) S , (12)where j = ( j , . . . , j s ) (cid:124) , | j (cid:105) S = | j (cid:105) S ⊗ | j (cid:105) S ⊗ · · · ⊗| j s (cid:105) S s , and similarly for (cid:12)(cid:12) j (cid:48) (cid:11) S . Once we obtain ˆ ρ S , whosedimension is much smaller than the original Hilbertspace, we can study detailed properties of the state,such as Wigner function and entanglement among super-modes with standard techniques for continuous variablesystems [29].To obtain ˆ ρ S , we ideally would like to have a low-rankMPO representation of the operatorsˆ µ jj (cid:48) = (cid:12)(cid:12) j (cid:48) (cid:11) S (cid:104) j | S ⊗ ˆ E , (13)where ˆ E is an identity operator on E . The expectationvalues of these operators would directly give the matrixelements of ˆ ρ S via (ˆ ρ S ) jj (cid:48) = (cid:104) Ψ | ˆ µ jj (cid:48) | Ψ (cid:105) , but the highlynon-local structure of ˆ µ jj (cid:48) makes it nontrivial to find sucha low-rank MPO expression.In the following, as a main result of this research, wedescribe a procedure to efficiently calculate the reduceddensity matrix ˆ ρ S of an MPS | Ψ (cid:105) for an arbitrary setof supermodes comprising S . As shown schematically inFig. 2, our approach is based on constructing a linearunitary operation ˆ V which “demultiplexes” the s (nonlo-cal) supermodes into the leftmost s (local) spatial bins.In other words, by computing | Φ (cid:105) = ˆ V | Ψ (cid:105) , we would liketo calculate the matrix elements of ˆ ρ S via(ˆ ρ S ) jj (cid:48) = (cid:104) Φ | ˆ π jj (cid:48) | Φ (cid:105) (14)where the operators ˆ π jj (cid:48) have a local formˆ π jj (cid:48) = (cid:32) s (cid:79) m =1 | j (cid:48) m (cid:105)(cid:104) j m | (cid:33) ⊗ (cid:32) N (cid:79) m = s +1 ˆ m (cid:33) , (15)where ˆ m is an identity operator on the m th spatial bin.Explicitly, the MPO representation of ˆ π jj (cid:48) is O [ m ] i m i (cid:48) m = (cid:26) δ i m j (cid:48) m δ i (cid:48) m j m m ≤ sδ i m i (cid:48) m s < m , (16)with all other tensor elements zero; because (16) only re-quires a bond dimension of 1, its expectation value can beefficiently computed. Then the demultiplexing problemposed by (14) is to obtainˆ V † ˆ π jj (cid:48) ˆ V = ˆ µ jj (cid:48) . (17)While solutions to (17) are not generally unique, wechoose to construct ˆ V as a product of local one-modeand two-mode linear operations to allow for the compu-tation of | Φ (cid:105) = ˆ V | Ψ (cid:105) using the standard MPS operations a a a a N-2 a N-1 a N R (1) R (1) R (1) R N-2 (1) R N-1 (1) R (2) R (2) R N-2 (2) R N-1 (2) R (2) A = Σ f m a m R (1) R N (1) R N (2) a a N-2 a N-1 a N (2)(2)(2)(2) (1) A (2) A (1) (1) (1) A = Σ f m a m (2) (2) (2) FIG. 2. Illustration of our supermode demultiplexing scheme.Application of a unitary transformation ˆ R ( r ) , which is com-posed of one-mode and two-mode gates { ˆ R ( r ) r , . . . , ˆ R ( r ) N } , de-multiplexes r th supermode ˆ A ( r ) = (cid:80) m f ( r ) m ˆ a m to the r th (lo-cal) spatial bin. We show operations up to r = 2, wherelocal annihilation operators of the first and second bins aretransformed to ˆ A (1) and ˆ A (2) , respectively. For the purpose ofillustration, spatial bins are labeled by associated annihilationoperators. also used in TEBD. Such construction of a linear unitarygate with one-mode and two-mode operations may beseen as a family of general multiport linear interferom-eters [57, 58], where each “port” is a discretized spatialbin of an MPS in our setup.Our scheme to construct ˆ V works in an iterative man-ner with iteration index r ∈ { , , . . . , s } . In the r thiteration, we construct a transformation ˆ R ( r ) that de-multiplexes the r th supermode into the r th spatial bin.Note that, as a consequence, the construction of ˆ R ( r ) isdependent on the partial operationsˆ V ( r − = ˆ R ( r − ˆ R ( r − · · · ˆ R (1) (18)which have already been applied. If the previous( r −
1) iterations were valid, then we can supposeˆ V ( r − implements the transformations ˆ a m (cid:55)→ ˆ a ( r − m =ˆ V ( r − † ˆ a m ˆ V ( r − , whereˆ a ( r − m = (cid:40) ˆ A ( m ) m ≤ r − (cid:80) N(cid:96) = m − r +1 c ( r − m(cid:96) ˆ a (cid:96) m ≥ r . (19)Here, the first line of (19) means the previous ( r − r −
1) into the leftmost ( r −
1) spatial bins.The second line indicates that for spatial bins with index m ≥ r , the effect of the transformation ˆ V ( r − has beento only “mix in” components from bins m − r + 1 up to N ; more precisely, for m ≥ r , ˆ a ( r − m is independent ofany ˆ a (cid:96) such that (cid:96) ≤ m − r . (19) is fulfilled for the basecase of ˆ V (0) = ˆ and c (0) m(cid:96) = δ m(cid:96) at r = 0, so we only needto ensure that our construction of ˆ R ( r ) preserves (19) forany r >
0. At the end of the final s th iteration, we shouldthen find that ˆ a ( s ) m = ˆ A ( s ) m for 1 ≤ m ≤ s , which satisfiesour goal of (17) and allows us to take ˆ V = ˆ V ( s ) .We parametrize ˆ R ( r ) with two sets of angles { θ ( r ) r , . . . , θ ( r ) N } and { ϕ ( r ) r , . . . , ϕ ( r ) N − } (i.e., the r th iter-ation requires 2( N − r ) + 1 parameters). For m < N ,these parameters are taken to define a set of phase-shifter/beamsplitter operationsˆ R ( r ) m = exp (cid:104) i θ ( r ) m (ˆ a † m ˆ a m − ˆ a † m +1 ˆ a m +1 ) (cid:105) (20a) × exp (cid:104)(cid:16) π − ϕ ( r ) m (cid:17) (cid:16) e − i θ ( r ) m ˆ a † m ˆ a m +1 − e i θ ( r ) m ˆ a m ˆ a † m +1 (cid:17)(cid:105) , which implement the operator transformations [59]ˆ R ( r ) † m ˆ a m ˆ R ( r ) m = e i θ ( r ) m sin ϕ ( r ) m ˆ a m + cos ϕ ( r ) m ˆ a m +1 (20b)ˆ R ( r ) † m ˆ a m +1 ˆ R ( r ) m = e − i θ ( r ) m sin ϕ ( r ) m ˆ a m +1 − cos ϕ ( r ) m ˆ a m . For m = N , we fix ˆ R ( r ) N = e i θ ( r ) N ˆ a † N ˆ a N . We then constructˆ R ( r ) by cascading these one- and two-mode operationsaccording to ˆ R ( r ) = ˆ R ( r ) r ˆ R ( r ) r +1 · · · ˆ R ( r ) N − ˆ R ( r ) N . (20c)Because ˆ R ( r ) should demultiplex the r th supermodeinto the r th spatial bin, we requireˆ a ( r ) r = ˆ A ( r ) = N (cid:88) (cid:96) =1 f ( r ) (cid:96) ˆ a (cid:96) . (21)On the other hand, based on (19) and (20), we haveˆ a ( r ) r = N (cid:88) m = r g ( r ) m ˆ a ( r − m = N (cid:88) (cid:96) =1 min( r − (cid:96),N ) (cid:88) m = r g ( r ) m c ( r − m(cid:96) ˆ a (cid:96) , (22)where g ( r ) m = e i θ ( r ) m sin ϕ ( r ) m (cid:81) m − k = r cos ϕ ( r ) k . By equating(21) and (22), we can solve for the angles via θ ( r ) m = arg I m ϕ ( r ) m = sin − |I m | , (23a)where I m = f ( r ) m − r +1 − (cid:80) m − k = r g ( r ) k c ( r − k,m − r +1 c ( r ) m,m − r +1 (cid:81) m − k = r cos ϕ ( r ) k (23b)for m ∈ { r, r + 1 , . . . , N } . Notice that the right handside of (23a) only depends on θ ( r ) m (cid:48) and ϕ ( r ) m (cid:48) with m (cid:48) < m ,and thus, we can solve the equations starting from m = r towards n iteratively to straightforwardly determine allof θ ( r ) m and ϕ ( r ) m .For the first supermode, we have analytic solutions θ (1) m = arg f (1) m , ϕ (1) m = sin − (cid:12)(cid:12)(cid:12) f (1) m (cid:12)(cid:12)(cid:12)(cid:114) − (cid:80) m − k =1 (cid:12)(cid:12)(cid:12) f (1) k (cid:12)(cid:12)(cid:12) , (24)while we generally need to employ numerical methodsto demultiplex further supermodes. After solving for all θ ( r ) m and ϕ ( r ) m and determining ˆ R ( r ) , it is straightforwardto obtain c ( r ) (cid:96)m and confirm that (19) is fulfilled by theabove construction. The iterative procedure culminatingin ˆ V = ˆ V ( s ) thus demultiplexes all s supermodes intothe leftmost s spatial bins. We can then compute thereduced density matrix ˆ ρ S via (14). IV. QUANTUM PROPAGATION OF KERRSOLITONS
In this section, we study quantum propagation of aKerr soliton. A classical soliton solution for (2) withan average photon number of ¯ n takes a well-known sechform [40, 46] φ (sech) z ( t ) = ¯ n e i¯ n t/ sech ¯ nz . (25)While this solution is exact in the realm of classical opticswhere it may be thought of as specifying the pulse am-plitudes of a coherent state, quantum mechanical effectssuch as squeezing [60, 61], quantum-induced dispersionof pulse envelope [62], and soliton evaporation [63] canbecome pronounced in the regime of large nonlinearitywhere ¯ n becomes small.Conventionally, quantum mechanical soliton solitonscan be treated using the time-dependent Hartree-Fockapproximation (TDHF) [44]. As long as ¯ n (cid:29) φ z (0) approximatelyevolves according to [45] | Ψ( t ) (cid:105) ≈ e − ¯ n (cid:88) n =0 exp (cid:18) i(2 n − ¯ n )¯ nnt (cid:19) α n ˆ A † n n ! | (cid:105) , (26)where α = √ ¯ n , and ˆ A is the annihilation operatorof the classical soliton-pulse supermode [45]. In dis-cretized coordinates, ˆ A = (cid:80) f (sech) m ˆ a m , where f (sech) m ∝ φ (sech) m ∆ z − L/ (0), up to a proportionality constant indepen-dent of m such that f (sech) is normalized. The mainfeature of the approximate TDHF solution is that it isclosed within the subspace S of the soliton supermode f (sech) , so that, e.g., the reduced density matrix of | Ψ( t ) (cid:105) in the subspace S has unit purity throughout the dynam-ics of (26). Nevertheless, due to the Kerr-type nonlinearphase shifts, (26) deviates from a coherent state as itevolves, leading to a variety of interesting phase-spacedynamics [16, 45, 64, 65]. On the other hand, it is diffi-cult to quantify the accuracy or regime of validity of theTDHF due to its non-perturbative nature, and phase-space dynamics of quantum solitons beyond TDHF inthe few-photon regime remain largely unexplored. In thefollowing, we show that our MPS-based scheme serves asa powerful numerical tool to explore this latter regime.To simulate the propagation of a soliton, we initial-ize a coherent-state MPS according to (6) with envelope f (sech) and with α = √ ¯ n . Using TEBD, we simulatethe time evolution of the state under the Hamiltonian(3) with and without linear loss, where the loss dynam-ics are simulated via MCWF with quantum jump opera-tors {√ κ ˆ a , . . . , √ κ ˆ a N } , where κ is the power decay rate.Fig. 3(a) shows the time evolution of the photon densitydistribution (i.e., g (1) correlation function) (cid:104) ˆ φ † z ˆ φ z (cid:105) , foran initial pulse amplitude ¯ n = 3. We see that even onthe level of g (1) , the pulse envelope exhibits dispersionas a function of time [62], reflecting the fact that theinitial classical soliton is not an exact eigenstate of thequantum Hamiltonian; we sometimes refer to such non-classical dispersion as being “quantum induced”.We next consider the reduced density matrix ˆ ρ S for thesech-pulse supermode S using the demultiplexing schemedeveloped in Sec. III with the envelope function f (sech) .Consider first the case without loss, or κ = 0. Fig. 3(b)shows snapshots of the Wigner functions of ˆ ρ S , whichexhibit a substantial amount of Wigner function nega-tivity and signify non-classicality in the state as mightbe expected for single-mode Kerr evolution. However,Fig. 3(c) shows that the purity Tr(ˆ ρ S ) exhibits a mono-tonic decay in time, indicating that the sech-pulse sub-space S is not closed under the dynamics and in fact,coupling between S and the rest of the system E can actas an effective decoherence channel. Actually, due to thenonlinear nature of the dynamics, ˆ ρ S may not be pure forany choice of a supermode f in general. Also shown inFig. 3(c) is the volume of the Wigner function negativity,which serves as a measure of the non-classicality of thestate [66]. Following an initial increase, the volume ofWigner function negativity starts to decrease after sometime, again due to competition between the nonlineardynamics and the effective decoherence caused by entan-glement with E . These features are in stark contrast tothe single-mode dynamics predicted by TDHF.We can also contrast this effective decoherence due toentanglement between S and E with standard dissipationdue to linear loss. When loss is incorporated to the sim-ulation, we obtain an ensemble of quantum trajectories {| Ψ ( t ) (cid:105) , . . . , | Ψ M ( t ) (cid:105)} via the MCWF method [43]. Thefinal reduced density matrix is calculated by averagingthe reduced density matrices over the ensemble accordingto ˆ ρ S = M − (cid:80) Mi =1 ˆ ρ S ,i , where ˆ ρ S ,i is the reduced den-sity matrix of | Ψ i (cid:105) . As expected, the linear loss causes Position z P h o t o n d e n s i t y = . (a) (b) (c)(d) t = 0.0 t = 0.5 t = 1.0 t = 1.5 z P h o t o n d e n s i t y = . p t = 0.5 t = 1.0 x p x t S t a t e p u r i t y i n s o li t o n m o d e = 0.0= 0.5 t W F n e g a t i v i t y n W F n e g a t i v i t y t = 0.1 t = 0.2 t = 0.3 FIG. 3. (a)–(c) Quantum simulations of a pulse instantiated as a coherent-state sech soliton using MPS with bond dimension χ = 40. Average photon number of ¯ n = 3 . κ = 0 .
5) and without ( κ = 0) linear loss. MCWF with M = 100quantum trajectories is used for the simulation with finite loss. (a) Time evolution of photon density distribution (cid:104) ˆ φ † z ˆ φ z (cid:105) . (b)Wigner functions for the reduced density matrix ˆ ρ S for the sech supermode f sech . Upper and lower rows are for κ = 0 . κ = 0 .
5, respectively. (c) Time evolution of the purity Tr(ˆ ρ S ) (left figure) and the doubled volume of the Wigner function(WF) negativity [66] (lower figure) of the soliton pulse state. (d) WF negativity of the soliton pulse state at various times asa function of ¯ n . MPS-based simulation with χ = 50 is used. a decay in the amplitude of the photon density distribu-tion as shown in Fig. 3(a), while quantum mechanically,Fig. 3(b) and (c) show that the non-classical features ofthe state are critically diminished by the presence of thelinear loss.We additionally investigate how the amount of exci-tation affects the phase-space dynamics of pulse propa-gation. Classically, due to the nonlinear nature of theinteractions, the effective nonlinear rate is expected tobe enhanced when the peak pulse intensity is larger. InFig. 3(d), we show the volume of the Wigner functionnegativity as a function of the average photon number ¯ n in the soliton pulse, where we observe that larger pulseexcitations increase the rate at which non-classical fea-tures are formed, confirming the classical intuition in thefew-photon regime.Finally, to highlight the difference between full numer-ical simulation and TDHF, Fig. 4 compares the Wignerfunctions of a soliton pulse initialized with ¯ n = 6 ob-tained by the MPS-based simulation against TDHF (26).While they exhibit qualitatively similar interference pat-terns, the discrepancies in the time at which a similarphase-space structure is reached indicate that TDHF isoverestimating the rate of the nonlinear phase-shift (see t = 0 . t = 0 .
15 forTDHF, for instance). Additionally, full quantum sim-ulation results exhibit visible reduction in the magnitudeof Wigner function negativity compared to TDHF, high-lighting the effects of aforementioned decoherence dueto entanglement between S and E . These discrepan-cies point to the presence of broadband physics beyondTDHF in soliton propagation.While the quasi-stationary nature of the quantum dy- p M P S s i m u l a t i o n t = 0.05 t = 0.15 t = 0.3 x p T D H F a pp r o x . x x FIG. 4. Wigner functions of soliton pulse state with ¯ n =6 obtained using MPS-based full-quantum simulation with χ = 50 (upper row) and the time-dependent Hartree-Fockapproximation (TDHF) following (26) (lower row). namics of fundamental solitons is in qualitative agree-ment with their classical behavior, it is in fact possible forhighly-quantum photon dynamics to exhibit much morestriking deviations from classical behavior. This is thecase, for example, if we apply our simulation method tostudy the quantum propagation of a second-order solitonwith classical waveform [46] φ (2sech) z ( t ) = 2 e i¯ n t/ ¯ n (cid:16) e i¯ n t cosh(¯ nz/
2) + cosh(3¯ nz/ (cid:17) n t ) + 4 cosh(¯ nz ) + cosh(2¯ nz ) , (27)which, classically, is a periodic “breather” solution of the (b) Full quantum simulation(a) Classical solution FIG. 5. Time evolution of the photon density distribution ofa pulse instantiated as a coherent-state second-order solitonwith mean photon number ¯ n (2sech) = 8. (a) Classically ex-pected dynamics given by (27). (b) Full quantum simulationunder the Hamiltonian (1) using MPS with bond dimension χ = 60. NLSE (2). The average photon number of the second-order soliton is ¯ n (2sech) = 4¯ n , where ¯ n is the averagephoton number in its corresponding fundamental soliton.At t = 0, this means the waveform field amplitude of thesecond-order soliton is twice that of the fundamental soli-ton, i.e., φ (2sech) z (0) = 2 φ (sech) z (0). After t = 0, as shownin Fig. 5(a), the classical waveform of the second-ordersoliton exhibits significant narrowing and a characteris-tic triplet structure. On the other hand, as shown inFig. 5(b), full quantum evolution of a second-order soli-ton instantiated in a few-photon coherent state of (27)at t = 0 exhibits qualitatively different dynamics. Whilethe photon density distribution exhibits some narrow-ing of the pulse width initially, the peak pulse intensityfails to reach the level expected from the classical solu-tion. Moreover, no signature of the triplet structure isobserved. We attribute these features to the quantum-induced dispersion of the pulse envelope that we also ob-served in the quantum propagation of fundamental soli-tons [62], which appears to play a more critical role inthe evolution of this higher-order soliton. V. QUANTUM PROPAGATION OFSIMULTONS
In this section, we apply our technique to a 1D χ (2) nonlinear waveguide in which interactions occur betweena fundamental harmonic (FH) band and a second har-monic (SH) band. After normalization with respectto time and space, the system Hamiltonian takes theform [67, 68]ˆ H = − (cid:90) d z (cid:16) ˆ φ † z ∂ z ˆ φ z + β ˆ ψ † z ∂ z ˆ ψ z (cid:17) + 12 (cid:90) d z (cid:16) ˆ φ † z ˆ ψ z + ˆ φ z ˆ ψ † z (cid:17) , (28)where ˆ φ z and ˆ ψ z are respectively FH and SH local fieldannihilation operators with commutation relationships a b a b a b a b a b U NL,1 U NL,2 U NL,3 U NL,4 U NL,5 U SWAP U SWAP U SWAP a b a b a b a b a b U A,1 U A,3 U B,2 U B,4 U SWAP U SWAP U SWAP U SWAP U SWAP a b a b a b a b a b U B,1 U B,3 U A,2 U A,4 U SWAP U SWAP a b a b a b a b a b FIG. 6. An example implementation of a short-time evolu-tion under the χ (2) Hamiltonian (28) using TEBD. DiscretizedFH modes { ˆ a , . . . , ˆ a N } and SH modes { ˆ b , . . . , ˆ b N } are en-coded on an MPS in an alternating manner. Swap operationsˆ U SWAP are used to bring distant modes together. Two-modeoperations ˆ U a ,m , ˆ U b ,m , and ˆ U NL ,m are for FH dispersion, SHdispersion, and nonlinear three-wave mixing interaction, re-spectively. [ ˆ φ z , ˆ φ † z (cid:48) ] = [ ˆ ψ z , ˆ ψ † z (cid:48) ] = δ ( z − z (cid:48) ). We have assumed thatthe FH and SH carriers are group-velocity matched, while β represents the group velocity dispersion of SH relativeto FH. (28) can be discretized in space to giveˆ H = (cid:88) m (cid:16) ˆ H a ,m + ˆ H b ,m + ˆ H NL ,m (cid:17) , (29a)withˆ H a ,m = − z (cid:16) ˆ a † m +1 ˆ a m + ˆ a m ˆ a † m +1 − a † m ˆ a m (cid:17) (29b)ˆ H b ,m = − β z (cid:16) ˆ b † m +1 ˆ b m + ˆ b m ˆ b † m +1 − b † m ˆ b m (cid:17) (29c)ˆ H NL ,m = 12 √ ∆ z (cid:16) ˆ a † m ˆ b m + ˆ a m ˆ b † m (cid:17) , (29d)where ˆ a m and ˆ b m are the FH and SH field annihilationoperators for the m th spatial bin, respectively.The presence of both FH and SH fields requires somecare in applying the two-mode operations for implement-ing TEBD for the Hamiltonian (29). One approach isshown schematically in Fig. 6, where we prepare 2 N MPS bins to represent N spatial bins, with FH andSH modes encoded in an alternating manner. To ap-ply Trotterization as in Sec. II, a short-time unitary evo-lution e − i ˆ Hδt is decomposed into two-mode operationsˆ U λ,m = e − i ˆ H λ,m δt ( λ = a , b , NL) and applied as shown inFig. 6. Since ˆ a m and ˆ a m +1 are no longer next to eachother in this representation, we also utilize a swap op-eration ˆ U SWAP [41] to bring these modes together in analternating manner for the application of ˆ U a ,m (and sim-ilarly for ˆ U b ,m ).The Hamiltonian (28) supports various classical solitonsolutions [69, 70]. Specifically, an analytic “simulton”solution exists for β = 2 with the form [47] φ z ( t ) = φ sech (cid:16)(cid:112) φ / z (cid:17) e i φ t/ (30a) ψ z ( t ) = − φ (cid:16)(cid:112) φ / z (cid:17) e φ t/ , (30b)where φ is related to the average FH photon number ¯ n via φ = (cid:112) n /
32. As was done for (25), we discretizeand normalize (30) at t = 0 to construct FH and SHsupermodes with annhilation operators denoted ˆ A andˆ B , respectively. The initial simulton MPS state is thatof a coherent state in ˆ A and ˆ B with displacements √ ¯ n and −√ ¯ n/
2, respectively. In Fig. 7(a), we show the timeevolution via TEBD of the photon density distributionof a pulse initialized in the classical simulton supermode.As for the case of the Kerr soliton, the dynamics showa quantum-induced dispersion of the pulse envelope notpredicted by classical dynamics; as part of this process,we also numerically observe a slight exchange of excita-tions between FH and SH.To investigate the joint phase-space properties of theFH and SH pulse supermodes, we calculate two-modereduced density matrix ˆ ρ S where S = A ⊗ B , the jointHilbert space of FH and SH, respectively. In general,the state of the system features entanglement between A and B ; to quantify this entanglement, we utilize thenegativity measure N (ˆ ρ ) (not to be confused with theWigner function negativity used earlier) [71], defined fora bipartite density matrix ˆ ρ as N (ˆ ρ ) = (1 − (cid:107) ˆ ρ T A (cid:107) ) , (31)where T A is partial transposition with respect to the firstmode, and (cid:107)·(cid:107) is the trace norm. Here, N (ˆ ρ ) > N by con-sidering hybrid mixtures of the modes A and B . Findingthe linear combination that best “disentangles” A and B therefore provides insight into their entanglement struc-ture. More specifically, we introduce a two-mode linearoperation ˆ W (Φ , Θ) = exp (cid:110) Φ (cid:16) e iΘ ˆ A † ˆ B − e − iΘ ˆ A ˆ B † (cid:17)(cid:111) with − π/ ≤ Φ ≤ π/ − π/ ≤ Θ ≤ π/
2, whichimplements the transformation ˆ ρ (cid:48)S = W † ˆ ρ S ˆ W . Impor-tantly, the transformation ˆ W = ˆ W (Φ , Θ ) that mini-mizes N (ˆ ρ (cid:48)S ) is expected to maximize the amount of in-formation available in the single-mode reduced densitymatrices ˆ ρ (cid:48)A = Tr B (ˆ ρ (cid:48)S ) and ˆ ρ (cid:48)B = Tr A (ˆ ρ (cid:48)S ).In Fig. 7(b), for the final state of the pulse prop-agation at t = 4, we map the negativity N (ˆ ρ (cid:48)S ) forvarious (Φ , Θ), which shows a clear minimum at themarked position of (Φ , Θ ). In Fig. 7(c), we showthe Wigner functions of the single-mode reduced densitymatrices before and after the transformation ˆ W . Be-fore applying the transformation, the Wigner functions (a) (b)(c) FIG. 7. Quantum simulations of pulses instantiated as acoherent-state simulton with FH photon number ¯ n = 4. MPSwith bond dimension χ = 50 is used. (a) Time evolutionof the FH and SH photon density distribution. (b) Negativ-ity N ( ρ (cid:48)S ) calculated for the two-mode reduced density ma-trix under linear transformations ˆ ρ (cid:48)S = ˆ W † ˆ ρ S ˆ W , where ˆ W isparametrized by a pair of angles (Φ , Θ). White cross indicates(Φ , Θ ) ≈ (0 . π, . π ) that minimizes N (ˆ ρ (cid:48)S ) for t = 4. (c)Wigner function of the single-mode reduced density matrixbefore (ˆ ρ A , B ) and after (ˆ ρ (cid:48)A , B ) applying ˆ W (Φ , Θ ) at t = 4.Labels A and B represent FH and SH, respectively. of both ˆ ρ A and ˆ ρ B are crescent-shaped with no nega-tivity, indicating that they are highly mixed states dueto the entanglement between FH and SH. On the otherhand, after application of ˆ W , the Wigner function of ˆ ρ (cid:48)A remarkably exhibits considerable non-classicality, whilethe Wigner function of ˆ ρ (cid:48)B resembles that of a coher-ent state. This reveals a somewhat surprising featureof the simulton quantum dynamics: a hybrid supermodeˆ A = cos Φ ˆ A + e iΘ sin Φ ˆ B , composed of both FH andSH components, is the one which predominantly experi-ences strongly nonlinear dynamics, while the other hy-brid supermode ˆ B = cos Φ ˆ B − e − iΘ sin Φ ˆ A experi-ences little nonlinearity. A similar analysis can, in prin-ciple, be applied to a broader class of solitons and generalpulse propagation [46].0 VI. CONCLUSION
In this research, we have motivated the use of MPStechniques to efficiently represent and simulate a quan-tum optical pulse as it dynamically propagates througha nonlinear 1D waveguide. In doing so, we have devel-oped a numerical method to overcome the problem ofefficiently accessing and manipulating nonlocal pulse su-permodes of the local MPS representation, allowing usto view for the first time the full quantum dynamics ofthe pulse in a phase-space picture. As a demonstration,we have performed quantum simulations of Kerr solitonpropagation and observed that the phase-space portraitsof an initially classical sech-pulse supermode can evolvehighly non-classical features, i.e., Wigner function nega-tivity. These results have been contrasted with predic-tions based on TDHF, highlighting the presence of richquantum dynamics beyond conventional approximationsfor quantum Kerr solitons. We have also extended ouranalysis to the quantum propagation of a χ (2) simultonand have revealed unexpected entanglement structure be-tween the FH and SH pulses of the simulton, identifyinga hybrid supermode that predominantly exhibits non-classical features. Our scheme is compatible with localdissipation associated with, e.g., waveguide losses, and,more generally, could be applied to any one-dimensionalphotonic system in principle. Considering the rapidrecent progress towards single-photon nonlinearities indispersion-engineered and highly nonlinear nanophotonicplatforms, it is of imminent interest to establish a uni-fied theoretical framework in which to understand thequantum dynamics of photons in such devices. To this end, our work takes a step towards bridging the signifi-cant conceptual gaps between classical wave optics, CVphotonic quantum information, and strongly-interactingquantum many-body physics, all of which are expected toplay important roles in conceptualizing and engineeringthe future of broadband quantum optics. FUNDING
Army Research Office (W911NF-16-1-0086); NationalScience Foundation (CCF-1918549, PHY-2011363).
ACKNOWLEDGMENTS
R. Y. developed the numerical techniques, performedthe simulations, and generated the figures. E. N. andH. M. advised and directed the project. R. Y. and E. N.wrote the manuscript with detailed input and feedbackfrom all authors. All authors contributed significantly tothe conception of the project.The authors wish to thank NTT Research for theirfinancial and technical support. R. Y. would like tothank Tomohiro Soejima for helpful discussions. R. Y.is supported by Stanford Q-FARM Ph.D. Fellowship andMasason Foundation.
DISCLOSURES
The authors declare no conflicts of interest. [1] The LIGO Scientific Collaboration, Enhanced sensitiv-ity of the LIGO gravitational wave detector by usingsqueezed states of light, Nat. Photon. , 613 (2013).[2] M. Tsang, R. Nair, and X.-M. Lu, Quantum Theoryof Superresolution for Two Incoherent Optical PointSources, Phys. Rev. X , 031033 (2016).[3] N. Gisin and R. Thew, Quantum communication, Nat.Photon. , 165 (2007).[4] J. L. O’Brien, A. Furusawa, and J. Vuˇckovi´c, Photonicquantum technologies, Nat. Photon. , 687 (2009).[5] D. E. Chang, V. Vuleti´c, and M. D. Lukin, Quantumnonlinear optics - photon by photon, Nat. Photon. , 685(2014).[6] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C.Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu,X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan,G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu,and J.-W. Pan, Quantum computational advantage usingphotons, Science (2020).[7] M. Zhang, C. Wang, R. Cheng, A. Shams-Ansari, andM. Lonˇcar, Monolithic ultra-high-Q lithium niobate mi-croring resonator, Optica , 1536 (2017).[8] C. Wang, C. Langrock, A. Marandi, M. Jankowski,M. Zhang, B. Desiatov, M. M. Fejer, and M. Lonˇcar, Ultrahigh-efficiency wavelength conversion in nanopho-tonic periodically poled lithium niobate waveguides, Op-tica , 1438 (2018).[9] J. Lu, M. Li, C.-L. Zou, A. Al Sayem, and H. X. Tang,Towards 1% single photon nonlinearity with periodically-poled lithium niobate microring resonators, Optica ,1654 (2020).[10] M. Placke and S. Ramelow, Engineering AlGaAs-on-insulator towards quantum optical applications, Opt.Lett. , 6763 (2020).[11] S. Ramelow, A. Farsi, Z. Vernon, S. Clemmen, X. Ji,J. E. Sipe, M. Liscidini, M. Lipson, and A. L. Gaeta,Strong Nonlinear Coupling in a Si N Ring Resonator,Phys. Rev. Lett. , 153906 (2019).[12] M. Heuck, K. Jacobs, and D. R. Englund, Photon-photoninteractions in dynamically coupled cavities, Phys. Rev.A , 042322 (2020).[13] A. W. Bruch, X. Liu, J. B. Surya, C.-L. Zou, and H. X.Tang, On-chip χ (2) microring optical parametric oscilla-tor, Optica , 1361 (2019).[14] M. Jankowski, C. Langrock, B. Desiatov, A. Marandi,C. Wang, M. Zhang, C. R. Phillips, M. Lon˘car, and M. M.Fejer, Ultrabroadband nonlinear optics in nanophotonicperiodically poled lithium niobate waveguides, Optica ,
40 (2020).[15] L. Zhang, Q. Lin, Y. Yue, Y. Yan, R. G. Beausoleil,and A. E. Willner, Silicon waveguide with four zero-dispersion wavelengths and its application in on-chipoctave-spanning supercontinuum generation, Opt. Ex-press , 1685 (2012).[16] R. Yanagimoto, T. Onodera, E. Ng, L. G. Wright, P. L.McMahon, and H. Mabuchi, Engineering a Kerr-BasedDeterministic Cubic Phase Gate via Gaussian Opera-tions, Phys. Rev. Lett. , 240503 (2020).[17] K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E.Northup, and H. J. Kimble, Photon Blockade in an Op-tical Cavity With One Trapped Atom, Nature , 87(2005).[18] J. Javanainen and J. Ruostekoski, Light propagation be-yond the mean-field theory of standard optics, Opt. Ex-press , 993 (2016).[19] R. Yanagimoto, E. Ng, M. P. Jankowski, T. Onodera,M. M. Fejer, and H. Mabuchi, Broadband ParametricDownconversion as a Discrete-Continuum Fano Interac-tion (2020), arXiv:2009.01457.[20] P. D. Drummond and M. Hillery, The Quantum Theoryof Nonlinear Optics (Cambridge University Press, 2014).[21] G. Vidal, Efficient Simulation of One-Dimensional Quan-tum Many-Body Systems, Phys. Rev. Lett. , 040502(2004).[22] U. Schollw¨ock, The density-matrix renormalization groupin the age of matrix product states, Ann. Phys. , 96(2011).[23] R. Or´us, A practical introduction to tensor networks:Matrix product states and projected entangled pairstates, Ann. Phys. , 117 (2014).[24] D. Muth and M. Fleischhauer, Dynamics of Pair Correla-tions in the Attractive Lieb-Liniger Gas Dominik, Phys.Rev. Lett. , 150403 (2010).[25] A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller,Measuring Entanglement Growth in Quench Dynamicsof Bosons in an Optical Lattice, Phys. Rev. Lett. ,020505 (2012).[26] M. T. Manzoni, D. E. Chang, and J. S. Douglas, Simu-lating quantum light propagation through atomic ensem-bles using matrix product states, Nat. Commun. , 1743(2017).[27] S. Mahmoodian, G. Calaj´o, D. E. Chang, K. Hammerer,and A. S. Sørensen, Dynamics of Many-Body PhotonBound States in Chiral Waveguide QED, Phys. Rev. X , 031011 (2020).[28] K. E. Cahill and R. J. Glauber, Density Operators andQuasiprobability Distributions, Phys. Rev. , 1882(1969).[29] S. L. Braunstein and P. van Loock, Quantum informa-tion with continuous variables, Rev. Mod. Phys. , 513(2005).[30] D. Gottesman, A. Kitaev, and J. Preskill, Encoding aqubit in an oscillator, Phys. Rev. A , 012310 (2001).[31] M. A. Nielsen and I. L. Chuang, Quantum Computationand Quantum Information (Cambridge University Press,2000).[32] P. C. Humphreys, B. J. Metcalf, J. B. Spring, M. Moore,X.-M. Jin, M. Barbieri, W. S. Kolthammer, and I. A.Walmsley, Linear Optical Quantum Computing in a Sin-gle Spatial Mode, Phys. Rev. Lett. , 150501 (2013).[33] B. Brecht, D. V. Reddy, C. Silberhorn, and M. G.Raymer, Photon Temporal Modes: A Complete Frame- work for Quantum Information Science, Phys. Rev. X ,041017 (2015).[34] W. Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoen-sombutamon, H. Emura, R. N. Alexander, S. Takeda,J. Yoshikawa, N. C. Menicucci, H. Yonezawa, and A. Fu-rusawa, Generation of time-domain-multiplexed two-dimensional cluster state, Science , 373 (2019).[35] V. Ansari, J. M. Donohue, B. Brecht, and C. Silberhorn,Tailoring nonlinear processes for quantum optics withpulsed temporal-mode encodings, Optica , 534 (2018).[36] J. M. Lukens and P. Lougovski, Frequency-encoded pho-tonic qubits for scalable quantum information processing,Optica , 8 (2017).[37] J. Roslund, R. M. de Ara´ujo, S. Jiang, C. Fabre, andN. Treps, Wavelength-multiplexed quantum networkswith ultrafast frequency combs, Nat. Photon. , 109(2014).[38] E. Chitambar and G. Gour, Quantum resource theories,Rev. Mod. Phys. , 025001 (2019).[39] F. Albarelli, M. G. Genoni, M. G. A. Paris, and A. Fer-raro, Resource theory of quantum non-Gaussianity andWigner negativity, Phys. Rev. A , 052350 (2018).[40] G. P. Agrawal, Nonlinear Fiber Optics, 6th edition (Aca-demic Press, 2019).[41] G. Vidal, Efficient Classical Simulation of Slightly Entan-gled Quantum Computers, Phys. Rev. Lett. , 147902(2003).[42] J. J. Garc´ıa-Ripoll, Time evolution of Matrix ProductStates, New J. Phys. , 305 (2006).[43] H. W. Wiseman and G. J. Milburn, Quantum Measure-ment and Control (Cambridge University Press, 2009).[44] A. Haus and Y. Lai, Quantum theory of solitons in opticalfibers. I.Time-dependent Hartree approximation, Phys.Rev. A , 844 (1989).[45] E. W. Wright, Quantum theory of soliton propagation inan optical fiber using the Hartree approximation, Phys.Rev. A , 3836 (1991).[46] Y. S. Kivshar and G. P. Agrawal, Optical Solitons (Aca-demic Press, 2003).[47] M. J. Werner and P. D. Drummond, Simulton solutionsfor the parametric amplifier, J. Opt. Soc. Am. B , 2390(1993).[48] E. H. Lieb and W. Liniger, Exact Analysis of an Interact-ing Bose Gas. I. The General Solution and the GroundState, Phys. Rev. , 1605 (1963).[49] D. Muth, B. Schmidt, and M. Fleischhauer, Fermioniza-tion dynamics of a strongly interacting one-dimensionalBose gas after an interaction quench, New J. Phys. ,083065 (2010).[50] D. Muth, M. Fleischhauer, and B. Schmidt, Discretizedversus continuous models of p -wave interacting fermionsin one dimension, Phys. Rev. A , 013602 (2010).[51] D. Jaschke, M. L. Wall, and L. D. Carr, Open source Ma-trix Product States: Opening ways to simulate entangledmany-body quantum systems in one dimension, Comput.Phys. Commun , 59 (2018).[52] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin,J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull,S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov,D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcol-let, G. Pawlowski, J. D. Picon, L. Pollet, E. Santos,V. W. Scarola, U. Schollw¨ock, C. Silva, B. Surer, S. Todo,S. Trebst, M. Troyer, M. L. Wall, P. Werner, and S. Wes-sel, The ALPS project release 2.0: open source software for strongly correlated systems, J. Stat. Mech.: TheoryExp. , 05001.[53] A. T. Sornborger and E. D. Stewart, Higher-order meth-ods for simulations on quantum computers, Phys. Rev.A , 1956 (1999).[54] T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L.Gorodetsky, Dissipative Kerr solitons in optical microres-onators, Science , 567 (2018).[55] R. El-Ganainy, K. G. Makris, D. N. Christodoulides,and Z. H. Musslimani, Theory of coupled optical PT -symmetric structures, Opt. Lett. , 2632 (2007).[56] N. V. Alexeeva, I. V. Barashenkov, A. A. Sukhorukov,and Y. S. Kivshar, Optical solitons in PT -symmetricnonlinear couplers with gain and loss, Phys. Rev. A ,063837 (2012).[57] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S.Kolthammer, and I. A. Walmsley, Optimal design for uni-versal multiport interferometers, Optica , 1460 (2016).[58] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani,Experimental realization of any discrete unitary opera-tor, Phys. Rev. Lett. , 58 (1994).[59] S. Olivares, Quantum optics in the phase space, Eur.Phys. J. Special Topics , 3 (2012).[60] H. A. Haus and Y. Lai, Quantum theory of soliton squeez-ing: a linearized approach, J. Opt. Soc. Am. B , 386(1990).[61] S. J. Carter, P. D. Drummond, M. D. Reid, and R. M.Shelby, Squeezing of quantum solitons, Phys. Rev. Lett. , 1841 (1987). [62] Y. Lai and H. A. Haus, Quantum theory of solitons inoptical fibers. II. Exact solution, Phys. Rev. A , 854(1989).[63] L. Di Mauro Villari, D. Faccio, F. Biancalana, andC. Conti, Quantum soliton evaporation, Phys. Rev. A , 043859 (2018).[64] N. Korolkova, R. Loudon, G. Gardavsky, M. W. Hamil-ton, and G. Leuchs, Time evolution of a quantum solitonin a kerr medium, J. Mod. Opt. , 1339 (2001).[65] F. Singer, M. J. Potasek, J. M. Fang, and M. C. Teich,Femtosecond solitons in nonlinear optical fibers: Classi-cal and quantum effects, Phys. Rev. A , 4192 (1992).[66] A. Kenfack and K. ˙Zyczkowski, Negativity of the Wignerfunction as an indicator of non-classicality, J. Opt. B:Quantum Semiclass. Opt. , 396 (2004).[67] P. D. Drummond and H. He, Optical mesons, Phys. Rev.A , R1107(R) (1997).[68] M. G. Raymer, P. D. Drummond, and S. J. Carter,Limits to wideband pulsed squeezing in a traveling-waveparametric amplifier with group-velocity dispersion, Opt.Lett. , 1189 (1991).[69] A. V. Buryak and Y. S. Kivshar, Solitons due to secondharmonic generation, Phys. Lett. A , 407 (1995).[70] A. V. Buryak, P. Trapani, D. V. Skryabin, and S. Trillo,Optical solitons due to quadratic nonlinearities: from ba-sic physics to futuristic applications, Phys. Rep. , 63(2002).[71] G. Vidal and R. F. Werner, Computable measure of en-tanglement, Phys. Rev. Lett.65