Degenerate Squeezing in Waveguides: A Unified Theoretical Approach

We consider pulsed-pump spontaneous parametric downconversion (SPDC) as well as pulsed single- and dual-pump spontaneous four-wave mixing processes in waveguides within a unified Hamiltonian theoretical framework. Working with linear operator equations in $k$-space, our approach allows inclusion of linear losses, self- and cross-phase modulation, and dispersion to any order. We describe state evolution in terms of second-order moments, for which we develop explicit expressions. We use our approach to calculate the joint spectral amplitude of degenerate squeezing using SPDC analytically in the perturbative limit, benchmark our theory against well-known results in the limit of negligible group velocity dispersion, and study the suitability of recently proposed sources for quantum sampling experiments.


I. INTRODUCTION
First observed in the 1980s [1][2][3][4], squeezed states of light have become more ubiquitous over time [5]. While part of this is due to fundamental interest, there are many practical applications of squeezed light as well, including computation [6], communication [7], and metrology [8].
GBS involves a number of Gaussian states (i.e. those with a Gaussian Wigner function) input to an interferometer and subsequently detected at the output of the interferometer with photon number resolving (PNR) detectors. For non-classical input Gaussian states, sampling from the output distribution of photon number events at the detectors is conjectured to be a #P hard problem and thus intractable by classical computers for a large enough number of input states [19][20][21][22]. This is interesting from a theoretical point of view, for it stands to provide evidence against the extended Church-Turing thesis [23]. It is also experimentally interesting, for, in contrast to heralded single photon states, non-classical Gaussian states such as squeezed states can be created deterministically and so medium-scale implementations of GBS could be possible in the near future.
However, it should be noted that not every squeezed state of light is useful for GBS. In particular, the #P hardness depends on collecting samples in the photon number basis from states in a single known mode, and PNR detectors are not fast enough to provide spectral mode selectivity [24]. To get around this, a well defined spectral mode can be arranged by employing a pulsed pump in an appropriate nonlinear medium [25]. Common nonlinear media include second-and thirdorder nonlinear optical waveguides, capable of producing squeezed states via spontaneous parametric downconversion (SPDC) or spontaneous four-wave mixing (SFWM), respectively. We note that this situation is different from both the perturbative photon pair generation regime [26], as well as that of a continuous wave pump followed by homodyne detection [27]. Additionally, degenerate, rather than non-degenerate or so called twin-beam squeezing, is considered in studies of the hardness of GBS [9,10] and its applications [11][12][13][14][15][16][17][18].
In this work we therefore develop theoretical tools to calculate degenerate squeezing generated by pulsed pumps in waveguides. In contrast to some previous approaches, we focus on pulsed, rather than CW [28][29][30], pumps, avoid phase-space methods [31], and work with operator coupled mode equations driven by first-order temporal, rather than spatial, [32][33][34][35][36][37][38][39][40], derivatives. In particular, we note that coupled mode equations in which a series expansion in ∂ ∂t is used to accommodate material and modal dispersion, often used in the nonlinear optics community for classical fields, are not derived from a canonical formalism and thus it is not obvious that they are valid in the quantum regime [41] To provide a treatment consistent with Maxwell's equations and quantum mechanics we start from an established Hamiltonian formalism in Section II, where we use the Heisenberg equations of motion to obtain operator coupled mode equations. In subsections contained therein, we massage these equations into a form useful for numerical evaluation, introduce a simple method to include loss, and work out expressions for the phase sensitive and phase insensitive moments, generated state, and the state's associated joint spectral amplitude (JSA). Having established these tools, in Section III we consider three example calculations: SPDC in the low-gain regime, homodyne detection of single-pump SFWM as a function of gain, and dual-pump SFWM in the highgain regime. The first two examples serve as reasonable checks of our formalism, for they have been studied previously [26,42], while the final example is most relevant for GBS. We conclude in section IV.

II. FORMALISM
We follow a canonical formalism [43][44][45][46][47][48], which accounts for material and modal dispersion and yields Heisenberg equations of motion for the operator fields D (r) and B (r) that reproduce the dynamical Maxwell equations. This approach allows results of calculations to be absolute, rather than relative to an unknown constant, and facilitates simple comparisons with both experimental results and classical calculations. Waveguide propagation is assumed to be in the positive z direction. We begin with linear (L) and nonlinear (NL) Hamiltonians where and with Note that the Γ (n) (r) are nth-rank tensors linearly related to the more common susceptibility tensors χ (n) (r) (see Appendix A for details). We associate relevant waveguide modes with center wavevectors k J and frequenciesω J = ω Jk J . Taylor expanding where and introducing Fourier transformed operators centered at k J such that they are slowly-varying in space allows us to rewrite Eq. (2) as Note that these new operators satisfy simple commutation relations to excellent approximation (as they are exact only assuming that the modes exist for all k) for the ψ J (z) of interest. However the real power of the ψ J (z) is in the approximations they allow for D (r) and thus the nonlinear Hamiltonians. In general, the displacement field is expanded in terms of waveguide modes J in the transverse plane d Jk (x, y) as a Jk +H.c. (11) Using Eq. (8), however, we may form a useful approximate expression for Eq. (11) as (12) Because the nonlinearity is weak, while we include dispersion in the linear Hamiltonian, we comfortably neglect dispersion in the nonlinear Hamiltonians (aside from in the normalization of the d Jk J (x, y) [43,44]) and use only the first term of Eq. (12) in their construction. This allows simple identification of phase matched terms in nonlinear optical processes of interest.
With Hamiltonians defined, we are now in a position to study system dynamics. One approach is to solve the Schrödinger equation for the time evolution operator satisfying U(t, t) = I, and subsequently use this operator to transform kets in the usual way |Ψ(t) = U(t, t 0 ) |Ψ(t 0 ) . However, we can also transform Schrödinger operators to time dependent Heisenberg operators and solve their equations of motion instead. In the remainder of this Section we derive equations for these Heisenberg operators and show how to find time evolved kets using these solutions in the important case where the initial state |Ψ(t 0 ) = |vac .

A. Lossless coupled-mode equations
Here we use the expressions from the start of the Section to develop nonlinear Hamiltonians for degenerate squeezing processes, and subsequently use these Hamiltonians to derive operator coupled-mode mode equations in a common form. Details can be found in Appendix A.
For SPDC, we take J → {F, SH} for "fundamental" and "second harmonic", consider only the second-order nonlinearity H N L → H (2) N L , and, using only the first term of Eq. (12), write where ζ (2) (z) is a nonlinear coupling parameter defined in Appendix A, and we have assumed that the only phase (energy) matched process occurs for 2k F = k SH (2ω F =ω SH ). Note that the z-dependence of the nonlinear coupling parameter simply defines the nonlinear region of the waveguide. Coupled mode equations then follow immediately from the Heisenberg equations of mo- To avoid propagating oscillations at optical frequencies in the dynamics, we calculate the Heisenberg equations of motion for the slowly oscillating operators obtaining and where we have considered the driving field at the second harmonic to be a large classical field and thus replaced all instances ofψ SH (z) by its mean field ψ SH (z) and neglected any back-action on it fromψ F (z). For clarity of presentation, we do not include higher-order dispersion or the effects of self-or cross-phase modulation (SPM and XPM, respectively) here. However, we note that our formalism admits consideration of higher-order dispersion simply by including more terms in Eq. (9), as well as consideration of SPM and XPM by allowing N L . For single-pump SFWM, we have only a single mode J → P for "pump", consider only the third-order nonlinearity H N L → H J1J2J3J4 (z) is a nonlinear coupling parameter defined in Appendix A. From the Heisenberg equations of motion we find and where we have consideredψ P (z, t) as its mean field plus fluctuationsψ P (z, t) → ψ P (z, t) + δψ P (z, t), and neglected terms quadratic in the fluctuations or higher. For dual-pump SFWM, we take J → {P 1 , P 2 , S} for "pump 1", "pump 2" and "signal", consider only the third-order nonlinearity H N L → H J1J2J3J4 (z) is a nonlinear coupling parameter defined in Appendix A, and we have assumed 2k S = k P1 + k P2 and 2ω S =ω P1 +ω P2 . Coupled mode equations then follow from the Heisenberg equations of motion as and where we have considered the driving fields labeled by P 1 and P 2 to be a large classical fields and thus replaced all instances ofψ P1 (z, t) andψ P2 (z, t) by their mean fields ψ P1 (z, t) and ψ P2 (z, t) as well as neglected all terms quadratic inψ S (z, t) or higher.

B. Solutions
We note that equations for mean fields [Eqs. (18), (21), (24), and (25)] can be found using standard approaches, such as split-step Fourier or finite difference methods [49]. The remaining operator equations [Eqs. (19), (22), and (26)] are all of the form whereM (z, t) is real and positive. Defining where M (κ, t) = M * (−κ, t) and using [recall Eq. (8)] where we have put b J (κ, t) = a J(k J +κ) (t), we can rewrite Eq. (27) as a form highly amenable to numeric evaluation. Again we note that higher-order dispersion can be considered by including more terms in the expansion of the linear Hamiltonian of Eq. (9), thus resulting in more terms in ω (κ). To solve Eq. (30), we discretize κ j = j∆κ and write Introducing where as well as the double vector where the vector b(t) has entries b j (t) = b(κ j , t), then allows us to rewrite Eq. (30) compactly as Note that the matrix S = S T is symmetric and R = R † is Hermitian. For a small enough propagation forward in time ∆t this has solution where the (single-time) infinitesimal propagator is defined as Propagation over a finite interval can then be calculated by concatenating infinitesimal propagators K(t) to form the (two-argument) Heisenberg picture propagator [48] where t j = t 0 + j∆t and ∆t = (t − t 0 )/ , as and the product in Eq. (38) is taken in chronological order. Reverting back to continuous notation, we note that we may also write the solution of Eq.
We note that similar input-output relations have been derived previously [32,40,50,51], albeit with different labels.

C. Including loss
For unitary evolution (i.e. ignoring scattering loss), the evolution of the b (κ, t) themselves is sufficient to characterize system dynamics. However, once we include loss it becomes extremely useful to consider the phase insensitive and phase sensitive second moments with entries Using Eq. (36) one can easily write equations that update these moments as where V ≡ V (t) and W ≡ W (t) are the (singleargument) blocks of the propagator K(t). Note that in these equations there are inhomogeneous terms that drive vacuum fluctuations into the system and cause the moments to acquire nonzero values.
To include loss in a Hamiltonian formalism we need to add extra degrees of freedom to the system that will account for the modes into which the lost photons go. To this end we consider t as the time the waveguide modes interact with non-guided modes c(κ, t) that are initially in vacuum via a hopping interaction that removes photons from the b(κ, t) and places them in the c(κ, t). The operators of the system will transform as b(κ, t + ∆t) = √ ηb(κ, t) + 1 − η c(κ, t).
For the moments one finds where recall that at time t the scattering modes c are in vacuum. This recipe allows us to update the moments of the waveguide after experiencing a small amount of loss in the time interval ∆t. For larger time intervals, we simply concatenate Eq. (45) for many ∆ts. The key assumption for the validity of this evolution is that at each time interval ∆t the system of interest interacts with "fresh" non-guided modes that always come prepared in vacuum. This is essentially the same Markov approximation used to derive the Lindblad form Master Equation [52].
The only remaining question is how to set the value of η. Since we expect loss to compound exponentially it is natural to take where γ can be easily obtained from the standard attenuation constant α [49]. Numerically, combining the rules in Eq. (43) and Eq. (45) one can propagate the moments of the state (cf. Appendix B for a detailed derivation). Note that, because both squeezing and loss are Gaussian operations, if the state in the waveguide is Gaussian at the beginning of the propagation it will remain Gaussian throughout the evolution and thus be fully characterized by the N and M moments as well as the expectation value b(κ, t) . Finally, it is easy to show, after some algebra, that if the state at t 0 is |Ψ(t 0 ) = |vac then the state at some later time t is the squeezed state where the joint spectral amplitude is written in terms of the Schmidt functions and coefficients of the phase sensitive moment This is because the Hamiltonian generating U(t, t 0 ) [cf. Eq. (13)], even though nonlinear in the fields, is, upon replacing the pump(s) by its (their) classical expectation value(s), quadratic in the quantum bosonic operators [53,54]. The expressions developed here constitute our main results. Starting from an equation of the form of Eq. (27), describing pulsed pump SPDC or SFWM, they show how to include dispersion to any order, as well as SPM, XPM, and loss, in determination of the output state, its joint spectral amplitude, and phase sensitive and phase insensitive moments for arbitrary input pump power. These results provide complete information about the degenerate squeezed state of light generated in a waveguide. In what follows, we provide examples of their utility.

III. EXAMPLE CALCULATIONS
A. SPDC in the low-gain regime As a first application of our formalism, we consider SPDC in the low-gain regime, where results are well known and can be obtained analytically [26]. This serves as a basic check on the results of the previous section, and provides a simple rule of thumb for the grid size one should use to represent the moments when solving the problem numerically.
In Appendix B we derive the continuous evolution of the moments in this regime including loss. Solving the equations of motion between t 0 = −T /2 and t 1 = T /2 perturbatively we find (53) and assumed that the pump mean field undergoes only free propagation. Furthermore, we have introduced the phase-matching function and the complex energy difference and sum Note that the loss coefficients appear in two different places in the expression for b F (κ, t 1 )b F (κ , t 1 ) : they appear as a sum + in an overall exponential damping factor, but, more importantly, they also appear as a difference − in the sinh loss matching term [56,57]. When loss is negligible one can let T → ∞ and write where ∆ω = ω F (κ) + ω F (κ ) − ω P (κ ). With this simplification one can rewrite Eq. (50) as dropping an irrelevant phase factor. Finally, moving to frequency space yields the well known expression showing that the low gain JSA (recall that the JSA and the phase sensitive moment are equal in the low gain regime) is a product of the pump and phase-matching functions.

B. Single-pump SFWM quadrature variance
As a second check, we move away from analytic results and the perturbative regime, but neglect both loss and dispersion. We consider a homodyne measurement, in which one measures with for a normalized function dκ |φ (κ)| 2 = 1, and t f = L/v the time it takes light to traverse the length of the nonlinear region L. Writing this out in detail, we find where the final term sets the so-called shot noise level. We note that when φ (κ) corresponds to a Schmidt function ρ (l) (κ) [recall Eq. (48)], maximizing and minimizing (±) over θ, this reduces to the more common V ± = 2 sinh 2 (r l ) ± sinh (2r l ) + 1 = e ±2r l .
For numeric evaluation, we discretize φ j = φ (κ j ) √ ∆κ and rewrite Eq. (61) as We can check results obtained using Eq. (63) in the limit of no loss and no second-order dispersion (i.e. (64) Ignoring the z dependence of ζ (3) P P P P and putting √ ω P v P ψ P (z, t) → A P (z, t) such that |A P (z, t)| 2 has units of power, we find with solution where we have introduced the nonlinear parameter γ (see Appendix A) [49]. Following Ref. [42] (cf. Eq. (5.8a) there) the quadrature variance minimum can then be written as where Φ (t) = |A (0, t)| 2 γL is assumed peaked at t = 0 and C is a constant specific to pump pulse shape. In Fig. 1 we plot this analytic V − (solid lines) as well as results from our numerical routines (points) in dB as functions of Φ (0) for Lorentzian (C = 3/4), Sech (C = 4/5), Gaussian (C = 2/3), and Rectangular (C = 1) pulse shapes.

C. Dual-pump SFWM in the high-gain regime
Finally, we consider dual-pump SFWM in the highgain regime, the situation that is most relevant for GBS, as achieving degenerate squeezing via SPDC can require significant dispersion engineering, and single-pump SFWM demands interferometric removal of the pump from the generated squeezed state. As noted earlier [25], an ideal degenerate squeezing source for GBS applications should satisfy four desiderata, namely: (i) scalability, (ii) single-mode operation, (iii) squeezing levels sufficient to enable a genuine quantum advantage in computation, and (iv) compatibility with single photon and photon number-resolving detection. We investigate the suitability of a dual pump source with respect to (ii-iv). We note that integrated platforms are some of the best suited platforms for scalability, given their ease of fabrication and integration with other optical components.
We consider a silicon waveguide where one aims to generate squeezing around 1550 nm (corresponding to k ≈ 2π/1550 nm). For simplicity, we neglect linear and nonlinear loss, noting that adding such decoherence mechanisms will only degrade the quality of the source. We take the two pumps to each have a bandwidth ∆k = 41.8 cm −1 (corresponding to ∆λ ≈ 1.6 nm), and have centres separated by approximately 3∆k, with shapes similar to a top hat function (this corresponds roughly to the parameters used in Ref. [58]). Note that it is inappropriate to describe the pumps and signal field as three distinct modes in this situation, as the wavevector support of the pumps and the generated squeezed light significantly overlap and thus it does not makes sense to attach different mode labels to the pump and squeezed light. Based on this observation we directly use the the- Before entering the nonlinear region the field is prepared in a coherent state with amplitude where we have used the function f (k) = e −k 4 as a smooth non-Gaussian approximation to a top-hat function of unit width, k 2 − k 1 = 3∆k, k1+k2 1550 nm , N is the mean photon number in the pumps, and C is a normalization constant. We truncate the dispersion relation at the quadratic level as in Eq. (6) with the parameters v P = 7.019 × 10 7 m/s and v P = 4.711 m 2 /s.
In Fig. 2 we show the photon density of the mean field for different values of the gain. We parametrize this gain in terms of the mean photon number of the squeezed fluctuations generated in the field where n j = sinh 2 r j is the mean photon number in the j th Schmidt mode and we have written the destruction operator at wavevector k as its mean plus its fluctuations b(k) = b(k) + δb(k). The photon density in the mean can be contrasted with the photon density in the squeezed fluctuations defined in Eq. (70) and shown in Fig. 3. Note that there is significant overlap between the mean and the fluctuations; filtering is thus necessary to prevent the large amounts of energy in the mean from blinding the detectors necessary for applications in quantum sampling. Note, however, that the filtering of non-classical light comes with its own complications. In particular this will introduce mixedness in the state; this has been investigated in twin-beam sources in [59][60][61] and will be explored for degenerate squeezing in a future communication. The problem of overlap between the field and the squeezed fluctuations is further compounded by the fact that the joint spectral amplitude [cf. Eq. (48)] of the squeezed fluctuations is spectrally non-separable. This can be seen by looking at its absolute value for different gains as shown in Fig. 4. We numerically confirm that the JSAs are not separable by calculating the Schmidt number, defined as This quantity takes the value 1 for a separable JSA and increases for more entangled JSAs. For the source we model we find K ∼ 3.1 both in the low-and high-gain regimes (cf. Fig. 5). Comparing against recently proposed degenerate squeezing sources using microrings, for which K ∼ 1.0 [25], suggests that sources with parameters similar to those we have considered in this section are ill-suited for quantum sampling applications, and that the parameter space should be explored in a more thorough way to find better operating configurations.

IV. CONCLUSIONS
In conclusion, we have presented a unified Hamiltonian formalism for squeezed states of light generated via either SPDC or SFWM in waveguides using pulsed pumps. Our straightforward techniques directly admit dispersion to any order, linear losses, and self-and cross-phase modulation. They can be used to calculate standard expectation values, such as moments and spectra, but also states (including joint spectral amplitudes) and thus any desired observables.
We have provided examples of the application of our formalism to the well-known examples of SPDC in the low-gain regime and homodyne detection of single-pump SFWM, demonstrating its utility. All of the Python code used in this work is available on GitHub [55]. Furthermore, we have explored the suitability of a standard waveguide dual-pump SFWM source of squeezed light for GBS applications, finding that it is difficult to obtain a Schmidt number near unity in such a source. This points to a need to continue to explore more exotic and engineered sources going forward, a task for which our formalism is suited.

Single-pump SFWM
Our third-order nonlinear coupling parameter is given by where we have introduced jlmn (x, y)n 4 d j and put ignoring the contribution of χ (2) to Γ (3) [48]. Here s (z) defines the nonlinear region of the waveguide, andn and χ 3 are, respectively, a refractive index and third-order susceptibility introduced solely for convenience. We point out that ζ JJJJ (z) is connected to the more common nonlinear parameter [49] γ with units of m −1 W −1 via The general operator equation resulting from Eqs. (9) and (20) of the main text is However, replacingψ P (z, t) by its mean field plus fluctuationsψ P (z, t) → ψ P (z, t) + δψ P (z, t), and neglecting terms quadratic in δψ P (z, t) leads to Eqs. (21) and (22) of the main text.
Appendix B: Equations of motion for the moments Starting from Eq. (30) of the main text and using the product rule one can easily find the equations of motion for the quadratics b(κ, t)b(κ , t) and b † (κ , t)b(κ, t). In the Heisenberg picture they are where, to simplify notation, we have not written the second argument (time) of the operators b(κ, t) and b † (κ, t) and the functions S(κ, t) and M(κ, t). Note that these equations are equally valid for the the moments b(κ)b(κ ) and b † (κ )b(κ) which are precisely the expectation values of the quadratics above. Note the last (inhomogenoeus) term in Eq. (B2) appears due to normal ordering b(κ)b † (κ ) = δ(κ − κ ) + b † (κ )b(κ). This is the term that drives vacuum fluctuations into the system and that is responsible for generating squeezed vacuum if the input of the waveguide is vacuum. At the moment level one can add loss by simply adding to the equations of motion derived before the rate of change of the moments due to loss. This is easily found from Eqs. (45) and (46) of the main text to be We can put the unitary evolution due to the nonlinearity together with the rate of loss for the moment to write Assuming that at some earlier time t 0 = −T /2 one had b(κ, t 0 )b(κ , t 0 ) = b † (κ, t 0 )b(κ , t 0 ) = 0 one can solve the last equation perturbatively to a final time t 1 = T /2 b(κ, t 1 )b(κ , t 1 ) = i √ 2π t1 t0 dt e −i(ω(κ)+ω(κ )−iγ)(t1−t) × S(κ + κ , t).
Using this last expression in Eq. (B6) together with the definitions of the phase matching function in Eq. (54) and the complex energies in Eq. (55) one arrives at the expression of Eq. (50) for the perturbative joint spectral amplitude in SPDC after performing a straightforward time integration.