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$.

 

Fig. 1. (a) Finite space interval ${-}L/2 \le z \le L/2$ is discretized to $n$ spatial bins with size $\Delta z = L/N$. Pulse envelope function $f(z)$ (red lines) is discretized to form a vector ${\boldsymbol f}=(f_1,\ldots,f_N)^\top$ (orange boxes). Self-phase modulation (SPM) acts locally on discretized photon field (blue arrows), while dispersion mediates interactions between neighboring bins (brown arrows). (b) Example of how time evolution under the Hamiltonian Eq. (3) can be decomposed into one-mode SPM operations (blue boxes) and two-mode operations corresponding to dispersions (brown boxes). Spatial bins are labeled by corresponding annihilation operators for the purpose of illustration.

Download Full Size | PDF

A generic system state can be written as

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

To obtain ${\hat \rho _{\cal S}}$, we ideally would like to have a low-rank MPO representation of the operators

 

Fig. 2. Illustration of our supermode demultiplexing scheme. Application of a unitary transformation ${\hat R^{(r)}}$, composed of one-mode and two-mode gates $\{\hat R_r^{(r)}, \ldots ,\hat R_N^{(r)}\}$, demultiplexes $r$th supermode ${\hat A^{(r)}} = \sum\nolimits_m f_m^{(r)}{\hat a_m}$ to the $r$th (local) spatial bin. We show operations up to $r = 2$, where local annihilation operators of the first and second bins are transformed to ${\hat A^{(1)}}$ and ${\hat A^{(2)}}$, respectively. For the purpose of illustration, spatial bins are labeled by associated annihilation operators.

Download Full Size | PDF

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

Here, 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.

 

Fig. 3. (a)–(c) Quantum simulations of a pulse instantiated as a coherent-state sech soliton using MPS with bond dimension $\chi = 40$. Average photon number of $\bar n = 3.0$ is used with ($\kappa = 0.5$) and without ($\kappa = 0$) linear loss. MCWF with $M = 100$ quantum trajectories is used for the simulation with finite loss. (a) Time evolution of photon density distribution $\langle \hat \phi _z^\dagger {\hat \phi _z}\rangle$. (b) Wigner functions for the reduced density matrix ${\hat \rho _{\cal S}}$ for the sech supermode ${{\boldsymbol f}^{\text{sech}}}$. Upper and lower rows are for $\kappa = 0.0$ and $\kappa = 0.5$, respectively. (c) Time evolution of the purity $\text{Tr}(\hat \rho _{\cal S}^2)$ (left figure) and the doubled volume of the Wigner function (WF) negativity [67] (lower figure) of the soliton pulse state. (d) WF negativity of the soliton pulse state at various times as a function of $\bar n$. MPS-based simulation with $\chi = 50$ is used.

Download Full Size | PDF

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.

 

Fig. 4. Wigner functions of soliton pulse state with $\bar n = 6$ obtained using MPS-based full-quantum simulation with $\chi = 50$ (upper row) and the time-dependent Hartree approximation (TDH) following Eq. (26) (lower row). Orange ellipses represent $1/{e^2}$-width of the Gaussian Wigner functions of the solitonic mode calculated with linearized equations of motions [68,69].

Download Full Size | PDF

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]:

 

Fig. 5. Time evolution of the photon density distribution of a pulse instantiated as a coherent-state second-order soliton with mean photon number ${\bar n^{\text{(2sech)}}} = 8$. (a) Classically expected dynamics given by Eq. (27). (b) Full quantum simulation under the Hamiltonian Eq. (1) using MPS with bond dimension $\chi = 60$.

Download Full Size | PDF

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}}$).

 

Fig. 6. Example implementation of a short-time evolution under the ${\chi ^{(2)}}$ Hamiltonian Eq. (28) using TEBD. Discretized FH modes $\{{\hat a_1}, \ldots ,{\hat a_N}\}$ and SH modes $\{{\hat b_1}, \ldots ,{\hat b_N}\}$ are encoded on an MPS in an alternating manner. Swap operations ${\hat U_{\text{SWAP}}}$ are used to bring distant modes together. Two-mode operations ${\hat U_{\text{a},m}}$, ${\hat U_{\text{b},m}}$, and ${\hat U_{\text{NL},m}}$ are for FH dispersion, SH dispersion, and nonlinear three-wave mixing interaction, respectively.

Download Full Size | PDF

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

 

Fig. 7. Quantum simulations of pulses instantiated as a coherent-state simulton with FH photon number $\bar n = 4$. MPS with bond dimension $\chi = 50$ is used. (a) Time evolution of the FH and SH photon density distribution. (b) Negativity ${\cal N}({\rho _{\cal S}}^\prime)$ calculated for the two-mode reduced density matrix under linear transformations ${\hat \rho _{\cal S}}^\prime = {\hat W^\dagger}{\hat \rho _{\cal S}}\hat W$, where $\hat W$ is parametrized by a pair of angles $(\Phi ,\Theta)$. White cross indicates $({\Phi _0},{\Theta _0}) \approx (0.21\pi ,0.28\pi)$, which minimizes ${\cal N}({\hat \rho _{\cal S}}^\prime)$ for $t = 4$. (c) Wigner function of the single-mode reduced density matrix before ($\hat{\rho}_{\cal{A}\cal{B}}$) and after ($\hat{\rho}^\prime_{\cal{A}\cal{B}}$) applying $\hat W({\Phi _0},{\Theta _0})$ at $t = 4$. Labels ${\cal A}$ and ${\cal B}$ represent FH and SH, respectively.

Download Full Size | PDF

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.