Efficient simulation of ultrafast quantum nonlinear optics with matrix product states

Ultra-short pulses propagating in nonlinear nanophotonic waveguides can simultaneously leverage both temporal and spatial field confinement, promising a route towards single-photon nonlinearities in an all-photonic platform. In this multimode quantum regime, however, faithful numerical simulations of pulse dynamics na\"ively require a representation of the state in an exponentially large Hilbert space. Here, we employ a time-domain, matrix product state (MPS) representation to enable efficient simulations by exploiting the entanglement structure of the system. In order to extract physical insight from these simulations, we develop an algorithm to unravel the MPS quantum state into constituent temporal supermodes, enabling, e.g., access to the phase-space portraits of arbitrary pulse waveforms. As a demonstration, we perform exact numerical simulations of a Kerr soliton in the quantum regime. We observe the development of non-classical Wigner-function negativity in the solitonic mode as well as quantum corrections to the semiclassical dynamics of the pulse. A similar analysis of $\chi^{(2)}$ simultons reveals a unique entanglement structure between the fundamental and second harmonic. Our approach is also readily compatible with quantum trajectory theory, allowing full quantum treatment of propagation loss and decoherence. We expect this work to establish the MPS technique as part of a unified engineering framework for the emerging field of broadband quantum photonics.


I. INTRODUCTION
The ability to manipulate photon-photon interactions at the quantum level holds the key to lifting conventional classical limits in a wide range of photonic technologies and applications [1][2][3][4][5][6].Recent efforts in the field of nonlinear nanophotonics have resulted in the development of ultra-low-loss and highly efficient platforms for nonlinear optics [7,8], with experimental numbers coming remarkably close to bridging the long-standing gap between classical optics and the "strong interaction regime" of quantum optics with single-photon-level nonlinearity [9][10][11][12][13].In particular, advances in dispersion engineering on these platforms enable ultra-short-pulse operation [14,15], where the available peak power further leverages the material nonlinearities by orders of magnitude, bringing the possibility of engineering highly nonclassical states of light into the foreseeable future [16].In the presence of such strong nonlinearities, the quantum nature of individual photons plays a critical role in the physical behavior of these systems [17], i.e., classical mean-field theories and semiclassical approximations are no longer valid in predicting the results of experiments [18,19].At the same time, controllably harnessing these exotic quantum phenomena for photonic quantum engineering requires the development of unified simulation and modeling methodologies faithful to the hardware level, posing a significant and imminent theoretical and modeling problem.In general, modeling non-classical photon dynamics in a broadband system like a strongly nonlinear propagating pulse is a nontrivial task due to the immense dimension of the Hilbert space that the quantum optical states in general occupy.For instance, when we discretize an optical pulse into 100 bins, as is conventionally done in classical pulse propagation, but with, say, < 10 interacting photons in each bin, we would naïvely need to compute the evolution of a ∼ 10 100 -dimensional quantum state vector, which clearly exceeds any available computational resource.Thus, to fully exploit the technological advantages offered by these emerging quantum-optical devices, it is essential to apply sophisticated model reduction techniques to the naïve full quantum model to obtain computationally tractable models that nevertheless retain quantum features expected to be important given the hardware specifications.
For modeling the propagation of a quantum pulse, some key features that we can utilize to reduce the complexity of the naïve quantum state are (1) the onedimensionality of the optical field, and (2) the locality of the 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 evolving under a Hamiltonian that mediates only local interactions, and intuitively, the classical envelope of a pulse (i.e., the g (1) correlation function) simply describes how the bosons are spatially distributed [20].Notably, it is known in the context of many-body physics that such a one-dimensional system of locally interacting particles often exhibits a limited amount of entanglement [21], indicating that a large portion of the naïvely large Hilbert space (i.e., the parts representing highly nonlocal entangled states) are not relevant.One technique extensively used in the field of many-body physics to exploit this feature is to cast the quantum state into the form of a matrix product state (MPS) [22,23], which yields an efficient state representation for numerical studies [24,25].
MPS has recently been extended to analyze optical systems, such as in the treatment of pulse propagation through a mesoscopic atomic cloud [26] and atom-coupled chiral waveguides [27].While these studies have uncovered rich many-body dynamics, the focus has predominantly been on the "particle" aspect of the physics, such as analyzing their mean-field distributions and n-particle correlation functions.On the other hand, relatively little attention has been paid to the phasespace properties of the quantum state [28] which characterize the "oscillator" aspect of the physics related to the continuous-variable (CV) nature of the optical field [29,30]; in many cases, it is this latter picture that conceptually connects more directly to the perspectives employed both in classical wave optics and CV quantum information.In the CV approach, it is natural to ask about, for example, the reduced density matrix [31] of a given pulse supermode, which is indispensable for photonic quantum state engineering in the time [32][33][34][35] and/or spectral [36,37] domains, but which is not straightforward to infer from the n-particle correlation functions.Other CV quantities such as the Wigner function and various entanglement measures computed from the reduced density matrix are routinely used to quantify photonic states as resources for quantum information processing [38,39].
In this article, we first review the application of MPS and MPS-based time evolution methods to the problem of quantum pulse propagation.In order to appreciate the optical phase-space dynamics of this novel nonlinear quantum regime, we develop a "demultiplexing" scheme which can be applied to an MPS representation of a quantum pulse to extract the reduced density matrix in a basis of pulse supermodes [33].In our scheme, a carefully chosen sequence of one-and two-mode (i.e., local) linear operations (phaseshifters and beamsplitters) are used to overcome the problem of manipulating nonlocal supermodes-which are not naturally accessible in the MPS framework-by mapping them onto local bins.
As a demonstration, we analyze the quantum propagation of a pulse initialized in a coherent-state Kerr soliton of a 1D χ (3) -nonlinear waveguide [40], using the MPS-based approach of time-evolving block decimation (TEBD) [21,41,42] which is compatible with standard quantum-optical treatments of linear loss [43].We show that the Wigner function of the state contained in the canonical sech-pulse supermode can exhibit a considerable amount of negativity with an enhancement correlated with the peak intensity.We highlight features of the pulse dynamics that indicate departures from established models, such as broadband corrections to the conventional time-dependent Hartree-Fock approximation for fundamental Kerr solitons [44,45], and, for higher-order solitons [46], stark qualitative deviations from classical breather dynamics.Then, we extend our analysis to the quantum propagation of a pulse initialized as a coherentstate simulton [47], a quadratic soliton of a 1D χ (2) nonlinear waveguide composed of two co-propagating pulses at the fundamental harmonic (FH) and second harmonic (SH).We compute a two-mode reduced density matrix for the FH and SH pulse supermodes and reveal their entanglement structure.As a result, we find that non-classicality in the quantum simulton primarily accumulates in a hybrid supermode consisting of a specific linear combination of the FH and SH supermodes.The numerical techniques highlighted in these examples can be generalized to address various questions about quantum pulse propagation expected to arise with the advent of strongly-interacting broadband quantum photonics.

II. MATRIX PRODUCT STATES FOR QUANTUM OPTICAL PULSE PROPAGATION
In this section, we introduce basic concepts of matrix product states (MPS), together with methods to implement time evolutions.While we base our discussions on quantum pulse propagation on a 1D χ (3) -nonlinear waveguide to keep discussions concrete, the basic concepts introduced in this section is general and can be applied to a broader class of systems.We extend our analysis to χ (2) nonlinear waveguides in Sec.V.
The Hamiltonian of a 1D χ (3) -nonlinear waveguide in the moving frame takes a form [20] where photon field annihilation operators fulfill commutation relationships [ φz , φ † z ] = δ(z − z ).The mean-field (i.e., c-number) equation under (1) takes a well-known form of nonlinear Schrödinger equation (NLSE) [40] i∂ where time t and space z have been normalized.In the context 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 analytic studies, it is often easier to work in a discretized coordinate for numerical evaluations.More importantly, discretization of the field allows us to encode system states on an MPS.We consider a finite space interval −L/2 ≤ z ≤ L/2 and discretize it into N spatial bins with size ∆z = L/N .As shown in Fig. 1 (a), we assign a photon annihilation operator âm to the mth the spatial bin for m ∈ {1, 2, . . ., N }.When N is large enough, (1)  Hamiltonian of the form [49,50] Similarly, as shown in the figure, a pulse waveform in the continuous coordinate f (z) is discretized as an Ndimensional vector f = (f 1 , . . ., f N ) .
A generic system state can be written as where Each local Hilbert space is spanned by Fock states |i m (i m ≥ 0).Generally, keeping and updating c i require computational resources that grow exponentially with respect to N .On the other hand, depending on the specific properties of the Hamiltonian, such as locality, much of the states in the entire Hilbert space might not be populated and thus could be excluded.To this end, matrix product states (MPS) allows an effective representation of quantum states given the amount of entanglement is limited and local.
Using an MPS representation with bond dimension χ, (4) can be approximated as [21,41] 1α1 λ [1]  α1 Γ [2]i2 α1α2 λ [2]  α2 where Γ [m] is a rank-3 tensor, λ [m] is a vector, and each α m runs through 1 to χ. Intuitively, ( 5) is a decomposition of a rank-N tensor c i into a product of low-rank tensors.For instance, a coherent pulse with a normalized envelope where all the other tensor elements are zero.Physically, ( 6) is obtained by displacing a supermode Â = m f m âm by α from the vacuum.Notice that parameters needed for (5) has a favorable polynomial scaling of O(χ 2 N ).
The bond dimension χ is related to the maximum amount of entanglement that (5) can support, and larger χ is needed to describe |Ψ with longer-range entanglement.Notably, it is heuristically known that entanglement in 1D quantum many-body systems is often limited [21], making MPS an ideal representation to study quantum pulse propagation.Similarly, a generic operator can also be expressed in the form of a matrix product as where O [m] is a rank-4 tensor.Operators expressed in the form of ( 8) are referred to as matrix product operators (MPO), and their expectation values with respect to MPS are computed via tensor contractions [22].Just in the same way as χ for an MPS is determined by how entangled the state is, bond dimension of an MPO is determined by how non-local the operator Ô is.
In his seminal paper [41], Vidal introduced an algorithm, which is often referred to as time-evolving block decimation (TEBD), to efficiently simulate the time evolution of an MPS.Multiple open-source packages have been developed for TEBD [51,52] to this date.TEBD utilizes the fact that updating an MPS for local onemode or two-mode unitary operations can be done efficiently.While this process involves a truncation of minor singular-value components, truncation error can in principle be arbitrary small by taking large enough bond dimension χ.For instance, as shown in Fig. 1(b), dynamics under (3) is simulated by Trotter decomposing a shorttime evolution Û = e −i Ĥδt into a sequence of one-mode and two-mode operations where Ŝm = exp iδtâ †2 m â2 m /2∆z implements Kerr selfphase modulation (SPM) on the mth bin, and represents hopping interactions between mth and (m + 1)th bins due to quadratic energy dispersions.Note that operations in each parentheses of ( 9) commute, and thus, they can be implemented in parallel numerically.Additionally, higher-order decomposition methods can reduce the Trotter discretization errors [53].For comprehensive overview on time-evolution methods for MPS, we lead readers to Ref. [42].
It is worth mentioning that we can readily include dissipations to the simulation, e.g., by the Monte-Carlo wavefunction method (MCWF) [43].The capability of including loss is particularly important for optical simulations, not only because realistic optical systems often have a non-negligible amount of loss, but also because dissipation plays critical roles in a host of emergent phenomena, such as dissipative Kerr solitons [54] and the physics of PT -symmetric systems [55,56].

III. DEMULTIPLEXING SUPERMODES FROM AN MPS
After a numerical simulation of a quantum pulse propagation using MPS, we wish to calculate various physical properties of the resultant quantum state |Ψ(t) .For instance, the two photon correlation function g (2) ( , m) = â † â † m âm â / â † â â † m âm can be calculated by expressing â † â † m âm â and â † j âj (j = , m) as MPOs and computing their expectation values.This procedure can be extended to compute n-particle correlation function as well.Generically, physical quantities associated to MPOs with larger bond dimension are more expensive to compute.
In the context of quantum engineering and information, it is of great importance to know quantum states populating certain supermodes of interest.To be more precise, let us consider s N supermodes, where the rth pulse supermode [33] is defined as for an arbitrary set of orthonormal vectors {f (1) , . . ., f (s) }.
We denote the Hilbert space for these supermodes as S = S 1 ⊗ S 2 ⊗ . . .S s , where S r is the space for the rth supermode, and the Hilbert space for the rest of the system is denoted as E. Here, we are interested in a reduced density matrix where j = (j 1 , . . ., j s ) , , and similarly for j S .Once we obtain ρS , whose dimension is much smaller than the original Hilbert space, we can study detailed properties of the state, such as Wigner function and entanglement among supermodes with standard techniques for continuous variable systems [29].
To obtain ρS , we ideally would like to have a low-rank MPO representation of the operators where 1E is an identity operator on E. The expectation values of these operators would directly give the matrix elements of ρS via (ρ S ) jj = Ψ|μ jj |Ψ , but the highly non-local structure of μjj makes it nontrivial to find such a low-rank MPO expression.
In the following, as a main result of this research, we describe a procedure to efficiently calculate the reduced density matrix ρS of an MPS |Ψ for an arbitrary set of supermodes comprising S. As shown schematically in Fig. 2, our approach is based on constructing a linear unitary operation V which "demultiplexes" the s (nonlocal) supermodes into the leftmost s (local) spatial bins.In other words, by computing |Φ = V |Ψ , we would like to calculate the matrix elements of ρS via where the operators πjj have a local form where 1m is an identity operator on the mth spatial bin.Explicitly, the MPO representation of πjj is with all other tensor elements zero; because (16) only requires a bond dimension of 1, its expectation value can be efficiently computed.Then the demultiplexing problem posed by ( 14) is to obtain While solutions to (17) are not generally unique, we choose to construct V as a product of local one-mode and two-mode linear operations to allow for the computation of |Φ = V |Ψ using the standard MPS operations R N-2 (1) FIG. 2. Illustration of our supermode demultiplexing scheme.Application of a unitary transformation R(r) , which is composed of one-mode and two-mode gates { R(r) r , . . ., m âm to the rth (local) spatial bin.We show operations up to r = 2, where local annihilation operators of the first and second bins are transformed to Â(1) and Â(2) , respectively.For the purpose of illustration, spatial bins are labeled by associated annihilation operators.
also used in TEBD.Such construction of a linear unitary gate with one-mode and two-mode operations may be seen as a family of general multiport linear interferometers [57,58], where each "port" is a discretized spatial bin of an MPS in our setup.
Our scheme to construct V works in an iterative manner with iteration index r ∈ {1, 2, . . ., s}.In the rth iteration, we construct a transformation R(r) that demultiplexes the rth supermode into the rth spatial bin.Note that, as a consequence, the construction of R(r) is dependent on the partial operations which have already been applied.If the previous (r − 1) iterations were valid, then we can suppose V (r−1) implements the transformations âm → â(r−1) , where â(r−1) Here, the first line of ( 19) means the previous (r − 1) iterations have successfully demultiplexed supermodes 1 through (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−1) has been to only "mix in" components from bins m − r + 1 up to N ; more precisely, for m ≥ r, â(r−1) m is independent of any â such that ≤ m − r. (19) is fulfilled for the base case of V (0) = 1 and c (0) m = δ m at r = 0, so we only need to ensure that our construction of R(r) preserves (19) for any r > 0. At the end of the final sth iteration, we should then find that â(s) m = Â(s) m for 1 ≤ m ≤ s, which satisfies our goal of (17) and allows us to take V = V (s) .
We parametrize R(r) with two sets of angles {θ We then construct R(r) by cascading these one-and two-mode operations according to Because R(r) should demultiplex the rth supermode into the rth spatial bin, we require On the other hand, based on ( 19) and ( 20), we have where g where for m ∈ {r, r + 1, . . ., N }.Notice that the right hand side of (23a) only depends on θ For the first supermode, we have analytic solutions θ (1)  m = arg f (1)  m , ϕ (1) while we generally need to employ numerical methods to demultiplex further supermodes.After solving for all θ m and confirm that ( 19) is fulfilled by the above construction.The iterative procedure culminating in V = V (s) thus demultiplexes all s supermodes into the leftmost s spatial bins.We can then compute the reduced density matrix ρS via (14).

IV. QUANTUM PROPAGATION OF KERR SOLITONS
In this section, we study quantum propagation of a Kerr soliton.A classical soliton solution for (2) with an average photon number of n takes a well-known sech form [40,46] While this solution is exact in the realm of classical optics where it may be thought of as specifying the pulse amplitudes of a coherent state, quantum mechanical effects such as squeezing [60,61], quantum-induced dispersion of pulse envelope [62], and soliton evaporation [63] can become pronounced in the regime of large nonlinearity where n becomes small.Conventionally, quantum mechanical soliton solitons can be treated using the time-dependent Hartree-Fock approximation (TDHF) [44].As long as n 1 holds true, TDHF predicts that a quantum soliton initialized as a coherent state of the envelope φ z (0) approximately evolves according to [45] where α = √ n, and Â is the annihilation operator of the classical soliton-pulse supermode [45].In discretized coordinates, Â = f (sech) m âm , where f m∆z−L/2 (0), up to a proportionality constant independent of m such that f (sech) is normalized.The main feature of the approximate TDHF solution is that it is closed within the subspace S of the soliton supermode f (sech) , so that, e.g., the reduced density matrix of |Ψ(t) in the subspace S has unit purity throughout the dynamics of (26).Nevertheless, due to the Kerr-type nonlinear phase shifts, (26) deviates from a coherent state as it evolves, leading to a variety of interesting phase-space dynamics [16,45,64,65].On the other hand, it is difficult to quantify the accuracy or regime of validity of the TDHF due to its non-perturbative nature, and phasespace dynamics of quantum solitons beyond TDHF in the few-photon regime remain largely unexplored.In the following, we show that our MPS-based scheme serves as a powerful numerical tool to explore this latter regime.
To simulate the propagation of a soliton, we initialize a coherent-state MPS according to (6) with envelope f (sech) and with α = √ n.Using TEBD, we simulate the time evolution of the state under the Hamiltonian (3) with and without linear loss, where the loss dynamics are simulated via MCWF with quantum jump operators { √ κâ 1 , . . ., √ κâ N }, where κ is the power decay rate.Fig. 3(a) shows the time evolution of the photon density distribution (i.e., g (1) correlation function) φ † z φz , for an initial pulse amplitude n = 3.We see that even on the level of g (1) , the pulse envelope exhibits dispersion as a function of time [62], reflecting the fact that the initial classical soliton is not an exact eigenstate of the quantum Hamiltonian; we sometimes refer to such nonclassical dispersion as being "quantum induced".
We next consider the reduced density matrix ρS for the sech-pulse supermode S using the demultiplexing scheme developed 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 , which exhibit a substantial amount of Wigner function negativity and signify non-classicality in the state as might be expected for single-mode Kerr evolution.However, Fig. 3(c) shows that the purity Tr(ρ 2 S ) exhibits a monotonic decay in time, indicating that the sech-pulse subspace S is not closed under the dynamics and in fact, coupling between S and the rest of the system E can act as an effective decoherence channel.Actually, due to the nonlinear nature of the dynamics, ρS may not be pure for any choice of a supermode f in general.Also shown in Fig. 3(c) is the volume of the Wigner function negativity, which serves as a measure of the non-classicality of the state [66].Following an initial increase, the volume of Wigner function negativity starts to decrease after some time, again due to competition between the nonlinear dynamics and the effective decoherence caused by entanglement with E. These features are in stark contrast to the single-mode dynamics predicted by TDHF.
We can also contrast this effective decoherence due to entanglement between S and E with standard dissipation due to linear loss.When loss is incorporated to the simulation, we obtain an ensemble of quantum trajectories {|Ψ 1 (t) , . . ., |Ψ M (t) } via the MCWF method [43].The final reduced density matrix is calculated by averaging the reduced density matrices over the ensemble according to ρS = M −1 M i=1 ρS,i , where ρS,i is the reduced density matrix of |Ψ i .As expected, the linear loss causes We additionally investigate how the amount of excitation affects the phase-space dynamics of pulse propagation.Classically, due to the nonlinear nature of the interactions, the effective nonlinear rate is expected to be enhanced when the peak pulse intensity is larger.In Fig. 3(d), we show the volume of the Wigner function negativity as a function of the average photon number n in the soliton pulse, where we observe that larger pulse excitations increase the rate at which non-classical features are formed, confirming the classical intuition in the few-photon regime.
Finally, to highlight the difference between full numerical simulation and TDHF, Fig. 4 compares the Wigner functions of a soliton pulse initialized with n = 6 obtained by the MPS-based simulation against TDHF (26).While they exhibit qualitatively similar interference patterns, the discrepancies in the time at which a similar phase-space structure is reached indicate that TDHF is overestimating the rate of the nonlinear phase-shift (see t = 0.3 for MPS-based simulation and t = 0.15 for TDHF, for instance).Additionally, full quantum simulation results exhibit visible reduction in the magnitude of Wigner function negativity compared to TDHF, highlighting the effects of aforementioned decoherence due to entanglement between S and E. These discrepancies point to the presence of broadband physics beyond TDHF in soliton propagation.
While the quasi-stationary nature of the quantum dy-  namics of fundamental solitons is in qualitative agreement with their classical behavior, it is in fact possible for highly-quantum photon dynamics to exhibit much more striking deviations from classical behavior.This is the case, for example, if we apply our simulation method to study the quantum propagation of a second-order soliton with classical waveform [46] φ = 2e in 2 t/8 n 3e in 2 t cosh(nz/2) + cosh(3nz/2) 3 cos(n 2 t) + 4 cosh(nz) + cosh(2nz) , which, classically, is a periodic "breather" solution of the NLSE (2).The average photon number of the secondorder soliton is n(2sech) = 4n, where n is the average photon number in its corresponding fundamental soliton.At t = 0, this means the waveform field amplitude of the second-order soliton is twice that of the fundamental soliton, i.e., φ (2sech) z (0) = 2φ (sech) z (0).After t = 0, as shown in Fig. 5(a), the classical waveform of the second-order soliton exhibits significant narrowing and a characteristic triplet structure.On the other hand, as shown in Fig. 5(b), full quantum evolution of a second-order soliton instantiated in a few-photon coherent state of ( 27) at t = 0 exhibits qualitatively different dynamics.While the photon density distribution exhibits some narrowing of the pulse width initially, the peak pulse intensity fails to reach the level expected from the classical solution.Moreover, no signature of the triplet structure is observed.We attribute these features to the quantuminduced dispersion of the pulse envelope that we also observed in the quantum propagation of fundamental solitons [62], which appears to play a more critical role in the evolution of this higher-order soliton.

V. QUANTUM PROPAGATION OF SIMULTONS
In this section, we apply our technique to a 1D χ (2)  nonlinear waveguide in which interactions occur between a fundamental harmonic (FH) band and a second harmonic (SH) band.After normalization with respect to time and space, the system Hamiltonian takes the form [67,68] where φz and ψz are respectively FH and SH local field annihilation operators with commutation relationships [ φz We have assumed that the FH and SH carriers are group-velocity matched, while β represents the group velocity dispersion of SH relative to FH. (28) with where âm and bm are the FH and SH field annihilation operators for the mth spatial bin, respectively.The presence of both FH and SH fields requires some care in applying the two-mode operations implementing TEBD for the Hamiltonian (29).One approach is shown schematically in Fig. 6, where we prepare 2N MPS bins to represent N spatial bins, with FH and SH modes encoded in an alternating manner.To apply Trotterization as in Sec.II, a short-time unitary evolution e −i Ĥδt is decomposed into two-mode operations Ûλ,m = e −i Ĥλ,m δt (λ = a, b, NL) and applied as shown in Fig. 6.Since âm and âm+1 are no longer next to each other in this representation, we also utilize a swap operation ÛSWAP [41] to bring these modes together in an alternating manner for the application of Ûa,m (and similarly for Ûb,m ).
The Hamiltonian (28) supports various classical soliton solutions [69,70].Specifically, an analytic "simulton" solution exists for β = 2 with the form [47] φ z (t) = φ 0 sech 2 φ 0 /6 z e iφ0t/3 (30a) where φ 0 is related to the average FH photon number n via φ 0 = 3 3n 2 /32.As was done for (25), we discretize and normalize (30) at t = 0 to construct FH and SH supermodes with annhilation operators denoted Â and B, respectively.The initial simulton MPS state is that of a coherent state in Â and B with displacements √ n and − √ n/2, respectively.In Fig. 7(a), we show the time evolution via TEBD of the photon density distribution of a pulse initialized in the classical simulton supermode.As for the case of the Kerr soliton, the dynamics show a quantum-induced dispersion of the pulse envelope not predicted by classical dynamics; as part of this process, we also numerically observe a slight exchange of excitations between FH and SH.
To investigate the joint phase-space properties of the FH and SH pulse supermodes, we calculate two-mode reduced density matrix ρS where S = A ⊗ B, the joint Hilbert 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 the negativity measure N (ρ) (not to be confused with the Wigner function negativity used earlier) [71], defined for a bipartite density matrix ρ as where T A is partial transposition with respect to the first mode, and • 1 is the trace norm.Here, N (ρ) > 0 serves as a sufficient condition for the entanglement.Generally, we can change the value of the N by considering hybrid mixtures of the modes A and B. Finding the linear combination that best "disentangles" A and B therefore provides insight into their entanglement structure.More specifically, we introduce a two-mode linear operation Ŵ (Φ, Θ) = exp Φ e iΘ Â † B − e −iΘ Â B † with −π/4 ≤ Φ ≤ π/4 and −π/2 ≤ Θ ≤ π/2, which implements the transformation ρ S = W † ρS Ŵ .Importantly, the transformation Ŵ0 = Ŵ (Φ 0 , Θ 0 ) that minimizes N (ρ S ) is expected to maximize the amount of information available in the single-mode reduced density matrices ρ A = Tr B (ρ S ) and ρ B = Tr A (ρ S ).
In Fig. 7(b), for the final state of the pulse propagation at t = 4, we map the negativity N (ρ S ) for various (Φ, Θ), which shows a clear minimum at the marked position of (Φ 0 , Θ 0 ).In Fig. 7(c), we show the Wigner functions of the single-mode reduced density matrices before and after the transformation Ŵ0 .Before applying the transformation, the Wigner functions of both ρA and ρB are crescent-shaped with no negativity, indicating that they are highly mixed states due to the entanglement between FH and SH.On the other hand, after application of Ŵ0 , the Wigner function of ρ A remarkably exhibits considerable non-classicality, while the Wigner function of ρ B resembles that of a coherent state.This reveals a somewhat surprising feature of the simulton quantum dynamics: a hybrid supermode Â0 = cos Φ 0 Â + e iΘ0 sin Φ 0 B, composed of both FH and SH components, is the one which predominantly experiences strongly nonlinear dynamics, while the other hybrid supermode B0 = cos Φ 0 B − e −iΘ0 sin Φ 0 Â experiences little nonlinearity.A similar analysis can, in principle, be applied to a broader class of solitons and general pulse propagation [46].

VI. CONCLUSION
In this research, we have motivated the use of MPS techniques to efficiently represent and simulate a quantum optical pulse as it dynamically propagates through a nonlinear 1D waveguide.In doing so, we have developed a numerical method to overcome the problem of efficiently accessing and manipulating nonlocal pulse supermodes of the local MPS representation, allowing us to view for the first time the full quantum dynamics of the pulse in a phase-space picture.As a demonstration, we have performed quantum simulations of Kerr soliton propagation and observed that the phase-space portraits of an initially classical sech-pulse supermode can evolve highly non-classical features, i.e., Wigner function negativity.These results have been contrasted with predictions based on TDHF, highlighting the presence of rich quantum dynamics beyond conventional approximations for quantum Kerr solitons.We have also extended our analysis to the quantum propagation of a χ (2) simulton and have revealed unexpected entanglement structure between the FH and SH pulses of the simulton, identifying a hybrid supermode that predominantly exhibits nonclassical features.Our scheme is compatible with local dissipation associated with, e.g., waveguide losses, and, more generally, could be applied to any one-dimensional photonic system in principle.Considering the rapid recent progress towards single-photon nonlinearities in dispersion-engineered and highly nonlinear nanophotonic platforms, it is of imminent interest to establish a unified theoretical framework in which to understand the quantum dynamics of photons in such devices.To this end, our work takes a step towards bridging the significant conceptual gaps between classical wave optics, CV photonic quantum information, and strongly-interacting quantum many-body physics, all of which are expected to play important roles in conceptualizing and engineering the future of broadband quantum optics.

N − 1 }
(i.e., the rth iteration requires 2(N − r) + 1 parameters).For m < N , these parameters are taken to define a set of phaseshifter/beamsplitter operations and(22), we can solve for the angles via

m
with m < m, and thus, we can solve the equations starting from m = r towards n iteratively to straightforwardly determine all of θ

FIG. 5 .
FIG. 5. Time evolution of the photon density distribution of a pulse instantiated as a coherent-state second-order soliton with mean photon number n(2sech) = 8.(a) Classically expected dynamics given by (27).(b) Full quantum simulation under the Hamiltonian (1) using MPS with bond dimension χ = 60.
can be approximated by a discretized Bose-Hubbard FIG.6.An example implementation of a short-time evolution under the χ(2)Hamiltonian (28) using TEBD.Discretized FH modes {â1, . . ., âN } and SH modes { b1, . . ., bN } are encoded on an MPS in an alternating manner.Swap operations ÛSWAP are used to bring distant modes together.Two-mode operations Ûa,m, Ûb,m , and ÛNL,m are for FH dispersion, SH dispersion, and nonlinear three-wave mixing interaction, respectively.