## Abstract

Ultrashort 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ïvely 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. 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 harmonics. 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.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. 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–6]. Recent efforts in the field of nonlinear nanophotonics have resulted in the development of ultralow-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–13]. In particular, advances in dispersion engineering on these platforms enable ultrashort-pulse operation [14,15], where the available peak power further leverages the material nonlinearities by orders of magnitude, bringing the possibility of engineering highly non-classical 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 such as 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, ${\lt}10$ interacting photons in each bin, we would naïvely need to compute the evolution of a ${\sim}{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 one-dimensionality 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 for the analysis of optical systems, such as time-multiplexed quantum optics [26] and pulse propagation through mesoscopic atomic clouds [27] and atom-coupled chiral waveguides [28]. 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 phase-space properties of the quantum state [29], which characterize the “oscillator” aspect of the physics related to the continuous-variable (CV) nature of the optical field [30,31]; in many cases, it is this latter picture that conceptually connects more directly to the perspectives employed in both classical wave optics and CV quantum information. In the CV approach, it is natural to ask about, for example, the reduced density matrix [32] of a given pulse supermode, which is indispensable for photonic quantum state engineering in time [33–36] and/or spectral [37,38] 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 [39,40].

In this paper, we first review the application of MPS and MPS-based time evolution methods to the problem of quantum pulse propagation. To appreciate the optical phase-space dynamics of this novel nonlinear quantum regime, we develop a “demultiplexing” scheme that can be applied to an MPS representation of a quantum pulse to extract the reduced density matrix in a basis of pulse supermodes [34]. In our scheme, a carefully chosen sequence of one- and two-mode (i.e., local) linear operations (phase shifters and beam splitters) 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 ${\chi ^{(3)}}$-nonlinear waveguide [41], using the MPS-based approach of time-evolving block decimation (TEBD) [21,42,43], which is compatible with standard quantum optical treatments of linear loss [44]. 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 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 approximation (TDH) for fundamental Kerr solitons [45,46] and, for higher-order solitons [47], stark qualitative deviations from classical breather dynamics. Then, we extend our analysis to the quantum propagation of a pulse initialized as a coherent-state simulton [48], a quadratic soliton of a 1D ${\chi ^{(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 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.

## 2. MATRIX PRODUCT STATES FOR QUANTUM OPTICAL PULSE PROPAGATION

In this section, we introduce basic concepts of MPS, together with methods to implement time evolutions. While we base our discussions on quantum pulse propagation on a 1D ${\chi ^{(3)}}$-nonlinear waveguide to keep discussions concrete, the basic concepts introduced in this section are general and can be applied to a broader class of systems. We extend our analysis to ${\chi ^{(2)}}$-nonlinear waveguides in Section 5.

The Hamiltonian of a 1D ${\chi ^{(3)}}$-nonlinear waveguide in the moving frame takes the form [20]

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 \le z \le L/2$ and discretize it into $N$ spatial bins with size $\Delta z = L/N$. As shown in Fig. 1 (a), we assign a photon annihilation operator ${\hat a_m}$ to the $m$th spatial bin for $m \in \{1,2, \ldots ,N\}$. When $N$ is large enough, Eq. (1) can be approximated by a discretized Bose–Hubbard Hamiltonian of the form [50,51]

Similarly, as shown in the figure, a pulse waveform in the continuous coordinate $f(z)$ is discretized as an $N$-dimensional vector ${\boldsymbol f}=(f_1,\ldots,f_N)^\top$.

A generic system state can be written as

where ${{\boldsymbol i}}= (i_1,i_2,\ldots,i_N)^\top$, and $|{\boldsymbol i}\rangle = |{i_1}\rangle \otimes |{i_2}\rangle \cdots |{i_N}\rangle$. Each local Hilbert space is spanned by Fock states $|{i_m}\rangle ({i_m} \ge 0)$. Generally, keeping and updating ${c_{\boldsymbol 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, many of the states in the entire Hilbert space might not be populated and thus could be excluded. To this end, MPS allows an effective representation of quantum states given that the amount of entanglement is limited and local.Using an MPS representation with bond dimension $\chi$, Eq. (4) can be approximated as [21,42]

Similarly, a generic operator,

In his seminal paper [42], Vidal introduced an algorithm, often referred to as TEBD, to efficiently simulate the time evolution of an MPS. Multiple open-source packages have been developed for TEBD [52,53] to this date. TEBD utilizes the fact that updating an MPS for local one-mode 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 arbitrarily small by taking large enough bond dimension $\chi$. For instance, as shown in Fig. 1(b), dynamics under Eq. (3) are simulated by Trotter decomposing a short-time evolution $\hat U = {e^{- {i}\hat H\delta t}}$ into a sequence of one-mode and two-mode operations:

In numerical implementations of TEBD, we represent quantum states and operators in the Fock basis with finite photon number cutoff, chosen to be large enough to accommodate the local photon population $\approx \langle \hat a_m^\dagger {\hat a_m}\rangle$. It is important to note that photons comprising the pulse are distributed over a large number of spatial bins, and thus, the photon number cutoff for each spatial bin can be much smaller than the total photon number in the system.

It is worth mentioning that we can readily include dissipations to the simulation, e.g., by the Monte Carlo wavefunction method (MCWF) [44]. 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 [55] and the physics of parity-time (${\cal P}{\cal T}$) symmetric systems [56,57].

## 3. DEMULTIPLEXING SUPERMODES FROM AN MPS

After a numerical simulation of a quantum pulse propagation using MPS, we calculate various physical properties of the resultant quantum state $|\Psi (t)\rangle$. For instance, the two-photon correlation function ${g^{(2)}}(\ell ,m) = \langle \hat a_\ell ^\dagger \hat a_m^\dagger {\hat a_m}{\hat a_\ell}\rangle /\langle \hat a_\ell ^\dagger {\hat a_\ell}\rangle \langle \hat a_m^\dagger {\hat a_m}\rangle$ can be calculated by expressing $\hat a_\ell ^\dagger \hat a_m^\dagger {\hat a_m}{\hat a_\ell}$ and $\hat a_j^\dagger {\hat a_j} (j = \ell ,m)$ as MPOs and computing their expectation values. This procedure can be extended to compute the $n$-particle correlation function as well. Generically, physical quantities associated with MPOs with larger bond dimensions 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 \ll N$ supermodes, where the $r$th pulse supermode [34] is defined as

for an arbitrary set of orthonormal vectors $\{{{\boldsymbol f}^{(1)}}, \ldots ,{{\boldsymbol f}^{(s)}}\}$. We denote the Hilbert space for these supermodes as ${\cal S} = {{\cal S}_1} \otimes {{\cal S}_2} \otimes \ldots {{\cal S}_s}$, where ${{\cal S}_r}$ is the space for the $r$th supermode, and the Hilbert space for the rest of the system is denoted as ${\cal E}$. Here, we are interested in a reduced density matrix:To obtain ${\hat \rho _{\cal S}}$, we ideally would like to have a low-rank MPO representation of the operators

In the following, as a main result of this research, we describe a procedure to efficiently calculate the reduced density matrix ${\hat \rho _{\cal S}}$ of an MPS $|\Psi \rangle$ for an arbitrary set of supermodes comprising ${\cal S}$. As shown schematically in Fig. 2, our approach is based on constructing a linear unitary operation $\hat V$ that “demultiplexes” the $s$ (nonlocal) supermodes into the leftmost $s$ (local) spatial bins. In other words, by computing $|\Phi \rangle = \hat V|\Psi \rangle$, we would like to calculate the matrix elements of ${\hat \rho _{\cal S}}$ via

While solutions to Eq. (17) are not generally unique, we choose to construct $\hat V$ as a product of local one-mode and two-mode linear operations to allow for the computation of $|\Phi \rangle = \hat V|\Psi \rangle$ using the standard MPS operations 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 [58,59], where each “port” is a discretized spatial bin of an MPS in our setup.

Our scheme to construct $\hat V$ works in an iterative manner with iteration index $r \in \{1,2, \ldots ,s\}$. In the $r$th iteration, we construct a transformation ${\hat R^{(r)}}$ that demultiplexes the $r$th supermode into the $r$th spatial bin. Note that, as a consequence, the construction of ${\hat 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 ${\hat V^{(r - 1)}}$ implements the transformations ${\hat a_m} \mapsto \hat a_m^{(r - 1)} = {\hat V^{(r - 1)\dagger}}{\hat a_m}{\hat V^{(r - 1)}}$, whereHere, the first line of Eq. (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 \ge r$, the effect of the transformation ${\hat V^{(r - 1)}}$ has been to “mix in” only components from bins $m - r + 1$ up to $N$; more precisely, for $m \ge r$, $\hat a_m^{(r - 1)}$ is independent of any ${\hat a_\ell}$ such that $\ell \le m - r$. Equation (19) is fulfilled for the base case of ${\hat V^{(0)}} = \hat {\unicode{x1D7D9}}$ and $c_{{m\ell}}^{(0)} = {\delta _{{m\ell}}}$ at $r = 0$, so we need to ensure only that our construction of ${\hat R^{(r)}}$ preserves Eq. (19) for any $r \gt 0$. At the end of the final $s$th iteration, we should then find that $\hat a_m^{(s)} = \hat A_m^{(s)}$ for $1 \le m \le s$, which satisfies our goal of Eq. (17) and allows us to take $\hat V = {\hat V^{(s)}}$.

We parametrize ${\hat R^{(r)}}$ with two sets of angles $\{\theta _r^{(r)}, \ldots ,\theta _N^{(r)}\}$ and $\{\varphi _r^{(r)}, \ldots ,\varphi _{N - 1}^{(r)}\}$ (i.e., the $r$th iteration requires $2(N - r) + 1$ parameters). For $m \lt N$, these parameters are taken to define a set of phase shifter/beam splitter operations:

For $m = N$, we fix $\hat R_N^{(r)} = {e^{{ i}\theta _N^{(r)}\hat a_N^\dagger {{\hat a}_N}}}$. We then construct ${\hat R^{(r)}}$ by cascading these one- and two-mode operations according to

Because ${\hat R^{(r)}}$ should demultiplex the $r$th supermode into the $r$th spatial bin, we require

On the other hand, based on Eqs. (19) and (20), we have

For the first supermode, we have analytic solutions

## 4. QUANTUM PROPAGATION OF KERR SOLITONS

In this section, we study quantum propagation of a Kerr soliton. A classical soliton solution to Eq. (2) with an average photon number of $\bar n$ takes a well-known sech form [41,47]:

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 [61,62], quantum-induced dispersion of pulse envelope [63], and soliton evaporation [64] can become pronounced in the regime of large nonlinearity where $\bar n$ becomes small.

Conventionally, quantum mechanical soliton solutions can be treated using the TDH [45]. As long as $\bar n \gg 1$ holds true, TDH predicts that a quantum soliton initialized as a coherent state of the envelope ${\phi _z}(0)$ approximately evolves according to [46]

To simulate the propagation of a soliton, we initialize a coherent-state MPS according to Eq. (6) with envelope ${{\boldsymbol f}^{\text{(sech)}}}$ and with $\alpha = \sqrt {\bar n}$. Using TEBD, we simulate the time evolution of the state under the Hamiltonian Eq. (3) with and without linear loss, where the loss dynamics are simulated via MCWF with quantum jump operators $\{\sqrt \kappa {\hat a_1}, \ldots ,\sqrt \kappa {\hat a_N}\}$, where $\kappa$ is the power decay rate. Figure 3(a) shows the time evolution of the photon density distribution [i.e., ${g^{(1)}}$ correlation function] $\langle \hat \phi _z^\dagger {\hat \phi _z}\rangle$, for an initial pulse amplitude $\bar n = 3$. We see that even on the level of ${g^{(1)}}$, the pulse envelope exhibits dispersion as a function of time [63], reflecting the fact that the initial classical soliton is not an exact eigenstate of the quantum Hamiltonian; we sometimes refer to such non-classical dispersion as being “quantum induced.”

We next consider the reduced density matrix ${\hat \rho _{\cal S}}$ for the sech-pulse supermode ${\cal S}$ using the demultiplexing scheme developed in Section 3 with the envelope function ${{\boldsymbol f}^{\text{(sech)}}}$. Consider first the case without loss, or $\kappa = 0$. Figure 3(b) shows snapshots of the Wigner functions of ${\hat \rho _{\cal 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 $\text{Tr}(\hat \rho _{\cal S}^2)$ exhibits a monotonic decay in time, indicating that the sech-pulse subspace ${\cal S}$ is not closed under the dynamics and in fact, coupling between ${\cal S}$ and the rest of the system ${\cal E}$ can act as an effective decoherence channel. Actually, due to the nonlinear nature of the dynamics, ${\hat \rho _{\cal S}}$ may not be pure for any choice of a supermode ${\boldsymbol 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 [67]. 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 ${\cal E}$. These features are in stark contrast to the single-mode dynamics predicted by TDH.

We can also contrast this effective decoherence due to entanglement between ${\cal S}$ and ${\cal E}$ with standard dissipation due to linear loss. When loss is incorporated into the simulation, we obtain an ensemble of quantum trajectories $\{|{\Psi _1}(t)\rangle , \ldots ,|{\Psi _M}(t)\rangle \}$ via the MCWF method [44]. The final reduced density matrix is calculated by averaging the reduced density matrices over the ensemble according to ${\hat \rho _{\cal S}} = {M^{- 1}}\sum\nolimits_{i = 1}^M {\hat \rho _{{\cal S},i}}$, where ${\hat \rho _{{\cal S},i}}$ is the reduced density matrix of $|{\Psi _i}\rangle$. As expected, the linear loss causes a decay in the amplitude of the photon density distribution as shown in Fig. 3(a), while quantum mechanically, Figs. 3(b) and 3(c) show that the non-classical features of the state are critically diminished by the presence of the linear loss.

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 Wigner function negativity as a function of the average photon number $\bar 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 TDH, Fig. 4 compares the Wigner functions of a soliton pulse initialized with $\bar n = 6$ obtained by the MPS-based simulation against TDH Eq. (26). While they exhibit qualitatively similar interference patterns, the discrepancies in the time at which a similar phase-space structure is reached indicate that TDH is overestimating the rate of the nonlinear phase shift (see $t = 0.3$ for MPS-based simulation and $t = 0.15$ for TDH, for instance). Additionally, full quantum simulation results exhibit visible reduction in the magnitude of Wigner function negativity compared to TDH, highlighting the effects of the aforementioned decoherence due to entanglement between ${\cal S}$ and ${\cal E}$. These discrepancies point to the presence of broadband physics beyond TDH in soliton propagation.

Aside from TDH, treatments based on linearization of quantum noise around the mean field have been used to concisely capture quantum dynamics of solitons [61,68,69], when the quantum fluctuations can be well approximated as Gaussian multimode squeezing around the mean field [60]. These linearized approaches exclude non-Gaussian features, such as Wigner function negativities, that can arise in the strongly nonlinear regime. In Fig. 4, we also superimpose the shape of the Gaussian phase-space distribution predicted by the linearized formalism. We observe reasonable agreement with MPS simulations at early times, while they differ significantly as Kerr SPM induces crescent-like structures with large non-Gaussianity in phase space.

While the quasi-stationary nature of the quantum dynamics 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 [47]:

## 5. QUANTUM PROPAGATION OF SIMULTONS

In this section, we apply our technique to a 1D ${\chi ^{(2)}}$-nonlinear waveguide in which interactions occur between an FH band and an SH band. After normalization with respect to time and space, the system Hamiltonian takes the form [70,71]

The presence of both FH and SH fields requires some care in applying the two-mode operations for implementing TEBD for the Hamiltonian Eq. (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 Section 2, a short-time unitary evolution ${e^{- { i}\hat H\delta t}}$ is decomposed into two-mode operations ${\hat U_{\lambda ,m}} = {e^{- { i}{{\hat H}_{\lambda ,m}}\delta t}} (\lambda = \text{a},\text{b},\text{NL})$ and applied as shown in Fig. 6. Since ${\hat a_m}$ and ${\hat a_{m + 1}}$ are no longer next to each other in this representation, we also utilize a swap operation ${\hat U_{\text{SWAP}}}$ [42] to bring these modes together in an alternating manner for the application of ${\hat U_{\text{a},m}}$ (and similarly for ${\hat U_{\text{b},m}}$).

The Hamiltonian (28) supports various classical soliton solutions [72,73]. Specifically, an analytic “simulton” solution exists for $\beta = 2$ with the form [48]

To investigate the joint phase-space properties of the FH and SH pulse supermodes, we calculate two-mode reduced density matrix ${\hat \rho _{\cal S}}$, where ${\cal S} = {\cal A} \otimes {\cal B}$, the joint Hilbert space of FH and SH, respectively. In general, the state of the system features entanglement between ${\cal A}$ and ${\cal B}$; to quantify this entanglement, we utilize the negativity measure ${\cal N}(\hat \rho)$ (not to be confused with the Wigner function negativity used earlier) [74], defined for a bipartite density matrix $\hat \rho$ as

Generally, we can change the value of the ${\cal N}$ by considering hybrid mixtures of modes ${\cal A}$ and ${\cal B}$. Finding the linear combination that best “disentangles” ${\cal A}$ and ${\cal B}$ therefore provides insight into their entanglement structure. More specifically, we introduce a two-mode linear operation $\hat W(\Phi ,\Theta) = \exp \{\Phi ({e^{{ i}\Theta}}{{\hat A}^\dagger}\hat B - {e^{- \text{i}\Theta}}\hat A{{\hat B}^\dagger}) \}$, with ${-}\pi /4 \le \Phi \le \pi /4$ and ${-}\pi /2 \le \Theta \le \pi /2$, which implements the transformation ${\hat \rho _{\cal S}^\prime} = {W^\dagger}{\hat \rho _{\cal S}}\hat W$. Importantly, the transformation ${\hat W_0} = \hat W({\Phi _0},{\Theta _0})$ that minimizes ${\cal N}({\hat \rho _{\cal S}}^\prime)$ is expected to maximize the amount of information available in the single-mode reduced density matrices ${\hat \rho _{\cal A}}^\prime = {\text{Tr}_{\cal B}}({\hat \rho _{\cal S}^\prime})$ and ${\hat \rho _{\cal B}^\prime} = {\text{Tr}_{\cal A}}({\hat \rho _{\cal S}}^\prime)$.

In Fig. 7(b), for the final state of pulse propagation at $t = 4$, we map the negativity ${\cal N}({\hat \rho _{\cal S}}^\prime)$ for various $(\Phi ,\Theta)$, which shows a clear minimum at the marked position of $({\Phi _0},{\Theta _0})$. In Fig. 7(c), we show the Wigner functions of the single-mode reduced density matrices before and after the transformation ${\hat W_0}$. Before applying the transformation, the Wigner functions of both ${\hat \rho _{\cal A}}$ and ${\hat \rho _{\cal 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 ${\hat W_0}$, the Wigner function of ${\hat \rho _{\cal A}}^\prime $ remarkably exhibits considerable non-classicality, while the Wigner function of ${\hat \rho _{\cal B}}^\prime $ resembles that of a coherent state. This reveals a somewhat surprising feature of the simulton quantum dynamics: a hybrid supermode ${\hat A_0} = \cos {\Phi _0}\hat A + {e^{{i}{\Theta _0}}}\sin {\Phi _0}\hat B$, composed of both FH and SH components, is the one that predominantly experiences strongly nonlinear dynamics, while the other hybrid supermode ${\hat B_0} = \cos {\Phi _0}\hat B - {e^{- {i}{\Theta _0}}}\sin {\Phi _0}\hat A$ experiences little nonlinearity. A similar analysis can, in principle, be applied to a broader class of solitons and general pulse propagation [47].

## 6. 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 TDH, 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 ${\chi ^{(2)}}$ simulton and have revealed unexpected entanglement structure between FH and SH pulses of the simulton, identifying a hybrid supermode that predominantly exhibits non-classical 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.

## Funding

Army Research Office (W911NF-16-1-0086); National Science Foundation (CCF-1918549, PHY-2011363).

## Acknowledgment

R. Y. developed the numerical techniques, performed the simulations, and generated the figures. E. N. and H. M. advised and directed the project. R. Y. and E. N. wrote the paper with detailed input and feedback from all authors. All authors contributed significantly to the conception of the project. The authors thank NTT Research for financial and technical support. R. Y. thanks Tomohiro Soejima for helpful discussions. R. Y. is supported by a Stanford Q-FARM Ph.D. Fellowship and the Masason Foundation.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## Supplemental document

See Supplement 1 for supporting content.

## REFERENCES

**1. **The LIGO Scientific Collaboration, “Enhanced sensitivity of the LIGO gravitational wave detector by using squeezed states of light,” Nat. Photonics **7**, 613–619 (2013). [CrossRef]

**2. **M. Tsang, R. Nair, and X.-M. Lu, “Quantum theory of superresolution for two incoherent optical point sources,” Phys. Rev. X **6**, 031033 (2016). [CrossRef]

**3. **N. Gisin and R. Thew, “Quantum communication,” Nat. Photonics **1**, 165–171 (2007). [CrossRef]

**4. **J. L. O’Brien, A. Furusawa, and J. Vučković, “Photonic quantum technologies,” Nat. Photonics **3**, 687–695 (2009). [CrossRef]

**5. **D. E. Chang, V. Vuletić, and M. D. Lukin, “Quantum nonlinear optics - photon by photon,” Nat. Photonics **8**, 685–694 (2014). [CrossRef]

**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 using photons,” Science **370**, 1460–1463 (2020). [CrossRef]

**7. **M. Zhang, C. Wang, R. Cheng, A. Shams-Ansari, and M. Lončar, “Monolithic ultra-high-Q lithium niobate microring resonator,” Optica **4**, 1536–1537 (2017). [CrossRef]

**8. **C. Wang, C. Langrock, A. Marandi, M. Jankowski, M. Zhang, B. Desiatov, M. M. Fejer, and M. Lončar, “Ultrahigh-efficiency wavelength conversion in nanophotonic periodically poled lithium niobate waveguides,” Optica **5**, 1438–1441 (2018). [CrossRef]

**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 **7**, 1654–1659 (2020). [CrossRef]

**10. **M. Placke and S. Ramelow, “Engineering AlGaAs-on-insulator towards quantum optical applications,” Opt. Lett. **45**, 6763–6766 (2020). [CrossRef]

**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_{3}N_{4} ring resonator,” Phys. Rev. Lett. **122**, 153906 (2019). [CrossRef]

**12. **M. Heuck, K. Jacobs, and D. R. Englund, “Photon-photon interactions in dynamically coupled cavities,” Phys. Rev. A **101**, 042322 (2020). [CrossRef]

**13. **A. W. Bruch, X. Liu, J. B. Surya, C.-L. Zou, and H. X. Tang, “On-chip χ^{(2)} microring optical parametric oscillator,” Optica **6**, 1361–1366 (2019). [CrossRef]

**14. **M. Jankowski, C. Langrock, B. Desiatov, A. Marandi, C. Wang, M. Zhang, C. R. Phillips, M. Lončar, and M. M. Fejer, “Ultrabroadband nonlinear optics in nanophotonic periodically poled lithium niobate waveguides,” Optica **7**, 40–46 (2020). [CrossRef]

**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-chip octave-spanning supercontinuum generation,” Opt. Express **20**, 1685–1690 (2012). [CrossRef]

**16. **R. Yanagimoto, T. Onodera, E. Ng, L. G. Wright, P. L. McMahon, and H. Mabuchi, “Engineering a Kerr-based deterministic cubic phase gate via Gaussian operations,” Phys. Rev. Lett. **124**, 240503 (2020). [CrossRef]

**17. **K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E. Northup, and H. J. Kimble, “Photon blockade in an optical cavity with one trapped atom,” Nature **436**, 87–90 (2005). [CrossRef]

**18. **J. Javanainen and J. Ruostekoski, “Light propagation beyond the mean-field theory of standard optics,” Opt. Express **24**, 993–1001 (2016). [CrossRef]

**19. **R. Yanagimoto, E. Ng, M. P. Jankowski, T. Onodera, M. M. Fejer, and H. Mabuchi, “Broadband parametric downconversion as a discrete-continuum Fano interaction,” arXiv:2009.01457 (2020).

**20. **P. D. Drummond and M. Hillery, *The Quantum Theory of Nonlinear Optics* (Cambridge University, 2014).

**21. **G. Vidal, “Efficient simulation of one-dimensional quantum many-body systems,” Phys. Rev. Lett. **93**, 040502 (2004). [CrossRef]

**22. **U. Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Ann. Phys. **326**, 96–192 (2011). [CrossRef]

**23. **R. Orús, “A practical introduction to tensor networks: matrix product states and projected entangled pair states,” Ann. Phys. **349**, 117–158 (2014). [CrossRef]

**24. **D. Muth and M. Fleischhauer, “Dynamics of pair correlations in the attractive Lieb-Liniger gas dominik,” Phys. Rev. Lett. **105**, 150403 (2010). [CrossRef]

**25. **A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, “Measuring entanglement growth in quench dynamics of Bosons in an optical lattice,” Phys. Rev. Lett. **109**, 020505 (2012). [CrossRef]

**26. **M. Lubasch, A. A. Valido, J. J. Renema, W. S. Kolthammer, D. Jaksch, M. S. Kim, I. Walmsley, and R. García-Patrón, “Tensor network states in time-bin quantum optics,” Phys. Rev. A **97**, 062304 (2018). [CrossRef]

**27. **M. T. Manzoni, D. E. Chang, and J. S. Douglas, “Simulating quantum light propagation through atomic ensembles using matrix product states,” Nat. Commun. **8**, 1743 (2017). [CrossRef]

**28. **S. Mahmoodian, G. Calajó, D. E. Chang, K. Hammerer, and A. S. Sørensen, “Dynamics of many-body photon bound states in chiral waveguide QED,” Phys. Rev. X **10**, 031011 (2020). [CrossRef]

**29. **K. E. Cahill and R. J. Glauber, “Density operators and quasi probability distributions,” Phys. Rev. **177**, 1882 (1969). [CrossRef]

**30. **S. L. Braunstein and P. van Loock, “Quantum information with continuous variables,” Rev. Mod. Phys. **77**, 513 (2005). [CrossRef]

**31. **D. Gottesman, A. Kitaev, and J. Preskill, “Encoding a qubit in an oscillator,” Phys. Rev. A **64**, 012310 (2001). [CrossRef]

**32. **M. A. Nielsen and I. L. Chuang, *Quantum Computation and Quantum Information* (Cambridge University, 2000).

**33. **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 single spatial mode,” Phys. Rev. Lett. **111**, 150501 (2013). [CrossRef]

**34. **B. Brecht, D. V. Reddy, C. Silberhorn, and M. G. Raymer, “Photon temporal modes: a complete framework for quantum information science,” Phys. Rev. X **5**, 041017 (2015). [CrossRef]

**35. **W. Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoensombutamon, H. Emura, R. N. Alexander, S. Takeda, J. Yoshikawa, N. C. Menicucci, H. Yonezawa, and A. Furusawa, “Generation of time-domain-multiplexed two-dimensional cluster state,” Science **366**, 373–376 (2019). [CrossRef]

**36. **V. Ansari, J. M. Donohue, B. Brecht, and C. Silberhorn, “Tailoring nonlinear processes for quantum optics with pulsed temporal-mode encodings,” Optica **5**, 534–550 (2018). [CrossRef]

**37. **J. M. Lukens and P. Lougovski, “Frequency-encoded photonic qubits for scalable quantum information processing,” Optica **4**, 8–16 (2017). [CrossRef]

**38. **J. Roslund, R. M. de Araújo, S. Jiang, C. Fabre, and N. Treps, “Wavelength-multiplexed quantum networks with ultrafast frequency combs,” Nat. Photonics **8**, 109–112 (2014). [CrossRef]

**39. **E. Chitambar and G. Gour, “Quantum resource theories,” Rev. Mod. Phys. **91**, 025001 (2019). [CrossRef]

**40. **F. Albarelli, M. G. Genoni, M. G. A. Paris, and A. Ferraro, “Resource theory of quantum non-Gaussianity and Wigner negativity,” Phys. Rev. A **98**, 052350 (2018). [CrossRef]

**41. **G. P. Agrawal, *Nonlinear Fiber Optics*, 6th ed. (Academic, 2019).

**42. **G. Vidal, “Efficient classical simulation of slightly entangled quantum computers,” Phys. Rev. Lett. **91**, 147902 (2003). [CrossRef]

**43. **J. J. García-Ripoll, “Time evolution of matrix product states,” New J. Phys. **8**, 305 (2006). [CrossRef]

**44. **H. W. Wiseman and G. J. Milburn, *Quantum Measurement and Control* (Cambridge University, 2009).

**45. **A. Haus and Y. Lai, “Quantum theory of solitons in optical fibers. I. Time-dependent Hartree approximation,” Phys. Rev. A **40**, 5729 (1989). [CrossRef]

**46. **E. W. Wright, “Quantum theory of soliton propagation in an optical fiber using the Hartree approximation,” Phys. Rev. A **43**, 3836 (1991). [CrossRef]

**47. **Y. S. Kivshar and G. P. Agrawal, *Optical Solitons* (Academic, 2003).

**48. **M. J. Werner and P. D. Drummond, “Simulton solutions for the parametric amplifier,” J. Opt. Soc. Am. B **10**, 2390–2393 (1993). [CrossRef]

**49. **E. H. Lieb and W. Liniger, “Exact analysis of an interacting Bose gas. I. The general solution and the ground state,” Phys. Rev. **130**, 1605 (1963). [CrossRef]

**50. **D. Muth, B. Schmidt, and M. Fleischhauer, “Fermionization dynamics of a strongly interacting one-dimensional Bose gas after an interaction quench,” New J. Phys. **12**, 083065 (2010). [CrossRef]

**51. **D. Muth, M. Fleischhauer, and B. Schmidt, “Discretized versus continuous models of p-wave interacting fermions in one dimension,” Phys. Rev. A **82**, 013602 (2010). [CrossRef]

**52. **D. Jaschke, M. L. Wall, and L. D. Carr, “Open source matrix product states: opening ways to simulate entangled many-body quantum systems in one dimension,” Comput. Phys. Commun. **225**, 59–91 (2018). [CrossRef]

**53. **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. Parcollet, G. Pawlowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwöck, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner, and S. Wessel, “The ALPS project release 2.0: open source software for strongly correlated systems,” J. Stat. Mech. **2011**, P05001 (2011). [CrossRef]

**54. **A. T. Sornborger and E. D. Stewart, “Higher-order methods for simulations on quantum computers,” Phys. Rev. A **60**, 1956 (1999). [CrossRef]

**55. **T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, “Dissipative Kerr solitons in optical microresonators,” Science **361**, eaan8083 (2018). [CrossRef]

**56. **R. El-Ganainy, K. G. Makris, D. N. Christodoulides, and Z. H. Musslimani, “Theory of coupled optical PT-symmetric structures,” Opt. Lett. **32**, 2632–2634 (2007). [CrossRef]

**57. **N. V. Alexeeva, I. V. Barashenkov, A. A. Sukhorukov, and Y. S. Kivshar, “Optical solitons in PT-symmetric nonlinear couplers with gain and loss,” Phys. Rev. A **85**, 063837 (2012). [CrossRef]

**58. **W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, “Optimal design for universal multiport interferometers,” Optica **3**, 1460–1465 (2016). [CrossRef]

**59. **M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, “Experimental realization of any discrete unitary operator,” Phys. Rev. Lett. **73**, 58 (1994). [CrossRef]

**60. **S. Olivares, “Quantum optics in the phase space,” Eur. Phys. J. Spec. Top. **203**, 3–24 (2012). [CrossRef]

**61. **H. A. Haus and Y. Lai, “Quantum theory of soliton squeezing: a linearized approach,” J. Opt. Soc. Am. B **7**, 386–392 (1990). [CrossRef]

**62. **S. J. Carter, P. D. Drummond, M. D. Reid, and R. M. Shelby, “Squeezing of quantum solitons,” Phys. Rev. Lett. **58**, 1841 (1987). [CrossRef]

**63. **Y. Lai and H. A. Haus, “Quantum theory of solitons in optical fibers. II. Exact solution,” Phys. Rev. A **40**, 854 (1989). [CrossRef]

**64. **L. Di Mauro Villari, D. Faccio, F. Biancalana, and C. Conti, “Quantum soliton evaporation,” Phys. Rev. A **98**, 043859 (2018). [CrossRef]

**65. **N. Korolkova, R. Loudon, G. Gardavsky, M. W. Hamilton, and G. Leuchs, “Time evolution of a quantum soliton in a Kerr medium,” J. Mod. Opt. **48**, 1339–1355 (2001). [CrossRef]

**66. **F. Singer, M. J. Potasek, J. M. Fang, and M. C. Teich, “Femtosecond solitons in nonlinear optical fibers: classical and quantum effects,” Phys. Rev. A **46**, 4192 (1992). [CrossRef]

**67. **A. Kenfack and K. Życzkowski, “Negativity of the Wigner function as an indicator of non-classicality,” J. Opt. B **6**, 396–404 (2004). [CrossRef]

**68. **L. G. Helt and N. Quesada, “Degenerate squeezing in waveguides: a unified theoretical approach,” J. Phys. Photon. **2**, 035001 (2020). [CrossRef]

**69. **M. Shirasaki and H. A. Haus, “Squeezing of pulses in a nonlinear interferometer,” J. Opt. Soc. Am. B **7**, 30–34 (1990). [CrossRef]

**70. **P. D. Drummond and H. He, “Optical mesons,” Phys. Rev. A **56**, R1107 (1997). [CrossRef]

**71. **M. G. Raymer, P. D. Drummond, and S. J. Carter, “Limits to wideband pulsed squeezing in a traveling-wave parametric amplifier with group-velocity dispersion,” Opt. Lett. **16**, 1189–1191 (1991). [CrossRef]

**72. **A. V. Buryak and Y. S. Kivshar, “Solitons due to second harmonic generation,” Phys. Lett. A **197**, 407–412 (1995). [CrossRef]

**73. **A. V. Buryak, P. Trapani, D. V. Skryabin, and S. Trillo, “Optical solitons due to quadratic nonlinearities: from basic physics to futuristic applications,” Phys. Rep. **370**, 63–235 (2002). [CrossRef]

**74. **G. Vidal and R. F. Werner, “Computable measure of entanglement,” Phys. Rev. Lett. **65**, 032314 (2002). [CrossRef]