1. Introduction

Spontaneous emission, where a single photon is emitted as an atom decays from an excited state to its ground state, was the earliest source of nonclassical light. Within a decade of the invention of the laser, however, researchers began to explore the possibility of using parametric fluorescence to generate other nonclassical states [1,2]. Here the state of a material medium is left unchanged after the passage of a strong pump pulse, but because of optical nonlinearities it is possible for a pair of new photons to be generated. Much of the interest has focused on two processes. In one of them, called spontaneous parametric downconversion (SPDC), a pump photon splits into two daughter photons, usually referred to as “signal” and “idler,” in what can be thought of as a “photon fission” process. This is sufficiently strong only in materials without inversion symmetry, where a $\chi _2$ nonlinear susceptibility is present. In the other process, called spontaneous four-wave mixing (SFWM), what can be thought of as an “elastic scattering” between two pump photons occurs to make two new photons, again usually referred to as “signal” and “idler,” with one of the new photons at a frequency higher than that of the pump photons, and the other at a frequency lower than that. This process is governed by a $\chi _3$ nonlinear susceptibility, and so in principle can occur in any material medium. In both cases the pairs of photons produced are typically entangled, and thus of interest in quantum information processing protocols.

It was soon realized that not only could a strong enough pump pulse generate a pair of photons with significant probability, but that multiple pairs might also be generated. The quantum superposition of these states with different numbers of pairs of photons leads to a state of light called squeezed, in that the uncertainty in the amplitude of one quadrature is suppressed below the usual quantum limit, whereas the uncertainty in the amplitude of the other is correspondingly increased so as not to violate the uncertainty relation [3,4]. These squeezed states are central to current efforts in quantum information processing based on continuous variables [5–12] and to quantum metrology [13–15].

One may think that the low efficiency of SPDC or SFWM is fundamentally related to their quantum nature. However, this is not the case, for it rather depends on the weakness of the nonlinear light–matter interaction. In this respect, there are at least two strategies to tackle this problem, either by increasing the nonlinear interaction length, e.g., by working with bigger and bigger crystals, or by enhancing the intensity of the electromagnetic field interacting with the matter.

In the case of SFWM, which is the weaker of the two above-mentioned nonlinear processes, significant progress was made at the beginning of this century with the use optical fibers [16,17]. These allow one to simultaneously enhance the electromagnetic field intensity by leveraging the transverse light confinement and, at the same time, increasing the nonlinear interaction lengths from a few millimeters to a few kilometers.

More recently, studies of squeezed light and its generation have intersected with impressive improvements in the design and fabrication of integrated photonic structures involving channel waveguides and resonators. In these systems, the spatial and temporal confinement of light can enhance the nonlinear interaction strength up to ten orders of magnitude with respect to what can be achieved in bulk systems, leading to efficient parametric fluorescence with continuous wave (CW) excitation at milliwatt pump powers.

One important consequence is that the generation of nonclassical light “beyond photon pairs” is becoming commonplace [18–26]. A second consequence is that a more detailed comparison of theory and experiment is now possible, because implementations of well-characterized integrated structures do not suffer from the uncertainties that can plague experiments with bulk systems, such as details of beam shape and propagation.

From a theoretical point of view, the description of the nonlinear light–matter interaction in micro- and nanostructures is far more complex than that required in bulk systems, starting with the quantization of the electromagnetic field and moving on to the construction of the Hamiltonian describing a system that can be composed of several optical elements. Finally, the strong enhancement of the nonlinear response of the system necessitates tools able to describe fast, nontrivial dynamics, which cannot always by approached using perturbation theories.

In this tutorial we present an overview of the physics of nonlinear quantum optics “beyond photon pairs.” We begin in Sec. 2 with a treatment of the quantization of light in integrated photonic structures. Although correct treatments of this subject can of course be found in the literature, there are pitfalls in adding to the linear Hamiltonian what one might think would be the obvious nonlinear contribution. We discuss the subtleties involved, and detail the mode expansion and nonlinear Hamiltonians for treating channel waveguides and microring resonators, the latter being the cavity structure we use as an example in this tutorial. In Sec. 3 we consider general “squeezing Hamiltonians” that result when the usual classical approximation is made for the pump pulse; these Hamiltonians then only involve sums of products of pairs of operators. We consider low-gain solutions, and then develop a general framework in which the solution of the Heisenberg equations of motion can be leveraged to build the full ket describing the squeezed state of light. This can be used to treat both spontaneous and stimulated parametric fluorescence, and we present results for both photon-number and homodyne statistics of the generated light. In Sec. 4 we focus on waveguides, and use a particular example to illustrate the low-gain regime and the joint spectral amplitude (JSA) of pairs of photons that can be generated, and how it is modified upon entering the high-gain region. In Sec. 5 we turn to microring resonators, and consider in detail both single- (SP-) and dual-pump (DP-)SFWM, as well as SPDC. With all this as background we turn in Sec. 6 to connections between classical and quantum nonlinear optics, and how the spontaneous processes we have considered can be understood as “stimulated” by vacuum power fluctuations. In Sec. 7 we move beyond the world of coherent states and squeezed states, both of which fall in the category of Gaussian states. We present a leading strategy, based on “matrix product states” (MPSs), to describe states more complicated than Gaussian. As moving beyond information processing that involves only Gaussian states and measurements is necessary to achieve universal quantum computing [27,28], the description and characterization of non-Gaussian states is a current area of active research, with many different approaches being explored [29–36]. Indeed, it can be hoped that at some point in the future a tutorial that could be considered a successor to this one, and might be called “Beyond Gaussian states,” will become available. Finally, we conclude in Sec.8.

2. Quantization in Integrated PhotonicStructures

In this section, we develop the quantization of nonlinear optics in an approach particularly suitable for treating integrated photonic structures. In linear quantum optics the usual treatments are naturally based on treating the electric and magnetic fields $\boldsymbol{E}(\boldsymbol{r},t)$ and $\boldsymbol{B}(\boldsymbol{r},t)$ as fundamental, but a straightforward extension of this to nonlinear optics is invalid [37]; it leads to incorrect expressions for the Hamiltonian that can be found even in well-respected references [38]. In fact, for nonlinear quantum optics a treatment going back to the early work of Born and Infeld [39], which treats $\boldsymbol{D}(\boldsymbol{r},t)$ and $\boldsymbol{B}(\boldsymbol{r},t)$ as the fundamental fields, seems the simplest way to develop the Hamiltonian treatment directly without first introducing an underpinning Lagrangian framework and building up the Hamiltonian treatment from that, as is done by, e.g., Drummond and Hillery [40]. In Sec. 2.1 we illustrate the problems that can arise if a Hamiltonian treatment of quantum nonlinear optics is based on $\boldsymbol{E}(\boldsymbol{r},t)$ and $\boldsymbol{B}(\boldsymbol{r},t)$ as the fundamental fields, and as a preamble to nonlinear quantum optics we develop linear quantum optics based on $\boldsymbol{D}(\boldsymbol{r},t)$ and $\boldsymbol{B}(\boldsymbol{r},t)$ as the fundamental fields in Sec. 2.2.

Of particular interest are channel waveguide and ring resonator structures, and after introducing the general quantization procedure we treat those special cases in Sec. 2.2.4. Initially we neglect dispersion, treating the relative dielectric constant $\varepsilon _{1}(\boldsymbol{r})$ as independent of frequency, but in Sec. 2.2.5 we indicate how dispersion can be included in frequency regions where absorption can be neglected, such as below the band-gap of a semiconductor, which is the usual regime of interest for integrated photonic structures. In Sec. 2.3 we present a treatment of the coupling between channel waveguides and ring resonators in the linear regime that will be useful in later sections, and in Sec. 2.4 we include nonlinearities by introducing the usual nonlinear susceptibilities. The particular results for the Hamiltonians describing self-phase modulation (SPM) and cross-phase modulation (XPM) in channels and rings are detailed in Secs. 2.4.1 and 2.4.2, respectively; in later sections of this tutorial we introduce the appropriate Hamiltonians for other nonlinear processes. Of course, other treatments are possible. One can start from the classical equations relating the nonlinear polarization and the electric field to each other, and “elevate” the classical fields to their quantum counterparts [41]. Or one could model a crystal as a lattice of isolated molecules and work out the nonlinear response of the crystal from that model [42]. However, we feel that the treatment given here provides both the trustworthiness and the generality appropriate for the problems at hand.

Throughout this tutorial we generally assume that the relative dielectric constant $\varepsilon _{1}(\boldsymbol{r})$, or, if dispersion effects are included, $\varepsilon _{1}(\boldsymbol{r};\omega )$, can be treated as a scalar. Structures of interest, such as those involving lithium niobate, are clearly exceptions to this. The generalization to a relative dielectric tensor is straightforward, but we do not do it explicitly here because it complicates the formulas and, we believe, makes the physics we want to highlight in this tutorial less accessible.

2.1 Trouble with $E$

The quantization of the electromagnetic field in vacuum is a standard exercise in elementary quantum optics, and is naturally formulated in terms of the electric and magnetic fields $\boldsymbol{E}(\boldsymbol{r},t)$ and $\boldsymbol{B}(\boldsymbol{r},t)$ [43]. The generalization to include a uniform, frequency-independent dielectric constant is just as easy, and including even a position dependent relative dielectric constant $\varepsilon _{1}(\boldsymbol{r})$, where the constitutive relations are taken to be

For our discussion in this section we just quote the results for a dielectric ring shown in Fig. 1; we give a derivation of the results in Sec. 2.2.4. Cylindrical coordinates $z$, $\rho =\sqrt {x^{2}+y^{2}}$, and the angle $\phi$ are the natural variables; it is convenient to introduce a nominal radius $R$ and use $\zeta =R\phi$ in place of $\phi$. Using the notation $\boldsymbol{r}_{\perp }=(\rho,z)$, a volume element is $\textrm{d}\boldsymbol{r}=\textrm{d}\boldsymbol{r}_{\perp }\textrm{d}\zeta$, where $\textrm{d}\boldsymbol{r}_{\perp }=R^{-1}\rho \textrm{d}\rho \textrm{d}z$, and $\zeta$ varies from $0$ to $\mathcal {L}=2\pi R$. With the relative dielectric constant $\varepsilon _{1}(\boldsymbol{r})=\varepsilon _{1}(\boldsymbol{r}_{\perp })$ in the ring large enough to confine modes, the modes will be labeled by integer quantum numbers $m_{J}$, where the wavenumber associated with propagation around the ring is $\kappa _{J}=2\pi m_{J}/\mathcal {L}$; for the moment we restrict ourselves to one relevant transverse field structure for each $\kappa _{J}$, identify its frequency by $\omega _{J}$, and consider only these bound modes. Then introducing raising and lowering operators $c_{J}^{\dagger }$ and $c_{J}$ for each mode, neglecting zero-point energies the Hamiltonian is

 

Figure 1. Schematic representation of an isolated ring resonator. The inset shows the typical field profile of the fundamental transverse-electric mode supported by the bent waveguide. In the simulation the waveguide is 800 nm wide and 400 nm high, is made of silicon nitride fully clad in silica, and has a bending radius of 10 $\mu$m. The simulation is performed at a wavelength of 1550 nm. The contour plot accounts for the $E_{\zeta }$ component, whereas the arrow field represents $\boldsymbol{E}_{\boldsymbol{r}_\perp }$

Download Full Size | PDF

This is all uncontroversial. All of the possible derivations of (2) and (3) rely on identifying the energy density $\mathfrak {h}(\boldsymbol{r})$ at a point $\boldsymbol{r}$ in space, and noting that for a change in the fields at this point the energy density changes according to [44]

However, what would seem to be a straightforward extension of this to include nonlinear effects, where then

However, this leads to nonsense. To see that, consider just one mode $P$ and its nonlinear interaction with itself. Neglecting normal ordering corrections in the nonlinear term in (9), we find that the equation for $c_{P}$ is

What has gone wrong? The Hamiltonian density (9) is certainly equal to the energy density, as it should be for canonical quantization. However, the use of the form (3) of the electric and magnetic fields in that Hamiltonian, together with the usual commutation relations for the operators $c_{P}$ and $c_{P}^{\dagger }$, does not lead to the correct equations, Maxwell’s equations, in the nonlinear regime. This is apparent because, very generally, the equation for the Heisenberg operator $\boldsymbol{B}(\boldsymbol{r},t)$ would be

A hint of the problem can be gleaned from the differential form (4) for the change in the energy density, which suggests that $\boldsymbol{B}$ and $\boldsymbol{D},$ rather than $\boldsymbol{B}$ and $\boldsymbol{E}$, should be thought of as the fundamental fields. This is in contrast to a Lagrangian approach, where the differential form for the change in the Lagrangian density is given [39] by

2.2 Linear Quantum Optics

We begin by adopting the constitutive relations (1), assuming a real relative dielectric constant $\varepsilon _{1}(\boldsymbol{r})$ that is dispersionless but arbitrarily position dependent. In Sec. 2.2.5 we generalize to include material dispersion, which can have important consequences for experiments in quantum photonics.

2.2.1 Modes

We begin with the classical Maxwell equations

We begin by looking at modes that are also chosen to satisfy periodic boundary conditions over a normalization volume $V$. Then the operator $\mathcal {M}$,

2.2.2 Canonical Formulation

It is useful to use a slightly more general canonical formulation than often employed, but we begin by reviewing the traditional approach. There one starts with a set of functions of the canonical coordinates and momenta, $q_{i}$ and $p_{j}$ respectively, and with the Poisson bracket of two such functions $f(q_{i},p_{j})$, $g(q_{i},p_{j})$ defined according to

As an example, suppose that instead of deriving (27) from (25) we simply assert (27) for the $q_{i}$ and the $p_{j}$, and take our set of quantities to be the $q_{i}$, the $p_{j}$, and all polynomial functions of them. Then, together with the definition that the Poisson bracket of two numbers vanishes, and the definition that the Poisson bracket of a number with any of the quantities vanishes, the Poisson brackets of all such quantities can be determined from (26) and (27), and the set of quantities is closed under those operations. If the Hamiltonian is taken as a polynomial function of the $q_{i}$ and the $p_{j}$, then the dynamics given by (28) are precisely what we would expect from the more traditional approach. This more general approach is often useful. One example is if one wants to treat the angular momentum components $J_{i}$ of an object as fundamental, as one does for the spin of a particle; then there are no canonical position and momentum in the description, but the Lie algebra involves the Poisson bracket

That approach is also useful here. Returning to our Eqs. (16), we note that the divergence equations can be taken as initial conditions, for if they are satisfied at an initial time and the curl equations are satisfied at all later times, the divergence equations will also be satisfied at all later times. Thus, we need only construct dynamical equations for the curl equations, and for ultimate use in quantization we seek to do this with a Hamiltonian that is numerically equal to the energy. From (5) we see that the Hamiltonian should then be

As first pointed out by Born and Infeld [39], to obtain the correct dynamical equations the fundamental Poisson brackets should be taken as

2.2.3 Quantization

We now quantize in the usual way, by taking

However, because those field operators must be Hermitian, the operators $\mathcal {C}_{\alpha }^{(1)}(t)$ and $\mathcal {C}_{\alpha }^{(2)}(t)$ are not all independent; we easily find the conditions

2.2.4 Special Cases

The usual first case of interest is a uniform medium. We treat that for completeness in Appendix A. Here we consider the form of the field operators for two simple structures from integrated optics, the first a channel waveguide as indicated schematically in Fig. 2, and the second a ring resonator, initially considered isolated and indicated schematically in Fig. 1; of interest in itself, the ring resonator can also be taken as a simple example of an integrated optical cavity.

 

Figure 2. Schematic representation of a channel waveguide. (a) Typical field intensity profile of the fundamental transverse-electric mode, computed by an eigenmode solver. The waveguide and simulation parameters are the same as in Fig. 1, but with no curvature. By comparison with Fig. 1 it is possible to appreciate the different field distribution produced by the waveguide curvature. The contour plot accounts for the $E_{x}$ component, whereas the arrow field represents $\boldsymbol{E}_{y,z}$. (b) Scanning electron microscopy (SEM) image of a typical high index contrast waveguide cross section. Reprinted from [54]. Copyright 2020 Optical Society of America.

Download Full Size | PDF

2.2.4.1 Channel Waveguide

Looking first at the waveguide structure, we take $\boldsymbol{\hat {z}}$ normal to the substrate and $\boldsymbol{\hat {x}}$ the direction along which the waveguide runs; then $\varepsilon _{1}(\boldsymbol{r})=\varepsilon _{1}(y,z)$ here.

We are interested in modes confined to the waveguide, and so we need only introduce periodic boundary conditions in the $\boldsymbol{\hat {x}}$ direction, with periodicity over a length $L$. Then because of the translational symmetry in the $x$ direction we can choose the modes $(\boldsymbol{D}_{\alpha }(\boldsymbol{r}),\boldsymbol{B}_{\alpha }(\boldsymbol{r}))$ satisfying (18) to be of the form

2.2.4.2 Ring Resonators

Next we turn to the ring structure shown in Fig. 1, and consider the form of the modes $\left (\boldsymbol{D}_{\alpha }\left (\boldsymbol{r}\right ),\boldsymbol{B}_{\alpha }\left (\boldsymbol{r}\right )\right )$ that can at least be approximated as confined to the ring. Here we use indices $J$ to denote the modes, and $\omega _{J}$ their frequencies; the modes can be taken to be of the form

The electromagnetic field operators associated with the modes (56) can then be written, using (38), as

The general normalization conditions (23) and (24) here become

For the curvature of the ring not too great, the fields $(\textsf{d}_{J}(\boldsymbol{r}_{\perp };\zeta ),\textsf{b}_{J}(\boldsymbol{r}_{\perp };\zeta ))$ will be related to those fields $(\boldsymbol{d}_{Ik}(y,z),\boldsymbol{b}_{Ik}(y,z))$ of a channel structure of the same cross section according to

2.2.5 Including Dispersion

Up to this point we have assumed that the relative dielectric constant $\varepsilon _{1}(\boldsymbol{r})$ is independent of the frequency of the optical field. More generally, however, if we have electric and displacement fields

We begin by noting that for $\varepsilon _{1}(\boldsymbol{r};\omega )$ purely real we can still consider “modes” for $\boldsymbol{D}(\boldsymbol{r},t)$ and $\boldsymbol{B}(\boldsymbol{r},t)$, as we did when dispersion was neglected. Seeking solutions of Maxwell’s Eqs. (16) of the form (17), but now using the constitutive relations

In Appendix A we consider the inclusion of dispersion in a description of light in a uniform medium. Here we begin with the inclusion of dispersion in the description of channel waveguide structures, first considering periodic boundary conditions. The mode frequencies $\omega _{Ik}$ must be found by solving (66) with modes of the form (43); following the discussion above, it is easy to show from (67) that the normalization condition (45) is then replaced by

Finally, for the isolated ring the inclusion of dispersion leads to using the mode forms (56) in (66) to find the mode frequencies $\omega _{J}$ and the fields $\textsf{d}_{J}(\boldsymbol{r}_{\perp };\zeta )$ and $\textsf{b}_{J}(\boldsymbol{r}_{\perp };\zeta )$. Instead of the normalization condition (61) we use

2.3 Channel Fields and Ring–Channel Coupling

We now turn to the coupling between the channel and ring, in a structure such as in Fig. 3. We identify the isolated ring resonances $\omega _{J}$, restricting ourselves here to resonances associated with one transverse mode; recall that the isolated ring fields are given by (57), with in general the normalization condition (75). Our electromagnetic fields associated with the channel are expanded as (54), with commutation relations (53) and a Hamiltonian (55); the mode normalization condition is generally given by (70). In those expressions we have used the dummy index $I$ to indicate the mode type. For the channel we now just consider one type, and use the dummy index $J$ to label a range of $k$ (with a corresponding range of $\omega _{Jk}$), introducing different wavenumber ranges, each centered at a ring mode resonance.

 

Figure 3. Schematic representation of a point-coupled ring resonator. The coupling is modeled by a simple $\begin {pmatrix} \sigma & i\kappa \\ i\kappa & \sigma \\ \end {pmatrix}$ matrix, where $(\sigma, \kappa ) \in \mathbb {R}$ are the self-coupling and cross-coupling coefficients, respectively. The coupler is assumed lossless, hence yielding $\sigma ^2 + \kappa ^2 = 1$.

Download Full Size | PDF

If we denote the range of $k$ allotted to $J$ by $\mathscr {R}(J)$, (see Fig. 4) we can then write

 

Figure 4. Typical transmission of a point-coupled ring resonator. Here $\mathscr {R}(J)$ represents the wave vector range associated with each resonance $J$, centered at $k_J$.

Download Full Size | PDF

Under certain conditions the Hamiltonian (79) and the derivation of the dynamics of $\psi _{J}(x,t)$ that follow from it can be simplified. We clearly have

To lowest order in $\Delta k/\Delta \kappa$ we can write the electromagnetic field operators (76) in terms of the channel field operators by neglecting the variation of $\omega _{Jk}$, and $\boldsymbol{b}_{Jk}(y,z)$ and $\boldsymbol{d}_{Jk}(y,z)$, over the $\Delta k$ characterizing the excitation bandwidth. Putting $\boldsymbol{b}_{J}(y,z)\equiv \boldsymbol{b}_{Jk_{J}}(y,z)$ and $\boldsymbol{d}_{J}(y,z)\equiv \boldsymbol{d}_{Jk_{J}}(y,z)$, we have

We now consider the channel and ring close enough that coupling can occur. The coupling arises because the evanescent fields associated with the ring and the channel overlap, and light can move from the channel to the ring and vice versa. The standard phenomenological model for this coupling [59] begins by treating the field in the ring in the same manner as one treats the field in the channel; in our notation this involves introducing a field $\psi _r(\zeta,t)$, where $\zeta =0$ at the nominal coupling point, to describe the propagation in the ring. If we let $x=0$ denote the nominal coupling point for the channel field $\psi _J(x,t)$ associated with the ring resonance of interest (see Fig. 4), we then take the Fourier components of $\psi _J(x,t)$ and $\psi _r(\zeta,t)$,

 

Figure 5. Enlarged view of the point coupler of Fig. 3 between a channel and a ring, with particular emphasis on the conventional coordinate systems used in the two waveguide.

Download Full Size | PDF

Here we assume that the waveguide properties of the ring and channel are the same, although the more general case can easily be considered. Conventionally both $\sigma$ and $\kappa$ are taken to be real and positive; then energy conservation yields

We are typically interested in frequencies $\omega$ such that $|\omega -\omega _{J}|$ is much less than the frequency difference between neighboring resonances; then we have $|\eta | \ll 1$, and we can approximate $e^{i\eta } \approx 1+i\eta$. We are also typically interested in small ring–channel coupling, such that $\kappa \ll \sigma$; then because $|\eta |\ll 1$, we can take $\eta \sigma \approx \eta$, and using these approximations in (103) we have

Although a whole range of linear and nonlinear optical phenomena in channel–ring structures can be considered at this level [59], for quantum nonlinear optics we want to construct a quantum treatment of the channel–ring coupling that can be built into a Hamiltonian framework. Beginning with the linear regime, we can do this by formally introducing a point coupling model [60], with a coupling constant $\gamma _{J}$. The full linear Hamiltonian is taken to be

The Heisenberg equations of motion following from the Hamiltonian (105) then are

If we then introduce Fourier components of $\psi _{J}(0^{\pm },t)$ as in the first of (97), we find [60]

 

Figure 6. Schematic representation of a point-coupled ring resonator with and additional fictitious coupler and waveguide, referred to as “phantom channel.”

Download Full Size | PDF

The solution of our general equations (118), again following the same procedure as used above, gives

To avoid dealing with the fields $\psi _J(x,t)$ that have a discontinuity at $x=0$, it is convenient to introduce fields $\psi _{J<}(x,t)$ and $\psi _{J>}(x,t)$ associated with each $\psi _J(x,t)$. The first of these is equal to $\psi _J(x,t)$ for $x<0$ and, for $x\ge 0$, follows from the evolution equation (84) that would hold were there no coupling to the ring. Similarly, $\psi _{J>}(x,t)$ is given by $\psi _J(x,t)$ for $x>0$ and, for $x\le 0$, follows from that evolution equation. Thus, the fields $\psi _{J<}(x,t)$ and $\psi _{J>}(x,t)$ are continuous even at $x=0$, and we have $\psi _J(0^-,t)=\psi _{J<}(0,t)$ and $\psi _J(0^+,t)=\psi _{J>}(0,t)$.

Finally, because the free spectral range (FSR) in the neighborhood of the $J$th resonance is given by

2.4 Including Nonlinearities

We now turn to including the effects of nonlinearities. Considering a very general dielectric structure, and beginning in the classical regime, the usual description of the nonlinear response is in terms of nonlinear response coefficients $\chi _{2}^{ijk}(\boldsymbol{r})$ and $\chi _{3}^{ijkl}(\boldsymbol{r}),$ etc., with the components of the nonlinear polarization $\boldsymbol{P}_{\text {NL}}(\boldsymbol{r},t)$ given by

Equation (135) is a usual starting point for investigating the propagation of light in nonlinear media. However, because our fundamental fields are taken to be $\boldsymbol{D}(\boldsymbol{r},t)$ and $\boldsymbol{B}(\boldsymbol{r},t)$, we want to rewrite this as an expression for $\boldsymbol{E}(\boldsymbol{r},t)$ in terms of $\boldsymbol{D}(\boldsymbol{r},t).$ To do that, write (135) as

Using an expansion of $\boldsymbol{D}(\boldsymbol{r})$ in terms of the linear modes of the system the form of ${\textrm{H}^{\textrm{NL}}}$ in terms of raising and lowering operators can be explicitly constructed. We do this later in this tutorial for SPDC and SFWM. In the rest of this section we treat the special cases of SPM and XPM. For simplicity we assume that $\Gamma _{3}^{ijkl}$ receives significant contributions only from $\chi _{3}^{ijkl}$, and not from $\chi _{2}^{ijk}$ [see (138)]; that usually holds, but if it does not the generalization of the formulas given in the following is straightforward.

2.4.1 SPM and XPM in Channels

Considering an isolated channel, we still assume that we can take our fields to be of the form (92), although the frequencies $\omega _{J}$ here are not linked to any ring resonances, but merely identify the center frequencies of different components of light that are well-separated from each other in frequency. Assuming that we have just one such component; taking $J=P$, its self-interaction is described by $\chi _{3}^{ijkl}(\boldsymbol{r})$, and using (92) in (142) the term describing the resulting SPM is

Taking as our Hamiltonian the sum of the linear Hamiltonian (90) and the nonlinear contribution (144) we find the Heisenberg equations of motion yield

To write ${\gamma_{\textrm{chan}}^{\textrm{SPM}}}$ in a more transparent form, divide the expression (145) for it by the square of

We now consider XPM, a consequence of the nonlinear interaction between a field $\psi _{S}(x,t)$ with frequency centered at $\omega _{S}$ and a field $\psi _{P}(x,t)$ with frequency centered at $\omega _{P}$. The contribution to the nonlinear Hamiltonian responsible for XPM is

2.4.2 SPM and XPM in Rings

In a ring coupled to a channel, the usual approximation is to assume that the important nonlinear interaction is in the ring, because in usual applications the field is concentrated there. Thus, to construct the appropriate Hamiltonian for SPM here we use the expression (57) for $\boldsymbol{D}(\boldsymbol{r})$ in the ring, keeping the contribution from just one mode $P$ with frequency $\omega _{P}$, in the general nonlinear Hamiltonian (142). Then keeping the contribution from $\chi _{3}$, the form of the nonlinear Hamiltonian closely follows that for a channel, and we have

The treatment of XPM also follows that of the channel. Considering ring fields at $\omega _{P}$ and $\omega _{S}$, the contribution to the nonlinear Hamiltonian associated with XPM is

We close this section by deriving the promised correct expression for the SPM term arising in Sec. 2.1. In the special case considered there, we can identify $\varepsilon _{3}(\boldsymbol{r})$ with $\chi _{3}^{zzzz}(\boldsymbol{r}),$ and from (163) we find

3. Obtaining the Ket from the Equations of Motion

In this section we develop mathematical tools to describe the quantum states that are typically generated in quantum nonlinear optical processes. The Hamiltonians describing these processes are polynomials with degree strictly greater than two in the electromagnetic fields, such as those developed in the previous section. However, it is often the case that some of the fields can be described classically, that is, replaced by their classical expectation values, leading to Hamiltonians that, while still nonlinear in the fields, are only quadratic in quantum operators. The validity and usefulness of this approximation is made apparent in subsequent sections, and thus for the moment we focus on the mathematical tools necessary to obtain the quantum state generated by these interactions. We stress that although there is a considerable reduction in complexity working with Hamiltonians at most quadratic in the operators, there are still significant complications stemming from the fact that these Hamiltonians are still time-dependent and have nonvanishing commutators at different times.

Before heading on directly to the most general problem one can tackle with the tools from Gaussian quantum optics and the symplectic formalism [5,64–67] we consider simple examples that help build intuition toward the most general results in Sec. 3.1. We assume the reader is familiar with basic quantum mechanics, including the Schrödinger and Heisenberg pictures [68,69]. In Sec. 3.2 we introduce the Dyson, Magnus, and Trotter–Suzuki expansions as methods for the solution of linear differential equations, including the Schrödinger and the Heisenberg equation. In Sec. 3.3 we use the first-order Magnus expansion to solve the low-gain multimode squeezing problem and to introduce the concept of Schmidt modes. In Sec. 3.4 we use the fact that the solutions to the Heisenberg equations must preserve equal-time commutation relations to derive important properties of the Heisenberg propagator. In Sec. 3.5 we introduce the characteristic function and the Heisenberg propagator to derive the form of the quantum state generated by a general quadratic Hamiltonian in the spontaneous regime when the initial state is vacuum. We specialize some of the results derived to the case of nondegenerate squeezing in Sec. 3.6. Then we generalize the results related to the spontaneous problem, when the input state is vacuum, to arbitrary input states in Sec. 3.7. In the last two subsections we study the effect of loss in Sec. 3.8 and finally show how to calculate the statistics of homodyne and photon-number resolved measurements in Sec. 3.9.

3.1 Squeezing Hamiltonians

To gain some insight into the Gaussian formalism we study a very simple and well-known problem in quantum optics. We solve, both in the Schrödinger and Heisenberg pictures, the dynamics generated by the Hamiltonian

For the very simple time-independent Hamiltonian of (174), we can immediately write

 

Figure 7. Wigner functions of the vacuum state and a squeezed vacuum state.

Download Full Size | PDF

One can also solve the dynamics of this problem in the Heisenberg picture. The operators are now dynamical quantities defined by $O(t) = \mathcal {U}^\dagger (t,t_{\text {i}})O(t_{\text {i}}) \mathcal {U}(t,t_{\text {i}})$ and their dynamics is determined by the Heisenberg equation of motion,

We can immediately find the solution of the last differential equation by exponentiation,

At this point it might seem like overkill to have derived the same expectation value using two different methods. However, as we show in a moment, it is useful to have two ways of looking at the problem because there are situations where solving the Schrödinger equation for the ket $| {\Psi (t)}$ is hopeless whereas solving the Heisenberg equations is practical. Moreover, once the Heisenberg equations are solved one can then write the sought after time-evolved Schrödinger picture ket.

3.2 Multimode Squeezing: Dyson, Magnus, and Trotter–Suzuki

Let us consider a slightly more complicated Hamiltonian given by

In this section we do not attach any particular meaning to the mode labels $i$ and $j$ in the last set of equations. In subsequent sections these labels will correspond to, for example, discretized wave vectors near a reference or central wave vector used to obtain a simplified dispersion relation. In this case, the term $\Delta _{ij}$ will correspond to a detuning between the different wave vectors (and possibly XPM of a classical pump on the quantum modes) whereas $\zeta _{kl}$ will correspond to wave-mixing induced by a nonlinear process.

As before, we would like to find the ket obtained by evolving the vacuum under this Hamiltonian, the so-called “spontaneous problem.” Note that the Hamiltonian describing the dynamics is now time-dependent and has three different terms that do not commute with each other, making a more complex dynamics than that of the single-mode problem analyzed in Sec. 3.1. In what follows we use boldface to refer to matrices $\boldsymbol{M}$ while referring to their entries as $M_{ij}$. In matrix notation, the constraints described in the previous paragraph are simply $\boldsymbol{\Delta } = \boldsymbol{\Delta }^\dagger$ and $\boldsymbol{\zeta } = \boldsymbol{\zeta }^T$ where $\dagger$ and $T$ indicate the conjugate transpose and the transpose, respectively.

We can formally solve the dynamics by writing the time-evolution operator associated with this Hamiltonian as

In practice, the expression (190) can be interpreted as a power series of the Dyson [72] or Magnus [73,74] type or can be approximated using Trotterization [75,76]. In the Dyson series the unitary evolution operator is written as an infinite power series

For the quadratic Hamiltonians we investigate here, the unitary evolution operator resulting from the Magnus expansion at any order respects the photon statistics that the exact operator is known to have, whereas that resulting from the Dyson expansion does not. Furthermore, in the Magnus expansion all terms beyond $\Omega _1$ are clearly time-ordering corrections, related to the fact that $\left [ H(t), H(t') \right ] \neq 0$ whereas in the Dyson series each successive term contains a part associated with time-ordering corrections as well as a part associated with the solution that would be obtained by ignoring time-ordering corrections and naïvely writing $\mathcal {U}(t,t_{\text {i}}) = \exp \left (-\frac {i}{\hbar } \int _{t_{\text {i}}}^t \textrm{d}t' H(t') \right )$.

Finally, a third approach is to use the Trotter–Suzuki expansion to approximate the time evolution operator as a product of time evolution operators $\mathcal {U}(t_i,t_i+\Delta t)$ where each Hamiltonian is assumed to be roughly constant over small intervals of duration $\Delta t$, thus writing

3.3 Low-Gain Solutions

The three strategies in the previous section have different strengths and weaknesses. If the Hamiltonian satisfies $[H(t),H(t')]=0$ for all the times $t,t'$ between $t_{\text {i}}$ and $t_{\text {f}}$ then the Magnus expansion at the first level immediately gives the exact solution. If the Hamiltonian $H(t)$ is “weak” we can treat the Magnus or Dyson series perturbatively and keep only the first few terms of either of them.

The simplest approximation is to keep only the first term of either of them. In this perturbative or low-gain limit one can obtain a simple solution to the spontaneous problem

Using the Takagi–Autonne decomposition lets us introduce the Schmidt (or broadband, or supermode) operators

We can also use the form of the low-gain time-evolution operator to transform the operators. Using the definitions of the Schmidt modes one can easily find

3.4 Properties of the Solution to the Heisenberg Equations of Motion: Bogoliubov Transformations

Let us now analyze this same problem in the Heisenberg picture for arbitrary gain. We now find

From the fact that the transformations connecting initial- and final-time Heisenberg operators are linear and the fact that these solutions must respect the canonical commutation relation at each time, one can obtain a number of useful properties. To simplify notation, from now on we drop initial and final times from matrices, writing, for example, $\boldsymbol{K}$ in place of $\boldsymbol{K}(t_{\text {f}},t_{\text {i}})$. Let us start by writing the Heisenberg propagator in block form

It is useful to note that the entries of the matrix

3.5 Obtaining the Ket from the Bogoliubov Transformations

In the previous section we showed how to solve the Heisenberg equations of motion, and, based on the fact that they correspond to unitary operations in Hilbert space, we derived some of their properties. In this section, we show how one can obtain the state that evolved under $\mathcal {U}$. Before doing this we introduce a useful object that uniquely describes a quantum state, its characteristic function [70]

For the vacuum state the characteristic function is easily calculated, and found to be a Gaussian in the complex amplitudes $\boldsymbol{\alpha }$

 

Figure 8. Block diagram describing how to calculate the state generated by a quadratic bosonic Hamiltonian.

Download Full Size | PDF

3.6 Nondegenerate Squeezing

As we show later, it is often useful to split our $\ell$-modes into two separate groups, called the signal and idler modes, for which we write ladder operators $b_i$ and $c_j$, respectively, satisfying canonical commutation relations

The Hamiltonian in (238) looks superficially different to the one in (188) but it is easy to show that the former is just a special case of the latter by identifying

We can now move to the high-gain regime, where we find that the solution to the spontaneous problem has exactly the same form as in the low-gain regime, but now the joint amplitude of the nondegenerate squeezed state will be inferred from the SVD of a second-order moment of the operators. As previously, we write the equations of motion for $\boldsymbol{b}$ and $\boldsymbol{c}^\dagger$ as

 

Figure 9. Experimental and theoretical prediction for the blocks of the transfer function $\boldsymbol{K}^{b,c}$ for nondegenerate squeezing in a periodically poled potassium titanyl phosphate SPDC source. The labels $U_{s,s}, U_{i,s}, U_{s,i}, U_{i,i}$ in the figure correspond to $V_{b,b}, W_{c,b}, W_{b,c}, W_{c,c}$ in the notation used in this paper. Figure 9 reprinted with permission from Triginer et al., Phys. Rev. X 10, 031063 (2019) [92]. Copyright 2019 by the American Physical Society.

Download Full Size | PDF

Using this input–output relation we can now calculate the (nonzero) second-order moments

Note that the output of the spontaneous problem is also a Gaussian state, characterized by its second-order moments and that we can write the output ket as

3.7 Solution to the Stimulated Problem

Having solved the spontaneous problem associated with a general quadratic Hamiltonian, we now solve the associated stimulated problem, i.e., we study the evolution of other states different from vacuum. One can write an arbitrary input state as

3.8 Passive Operations and Loss

An important set of processes in the evolution of bosonic modes are given by passive operations, i.e., linear operations where the energy is not increased. In these process, the modes of an environment prepared in the vacuum couple to the modes of interest via a beam-splitter-like interaction, which in the Heisenberg picture we write as [95–97]

3.9 Photon-Number and Homodyne Statistics

In this section we study properties of the photon number and quadrature statistics of the Gaussian states described in the previous sections. We first look at pure lossless states of the form given in (235) and ask what can be said about the statistical moments of the total photon number observable

Let us apply these results to study so-called coherence functions [98]. For a beam with destruction operators $b_i$ the second-order coherence function is defined as

The first interesting property of these quantities is that they are not affected by a loss. Assume, for example, that a beam is prepared in certain quantum state $\rho _{\boldsymbol{b}}$ and then undergoes uniform loss by energy transmission $\eta _{\boldsymbol{b}}$ as described in Sec. 3.8. It is straightforward to see that because the numerator and denominator of (291) are normal ordered in the creation and annihilation operators of the beam $\boldsymbol{b}$, then each of them will pick up a factor of $\eta _{\boldsymbol{b}}^2$ that will cancel out. Similarly, for the case of the cross-beam second-order correlation function, if beams $\boldsymbol{b}$ and $\boldsymbol{c}$ undergo uniform loss by energy transmissions $\eta _{\boldsymbol{b}}$ and $\eta _{\boldsymbol{c}}$, respectively, then both the numerator and denominator in the right-hand side of (292) will pick up a factor $\eta _{\boldsymbol{b}} \eta _{\boldsymbol{c}}$ that will cancel out. We emphasize that these results are true irrespective of the state prepared before loss. We now use these results to investigate the correlation functions of degenerate and nondegenerate squeezed light. For nondegenerate squeezed states it is easy to show [99]

The question of calculating individual photon-number probabilities of a multimode Gaussian state has led to the introduction of the Gaussian boson sampling problem [100–102]. In general the evaluation of these probabilities is computationally hard. Moreover, even the computational task of generating samples that distribute according to this probability distribution is believed to be hard (modulo reasonable complexity-theoretic assumptions [103]). The study of this problem is beyond the scope of this tutorial, however we point the reader to a number of recent theoretical [9,97,104] and experimental advances [105,106].

One can also study homodyne measurements where now one is interested in statistics of the quadrature observable:

Considering now nondegenerate squeezing, assuming the LOs for the signal and idler modes are given by the amplitudes $\boldsymbol{\mu }$ and $\boldsymbol{\nu }$ respectively and defining the quadrature observable

4. Heisenberg Equations for Waveguides

In this section, following [107] and [108], we present a general formalism for considering degenerate as well as nondegenerate squeezing in waveguides where either a second- or third-order nonlinearity is dominant. The six nonlinear optical processes that this includes are SPDC, SP-SFWM, and DP-SFWM for the creation of degenerate squeezed vacuum (DSV) states such as those in Sec. 3.5 as well as SPDC, SP-SFWM, and DP-SFWM for the creation of nondegenerate squeezed vacuum (NDSV) states such as those in Sec. 3.6. Beginning with Hamiltonians of the form presented in Sec. 2, suited to integrated photonic structures, we find that the Heisenberg equations of motion, in the limit of strong classical pumps and the undepleted pump approximation, lead to the expected classical coupled mode equations for the pump fields, as well as coupled operator equations for the generated fields. We then discretize coupled operator equations to arrive at compact equations of the form of (207), which can subsequently be solved according to the methods of Sec. 3.4. Perhaps somewhat surprisingly, despite having different physical origins, all of the coupled operator equations can be cast in the same general form.

All of the processes have the same linear Hamiltonian, which we write as (90) plus the GVD term, i.e.,

To simplify the equations resulting from the sum of the Hamiltonians $H^{\text {L}}$ and $H_{\text {NDSV}}^{\text {SP-SFWM}}$ above, we put

With this example in front of us, we note that its coupled operator equations (312) are of the form

In fact, as we show in Appendix  G, whether due to waveguide SPDC, SP-SFWM, or DP-SFWM, coupled operator equations for generating NDSV states can be placed in this general form. Similarly, as shown by Helt and Quesada [107], coupled operator equations for SPDC, SP-SFWM, or DP-SFWM generating DSV states in waveguides can be written in the form

Equations for mean fields, e.g., (311), can be solved using standard approaches, such as split-step Fourier or finite difference methods [109]. The operator equations, on the other hand, can be solved using methods first sketched by Vidrighin [110] and further elaborated by Helt and Quesada [107] and Quesada et al. [108], which we detail in the following.

4.1 Solving the Coupled Operator Equations

Defining

To solve these equations, we discretize $\kappa _{j}=j\Delta \kappa$ and write

Propagation over a finite interval can then be calculated by concatenating infinitesimal propagators $\textbf {K}^{J_{1}J_{2}}(t)$ to form the (two-argument) Heisenberg picture propagator [108]

4.2 Including Loss

With these expressions known, we may now form the phase-sensitive and phase-insensitive moments (cf. Secs. 3.5 and 3.6)

We note that, because the generated state is Gaussian, it is completely characterized by these two moments, and thus any experimentally accessible quantity, such as coherence functions [38,99], can be written in terms of them. Following Sec. 3.8, including loss is then simply a matter of properly updating the moments. In particular, for uniform waveguide loss

4.3 Example: “Separable” JSAs

The formalism presented above is quite general in that it allows for the calculation of squeezed vacuum states in waveguides where either a second- or third-order nonlinearity is dominant, correctly accounting for dispersion to any desired order, SPM and XPM, if present, and loss, in either the so-called low- or high-gain regime. However, it is also something of a black box, in that it relies on numerics. Thus, to help build intuition as one moves from the low- to the high-gain regime, here we provide an illustrative calculation of a lossless waveguide SP-SFWM process engineered to produce a separable JSA in the low-gain regime ignoring GVD, SPM, and XPM under three scenarios: (i) ignoring time-ordering corrections to arrive at an idealized analytic result; (ii) considering leading-order time-ordering corrections via the Magnus expansion; and (iii) using the full formalism presented previously. We note that we neglect GVD, SPM, and XPM in (iii) not because the numerical formalism cannot include them, but because doing so simplifies comparisons against the analytic results.

For completeness, using (79) and (303) we write the linear and nonlinear Hamiltonians as

4.3.1 Low-Gain Regime

In the low-gain regime, a common, and safe, approximation is to neglect time-ordering effects and, instead, using $\int _{-\infty }^{\infty }\textrm{d}{t^{\prime }}\,e^{i\Delta t^{\prime }}=2\pi \delta \left (\Delta \right )$, write

Here we choose to write the Schmidt decomposition of such a JSA as

Note that the Schmidt number of a separable JSA, $J_{SI}(\omega,\omega ')=\xi _0 f_{0;S}(\omega ) f_{0;I}(\omega ')$, is 1, whereas JSAs that are not separable have $K>1$. Further note that we have chosen to work with unnormalized JSAs. Although this marks a departure from more common usage, we find it useful here for comparing the results of this subsection against those that move beyond the approximations used in deriving $J_{1;SI}\left (\omega _{1},\omega _{2}\right )$. In the low-gain regime, where $K\approx \left (\sum _{\lambda } \xi ^{2}_{\lambda }\right )^2/\sum _{\lambda } \xi ^{4}_{\lambda }$, the Schmidt number is insensitive to a scaling transformation of $J_{SI}$. That is, $J_{SI}$ and $l J_{SI}$ have the same Schmidt number for any $l \neq 0$, though this does not hold in general.

Putting

We note, in passing, that the only square integrable function that is separable and has elliptical (or circular) contours is precisely a two-dimensional Gaussian, and that, moreover, only a combination of a Gaussian PMF and a Gaussian pump can give a separable JSA in the low-gain regime [121,122]. Under this condition (352), and introducing

We note that this JSA is a two-dimensional Gaussian with elliptical contours, becoming a two-dimensional Gaussian with circular contours for $\tau _{S}=\tau _{I}=\tau$. Using (352) to show that

4.3.2 Magnus Expansion

Such an analysis describes the state of generated photons well for small enough $\Phi$. However, as $\Phi$ increases, whether through an increase in $\gamma _{\text {chan}}^{SIPP}$, waveguide length $L$, or pump pulse energy via $N_{P}$, time-ordering corrections must be considered for a correct description. Using the Magnus expansion, discussed in Sec. 3.2, one can show that the leading-order corrections rewrite the JSA as [115,123]

The quantities $J_{3;SI}$ and $K_{3;SI}$ have been calculated previously [115,123,124] and we can write them in closed form when (352) holds as

This allows us to write the JSA up to these corrections as

Note that when $\Phi \ll 1$ the matrix $\textbf {L}\approx \left ( \begin {array}{cc} \overline {\xi } & 0\\ 0 & 0\\ \end {array} \right )$ and we recover the results obtained in (354). More interestingly, past this limit we can obtain the new time-ordering corrected Schmidt values and functions and verify that they are no longer linear functions of $\Phi$ and independent of $\Phi$, respectively.

 

Figure 10. Amplitudes of the orthonormal functions $\phi _{0}(\Omega )$ and $\phi _{1}(\Omega )$ needed to expand the time-ordering corrected JSA of (364).

Download Full Size | PDF

As we have written the JSA as quadratic form (the matrix $\textbf {L}$) in the orthonormal basis $\{\phi _{i;J}\}$, obtaining the Schmidt decomposition of $J_{SI}$ is equivalent to obtaining the SVD of the matrix $\textbf {L}$. In particular, because $\textbf {L}$ is Hermitian, its SVD is related to its eigendecomposition and we can rewrite (364) as

Written in this form, it is easy to see that now the Schmidt values are simply the absolute values of the eigenvalues of $\textbf {L}$ [cf. (356)], $\xi _0=|\xi _+|, \ \xi _1=|\xi _-|$ and the Schmidt functions are the eigenfunctions of $\textbf {L}$

Thus, the JSA that started as separable in the low-gain regime does not remain so as $\Phi$ increases sufficiently.

4.3.3 Comparisons

We are now in a position to compare the average number of generated pairs per pump pulse

We plot these, as well as a full numerical solution, as functions of $\Phi$, setting $\tau _{S}=\tau _{I}=\tau$, so that

 

Figure 11. Comparison of the generated mean photon number per pump pulse $\langle n \rangle$ and Schmidt number $K$ for no time-ordering corrections (Magnus 1), time-ordering corrections up to $\Phi ^{3}$ (Magnus 3), and a full numerical solution.

Download Full Size | PDF

 

Figure 12. Comparison of the absolute values of the generated JSAs $\left \vert J_{SI}(\Omega _{1},\Omega _{2})\right \vert$ with $\Phi =4$ (left column) for (a) no time-ordering corrections, (b) time-ordering corrections up to $\Phi ^{3}$, and (c) a full numerical solution, as well as $\Phi =11$ (right column) for (d) no time-ordering corrections, (e) time-ordering corrections up to $\Phi ^{3}$, and (f) a full numerical solution.

Download Full Size | PDF

5. Heisenberg Equations for Rings

In this section we build on top of Sec. 2 to continue the analysis of ring resonators in the presence of second- or third-order nonlinearities. Integrated optical resonators are some of the most efficient and successful sources of quantum states of light, thanks to the large field enhancements and the possibility to tailor independently the density of states of each mode involved in the nonlinear process, thereby controlling the statistical properties of the generated quantum light [125,126]. A wide range of integrated source geometries have been proposed and analyzed, including Fabry–Perot cavities, ring resonators, photonic crystal cavities, and multi-resonator designs, each with a variety of coupling and pumping schemes.

Here, to provide a simple and intuitive example, we restrict the discussion to a simple ring resonator point-coupled to a channel waveguide following the prescriptions of Sec. 2, and sketched in Fig. 13. Similar derivations can be carried out with more complex resonator or coupler designs [128,129], resulting in analogous conclusions.

 

Figure 13. Schematic representation of an integrated ring resonator of radius $R$ point-coupled to a bus waveguide. The interaction is regulated by the self-coupling and cross-coupling coefficients $\sigma$ and $\kappa$, linked by $\sigma ^2+\kappa ^2=1$ [127].

Download Full Size | PDF

When this idealized structure is pumped by bright classical fields, it proves an efficient, stable, and compact source of NDSV [130] or DSV states [131,132]. In the following, we discuss the generation of such states using second- and third-order nonlinearities in the most common pumping schemes, namely SP-SFWM, DP-SFWM, and SPDC (as illustrated in Fig. 14), and present some simplified cases.

 

Figure 14. Energy-level diagrams for the nonlinear processes considered in this section. From left to right: SP-SFWM (Sec. 5.4), SP-SFWM (Sec. 5.5), and SPDC (Sec. 5.6).

Download Full Size | PDF

5.1 SFWM Hamiltonian

Considering a ring with a third-order nonlinearity, first we can extend Sec. 2.4.2 to include SFWM, the driver of squeezing in such a ring. Starting from the form (57) of the $\boldsymbol{D}(\boldsymbol{r})$ operators, we restrict our attention to three resonant modes, with equally spaced wave vectors, conventionally labeled pump, signal, and idler, with frequencies $\omega _{\text {P}}$, $\omega _{\text {S}}$, and $\omega _{\text {I}}$, respectively. Here we have chosen to reduce the scope of this discussion to a three-resonance system, because it greatly simplifies the description of its dynamics. Indeed, with an unconstrained set of resonances, other processes, usually referred to as spurious, and that involve resonant modes beyond the three of interest here can lead to nonlinear frequency shifts and undesired parametric noise being generated completely or partly in the signal mode, affecting the quantum properties of the photons generated by SFWM [133,134]. Although the restriction to a three-resonance system is a practical expedient to simplify our discussion, it should be noted that a variety of interesting and complex quantum states can be generated across more than three resonances, and often extend over entire frequency combs [135–137].

The third-order SFWM contribution to the nonlinear Hamiltonian (141) is

In the expression of the effective area (381), as well as throughout the rest of the section, we assume that dispersion can be neglected within the frequency range of each ring resonance. Therefore, in the field normalization we can truncate the dispersion to the first derivative of the wave vectors. However, considering nonlinear processes that involve multiple resonances, we do not make any further assumptions on the effect of material and waveguide dispersion on the central frequencies of the resonant modes.

5.2 SPDC Hamiltonian

When the ring is characterized by a second-order nonlinear response, and we now consider the three relevant modes conventionally labeled second harmonic (SH), and fundamentals (F1 and F2) with frequencies $\omega _{\text {SH}}$, $\omega _{\text {F1}}$ and $\omega _{\text {F2}}$, respectively, the SPDC contribution to the nonlinear Hamiltonian (142) is

5.3 General Expression of the Dynamic Equation

Before diving into particular scenarios, it is useful to write down the expression of the Heisenberg equation

5.4 SP-SFWM

The resonant wavelengths of the idealized three-mode system introduced above are solutions of $\kappa _{m}\mathcal {L}=\frac {2\pi }{\lambda _{m}}n_{\text {P}}(\omega _{m})\mathcal{L}=2\pi m$ with $m\in \mathbb {Z}$ and are equally spaced in $\kappa$. Solutions are solely affected by the ring geometry, the material and the waveguide dispersion. Let the pump be the central resonance, surrounded by the signal and idler resonances, as sketched in Fig. 15. Although here we consider the interaction between nearest-neighbor resonances, any triple can be considered, as long as equally spaced in $\kappa$.

 

Figure 15. Ring resonances involved in a SP-SFWM process, and the associated SPM and XPM wavelength shifts.

Download Full Size | PDF

In Sec. 2.3 we have already introduced the channel field operators (80), the formalism for the channel–ring coupling (105)–(111), and the related derivation of the full coupling coefficients $\overline {\Gamma }_{J}$ including the scattering loss modeled by a phantom channel (112), (121), and (122). It is worth recalling at this point that such an approach is only suitable for high-finesse systems (131), but this is typically verified for integrated photonic structures at frequencies below the bandgap in most common material platforms such as silicon on insulator (SOI), silicon nitride ($\text {Si}_3\text {N}_4$), lithium niobate on insulator (LNOI), aluminum gallium arsenide (AlGaAs), aluminum nitride (AlN), and chalcogenide glasses.

We illuminate the ring with a single, bright light source, where it is implicit that its bandwidth cannot exceed that of the central resonance. We refer to this case as SP-SFWM. From the full Hamiltonian of the system by composition of the linear Hamiltonian (105) and the nonlinear response due to SPM (161), XPM (168), and SFWM (377), we obtain

In the interest of simplifying notation, throughout the rest of this section we drop the superscripts and subscripts unless necessary. Given the full Hamiltonian (388) we can now calculate the Heisenberg equations of motion for the ring operators following (387), and obtain

Let us then remove the rapidly oscillating part of the fields, by defining $\overline {c}_{J}(t)=e^{i\omega _{J}t}c_{J}(t)$ and $$\overline{\psi}_{J\lt}(0,t)=e^{i\omega_{J}t}\psi_{J\lt}(0,t)$$so that the equations of motion become

To solve Eqs. (391) let us restrict the discussion to the case of a coherent laser pump, thus treated classically. In this instance, the pump mode operator is replaced by its expectation value

We can reformulate the Heisenberg equations of motion (394) for a ring under SP-SFWM in matrix form as

The solution of (395) is exact within the assumptions introduced, and can, in general, be found numerically. One common approach to solve (395) is to use the Green function formalism, hence identifying a matrix $G(t,t')$ that satisfies $G(t_{0},t_{0})=\mathbb {I}$ and

The full solution can be eventually written as

Once the Green function has been calculated, it can be used to derive the $N$ and $M$ moments of the converted photon distribution, defined as

As discussed in Sec. 3.9 the $N$ and $M$ moments are instrumental in calculating key quantities such as the average photon number of the squeezed state, its variance and covariance matrix (286).

5.5 DP-SFWM

In a variety of applications, such as quantum sensing or photonic quantum computing, it is preferable to generate DSV rather than NDSV states; yet, the SP approach described in the previous section (where one neglects the bright squeezed state generated in the P mode), is not suitable for the task. Nonetheless, with the same high-finesse resonator and still restricting our discussion to three isolated modes, one can invert the pumping scheme and readily obtain a DSV state. Let us illuminate the ring with two individual pumps, named $\text {P1}$ and $\text {P2}$, addressing the side resonances as prescribed in the previous section and represented in Fig. 16.

 

Figure 16. Ring resonances involved in a DP-SFWM process, and the associated SPM and XPM wavelength shifts.

Download Full Size | PDF

Incidentally, we note that a two-lobed pump where each lobe is in correspondence of a ring resonance is a completely equivalent view of the pumping scheme. We can conveniently rename the three ring resonances as P1, S, and P2, respectively, where it is assumed that the DSV state is generated in the S resonance by SFWM. We refer to this pumping scheme as DP-SFWM.

The full Hamiltonian of the process is

As before, in the following we drop the unnecessary subscripts for the sake of clarity. We can proceed with the same reasoning as in the SP-SFWM case and work out the Heisenberg equations of motion, to obtain

As usual, we assume the pumps to be coherent lasers, thus treated classically. Each pump mode operator is then replaced by its expectation value

The analog of the Heisenberg equations of motion (395) for DP-SFWM is

We can finally calculate the $N$ and $M$ moments of the photon distribution as

5.5.1 Dynamics of the Cavity and Field Operators for DP-SFWM

So far, our main goal was to express the Heisenberg equations of motion in an analytic form, with little insight on how one could solve them and obtain the full dynamics of the cavity. That, in turn, allows one to deduce relevant physical quantities such as the squeezing level, the average photon numbers in the converted modes, the Schmidt number or the distribution of the Schmidt modes. Unsurprisingly, a general analytical solution of those equations cannot be provided. The pumping scheme adopted, the amplitude and phase profile of the pump(s), the strength of the nonlinear terms involved and the resonances detunings are some of the key parameters that affect the solution in a nontrivial way.

Yet, there exist some particular simplified scenarios where an analytical solution can indeed be written in a closed form. For simplicity, here we focus on the DP-SFWM case, but one can get to similar conclusions for SP-SFWM and nondegenerate SPDC. Furthermore, we imagine that XPM from the pump modes to the signal mode can be neglected and that the pumps are on-resonance (hence, $\Delta _{P1(P2)}=0$). Under these conditions we can rephrase (410) as

Finally, using (110) we obtain

Thus, in the idealized case of infinitely short pump pulses and neglecting both loss and XPM, we can analytically express the input–output relation for the squeezed mode and assess the single-Schmidt mode nature of the generated state [128,129].

 

Figure 17. Amplitude of the “ringdown” function $f_{0}(t)$ in (421). Note that we make explicit the fact that the Heaviside function takes the value 1/2 when evaluated at 0.

Download Full Size | PDF

5.5.2 Example: Squeezing Spectrum of a Nanophotonic Molecule

To illustrate the practical applications of the formalism presented in this section, we now consider the real-case scenario of Ref. [133]. The integrated resonator the authors propose is a “photonic molecule” made of two coupled rings, where the additional cavity serves the purpose of suppressing spurious nonlinear processes, allowing for an effective three-mode modeling of the squeezer, precisely as we have so far assumed in all the previous derivations. In this particular instance, the rings are fabricated in silicon nitride and have a $1.5\times 0.8\ \mu {\rm m}^2$ cross section. From the field distribution calculated by eigenmode simulation and the other geometrical parameters reported in Ref. [133] one can extract a nonlinear parameter $\Lambda \sim 5\ {\rm Hz}$. The Hamiltonian of the system is the same as (404), and therefore the Heisenberg equation of motion for the signal mode is (408).

In Ref. [133], the cavity is pumped with two quasi-CW lasers (pulses of constant peak power with duration much greater than the cavity lifetime) of equal intensity, and the squeezing spectrum is measured by balanced homodyne detection [3,139], that consists of interfering the signal with a LO on a 50–50 splitter and recording the difference photocurrent while varying the phase of the LO. Once the $N$ and $M$ moments are calculated as in Eqs. (413) and (414), it is more convenient to work in the frequency domain by taking their Fourier transform $N_{\text {S}}(\Omega, \Omega ')$ and $M_{\text {SS}}(\Omega, \Omega ')$ (given the quasi-CW nature of the pumps). Then, the squeezing spectrum can be expressed as

 

Figure 18. Squeezing and anti-squeezing spectra (blue and black, respectively) as a function of (a) the pump power and (b) the sideband frequency, according to (422). Reprinted by permission from Macmillan Publishers Ltd: Zhang et al., Nature Commun. 12, 2233 (2021) [133]. Copyright (2021).

Download Full Size | PDF

5.6 Parametric Downconversion

We can consider one last time the ring introduced previously, and focus on the second-order nonlinear response, specifically on SPDC. To this end, we can select three resonances conventionally labeled SH, F1, and F2, where the energy of the former approximately equals the sum of the energies of the latter, as depicted in Fig. 19. We note that, in contrast to the SFWM processes, it is more demanding to guarantee energy and phase matching for SPDC processes, given the large frequency difference between the second harmonic and fundamental modes and the optical dispersion of commonly used $\chi _2$ materials. However, several techniques to circumvent this issue have been developed, including quasi-phase matching via waveguide periodical poling [140–142], modal phase matching to higher-order modes (of either same or opposite polarization) [143,144], and cyclic phase matching [145–147].

 

Figure 19. Ring resonances involved in a nondegenerate SPDC process.

Download Full Size | PDF

Pumping the SH resonance with a bright field, the Hamiltonian governing the system becomes

Once more assuming an undepleted pump and ignoring all contributions from the phantom channel, the pump equation becomes

Finally, the downconverted photons Heisenberg equation can be expressed in matrix form as

It is worth noting that the nonlinear Hamiltonian in (423) drives the creation of a NDSV state [cf. (188)] in modes F1 and F2. However, in the limit $\text {F1}\rightarrow \text {F2}\equiv \text {F}$, the solution transitions to a DSV state in the F mode, and all of the apparatus for the calculation of the Heisenberg equations is unchanged, provided that the SPDC nonlinear coefficient is replaced by

6. Connections between Classical and Quantum Nonlinear Optics

In Sec. 3.7 we presented the evolution of a state different from vacuum in the presence of a quadratic Hamiltonian. We saw that even in the presence of a seed field, the joint amplitude appearing in the spontaneous case can act as a response function of the system. In other words, it seems as if there is a deep connection between spontaneous and stimulated nonlinear processes. Starting from these results, we now discuss this relation from a different point of view, with the aim of clarifying the connection between classical and quantum nonlinear optics, thus between those nonlinear phenomena that can be explained in the framework of a classical electromagnetic theory (e.g., second harmonic generation) and those that require the quantization of the electromagnetic field (e.g., SPDC). On the one hand, such a connection is interesting from a fundamental point of view, as it allows one to separate the features of a nonlinear phenomenon that are intrinsically quantum from those that are not. On the other hand, the link between classical and quantum nonlinear optics can also have important practical consequences. First, many systems that have been designed to enhance classical nonlinear processes can inspire new solutions to control and amplify quantum nonlinear interactions. Second, as we show in this section, one can gain important information about the quantum correlations of nonclassical light generated by parametric fluorescence by studying the same system and nonlinear interaction in a regime in which they can both be described classically, with advantages in terms of speed and accuracy.

Here we restrict our analysis to the generation of squeezed light via SPDC or SFWM. As we have seen, both of these processes result in the generation of photon pairs and are associated with either the second- or third-order nonlinear optical response of the system. The generic name with which SPDC and SFWM are often referred to is parametric fluorescence, which suggests a certain analogy with the spontaneous emission of light, for example, occurring in an atomic system. This analogy is quite strong. Indeed, likewise the emission of a photon by an atomic system can be either spontaneous or stimulated by the radiation field, there exist also for SPDC and SFWM two corresponding nonlinear processes, difference frequency generation (DFG) and stimulated four-wave mixing (FWM), in which the emission of photons pairs is stimulated by the presence of an additional field, usually labeled seed [148]. More interestingly, in the limit of a sufficiently weak seed field, thus neglecting SPM and XPM induced by it, both the stimulated and spontaneous emission of pairs are described by the very same Hamiltonian, and the probabilities of spontaneous and stimulated emission of photon pairs are related. Under these hypotheses, the study of the stimulated processes, i.e., DFG and FWM, allows one to gain information about the quantum properties and the generation rate of pairs that would be emitted spontaneously through SPDC or SFWM in the absence of the seed field [149].

6.1 Low-Gain Regime: Generation of Two-Photon States

The connection between stimulated and spontaneous emission of photon pairs is manifest in the context of the generation of two-photon states, i.e., in the low-gain regime, for which one has a small probability of generating a photon pair per pump pulse or, in the case of a CW pump, when such a probability is small within the pump coherence time. As an example, we consider the case of SPDC and DFG in a waveguide in the limit of a CW pump. In this case the photon pair generation rates for spontaneous and stimulated emission can be written as [148]

Here we assume phase matching between $\omega _P$ and $\omega _P/2$, with $P_P$ and $P_S$ the pump and seed powers, respectively. Finally, $\mathcal {P}$ is a characteristic power, which depends on the material nonlinearity, $\mathcal {A}$ is the nonlinear effective area (see (F4) in Appendix F), which we assume to be constant over the frequency range of interest here, $L$ is the length of the nonlinear interaction region, and $\beta _2$ is the GVD.

By taking the ratio of (433) and (432), one obtains a simple relation between the two rates

 

Figure 20. Ratio between the generated idler photons in spontaneous and stimulated FWM for several ring resonators. Reprinted with permission from [152]. Copyright 2012 Optical Society of America.

Download Full Size | PDF

Liscidini and Sipe [149] went even further and demonstrated that the biphoton wave function describing the spontaneously generated two-photon state acts as the response function characterizing how the seed pulse stimulates the emission of photon pairs in the corresponding classical process. Assuming the same pumping condition of the spontaneous process and a seed pulse sufficiently intense to be treated classically, but sufficiently weak to neglect any effect of SPM and XPM induced by the seed, it follows that

 

Figure 21. JSD (a) measured by SET and (b) calculated in the case of photon pairs generated under pulsed pumping in a microring resonator. Reprinted by permission from Macmillan Publishers Ltd: Grassani et al., Sci. Rep. 6, 23564 (2016) [159]. Copyright 2016.

Download Full Size | PDF

6.2 High-Gain Regime: Characterization of Squeezed Light

Although the vast majority of the results regarding the connection between classical and quantum nonlinear processes have been obtained in the low-gain regime, when only a few photon pairs are generated per pump pulse, a similar link is expected to hold also in the high-gain regime, e.g., in the case of the generation of bright squeezed light. Indeed, when effects associated with SPM and XPM induced by the seed field are negligible, stimulated and spontaneous processes are described by the same nonlinear Hamiltonian. This opens to the possibility of characterizing nonclassical light generated in a nonperturbative regime by exploiting the stimulated emission of photon pairs.

Recently, Triginer et al. [92] proposed that cascaded stimulated emission of photon pairs can be used to measure the strength of squeezed light generated by SPDC. Unlike in SET, where the connection between the biphoton wave function and the stimulated idler is direct (see (440)), in this case the authors used two cascaded DFG processes to have a self-reference field, which is made to interfere with the generated squeezed light. From such an interference, it is possible to gain information about the squeezing strength and phase information, which depends on a nontrivial interplay of SPM and XPM associated with the pump field.

7. World beyond Gaussian States: Tensor Networks for the Perplexed

7.1 Curse of Dimensionality in Many-Body Systems and MPSs

In Secs. 3, 4, and 5 we have introduced methods to describe the quantum states generated in quantum nonlinear optics whenever the Heisenberg equations of motion can be linearized. Under these circumstances we can use the methods from Sec. 3. In a nutshell, whenever we can linearize the equations of motion for the field operators we can always write the output of the spontaneous problem as a multimode squeezed state, which in the low-gain limit corresponds to pair production. Although these methods provide useful tools for many problems, they leave open the question of what can be done beyond linearized dynamics and Gaussian state generation. A naïve approach to this problem is to expand every possible time, position or frequency bin into an orthonormal basis, for example, the Fock basis. Thus, a state over $\ell$ bins would be written as

 

Figure 22. (a) Graphical representation of an arbitrary tensor $R$ with indices $i_0,i_1,\ldots, i_{\ell -1}$. (b) Representation of the same tensor using MPSs with auxiliary indices $\alpha _1,\alpha _2,\ldots,\alpha _{\ell -1}$.

Download Full Size | PDF

If a quantum state can be efficiently represented as an MPS, then one can easily calculate expectation values of local observables such as correlation functions. As recently shown by Yanagimoto et al. [30] one can also extract the density matrix of any Schmidt (or broadband) mode $A_\lambda ^\dagger = \sum _{k=0}^{\ell -1} F_{k,\lambda } a_k^\dagger$ defined as a linear combination of the bare modes $\{a_k \}$.

In the following section we show how to construct a constant bond-dimension MPS representation of a resonant system with an arbitrary Hamiltonian coupled to a waveguide with a linear dispersion relation. Note that very recently it was also shown that MPSs provide an efficient description of quantum states where the nonlinearity is distributed along a dispersive nonlinear waveguide [30]. Finally, note that the results we will derive are a qudit extension of the qubit results derived by Schön et al. [165,166], where they consider “the deterministic generation of entangled multiqubit states by the sequential coupling of an ancillary system to initially uncorrelated qubits” and “characterize all achievable states in terms of classes of matrix-product states”.

7.2 MPSs in Cavity QED

We start by writing a Hamiltonian for an isolated resonance coupled to a waveguide [167]

The nonzero canonical commutation relations of the field and resonance are

Note that $H_{\text {NL}}$ is a nonlinear Hamiltonian that is written in terms of resonator operators and that is possibly time-dependent. We leave this completely general, but it can be taken to be, for example, that of a classical field with frequency $\omega _c \approx 3\omega _0$ coupled to the resonance of interest; for this case we would write the Hamiltonian

We can now go to an interaction picture rotating at frequency $\omega _1$ by transforming the Hamiltonian using the following unitary (see Appendix C):

Note that by going into the interaction picture with respect to $\omega _1$ the nonlinear Hamiltonian is now evaluated at the time evolving operators $e^{i \omega _1 t}a^\dagger$ and $e^{-i \omega _1 t} a$. Finally, we can go to a second interaction picture with the unitary operation

Note that the last term on the right-hand side can be written as

This (interaction-picture) Hamiltonian has a very nice physical interpretation. The evolution of the resonator is given by the first two-terms in (456); then, at a given time $t$, the section of waveguide field located at $x=-v(t-t_0)$ gets to exchange some energy with the resonator by a beam splitter (hopping) interaction mediated by the last term in (456).

Now we are ready to discretize time and space. To this end we introduce $t_j = t_0 + j \Delta t$ and write

This operator satisfies $[\psi _j, \psi _k^\dagger ] = \delta _{j,k}$ and we can now write the time-dependent field-resonance Hamiltonian as

The Trotterized time-evolution operator after discretizing time in slices of size Δt and space into $\ell$ slices can be written as

 

Figure 23. Circuit representation of the unitary evolution operator in (460) in a discretization using four modes. Note that in this figure the first wire represents the resonator mode, whereas the rest of the wires represent the different discretized sections of the waveguide. We use vertical lines with dots at the ends to represent the beam splitter operation $\mathcal {B}_{0,j}$ between modes 0 (the resonator) and the $j$th part of the waveguide.

Download Full Size | PDF

 

Figure 24. Circuit representation of the unitary evolution operator in (460) and for a discretization using four modes. Note that in this figure we apply swap operations after each beam splitter. Thus, after the $i$th swap the resonator mode lies in the $i$th wire and, in particular, at the end of the interaction the resonator mode is represented by the last wire. We use vertical lines with dots at the end to represent the beam splitter operation $\mathcal {B}_{i,i+1}$ between wire $i$ and $i+1$ and similarly we use vertical lines with crosses to represent swap operations between wires $i$ and $i+1$.

Download Full Size | PDF

Introducing the following two-wire unitary

We start by writing

Next we introduce

We can proceed like this two more times, defining

Defining the simplified notation $(T^{[k]})^{x,y}_{w,0_p} \equiv (S^{[k]})_{x,y}^{w}$ for $0<k< \ell = 4$ (recall that $0_p$ is fixed from the initial ket (465)) and also $(T^{[0]})^{x,y}_{0_0,0_1} \equiv (S^{[0]})_{i_0,\alpha _1}$ and $(S^{[4]})_{i_4}^{j_4} = \langle i_4|\alpha _4 \rangle$ we can finally write the MPS form

As stated in the previous subsection, once the MPS form of a state is obtained, most useful properties of the state can be obtained in a straightforward manner. The methods to obtain expectation values from MPS states are presented and derived elsewhere [163,164]. We would just like to point out that standard computer libraries already exist for the calculation of expectation values of local observables from MPSs [169–172].

8. Summary, Conclusions, and Future Outlook

Since the initial experiments on squeezing carried out more than three decades ago, there has been a constant drive to move quantum optics from a collection of bulk optical elements spanning large tabletops to a single chip-scale device no larger than a fingernail. Such a transition not only reduces the system size and power requirements, but it also offers the promise of scalability and an unprecedented control over the device parameters, and the properties of the generated light. Although advances in the fabrication of integrated photonic devices have accelerated this evolution, they have also ushered in an era of “high-gain” quantum optics, for example, enabling single-photon generation via post-selection involving photon-number resolving detection, efficient photon frequency conversion, and state generation beyond the single-photon pair regime. For these and other high-gain applications, a first-order perturbation theory, the traditional tool of choice, is no longer sufficient to describe the nonlinear interaction. A key theoretical issue is that relevant Hamiltonian terms describing the nonlinear interaction do not commute with themselves at different times. In short, this means that weakly transformed photon amplitudes at early times influence their own transformation at later times.

Drawing on earlier works, in this tutorial we have shown how to overcome this difficulty, starting from a rigorous first-principles quantization procedure and employing a Trotter–Suzuki expansion of linearized Heisenberg picture dynamics to reduce the problem to the multiplication of matrices with a size proportional to the number of bosonic modes involved. In particular, in Sec. 2 we detailed the quantization of nonlinear optics in integrated photonic structures, paying particular attention to channel waveguides and ring resonators. In Sec. 3 we presented the tools necessary to calculate the output states of a few nonlinear optical processes, whether in the low-gain or high-gain regimes, and whether the initial quantum state is vacuum or not, resulting in a process that is spontaneous or stimulated, respectively. These tools correctly account for loss and material dispersion and, as we showed, they can be readily applied to calculate the results of homodyne and photon-number resolved measurements of the generated states. With these mathematical tools firmly established, we then turned to the particular structures singled out in Sec. 2, going into further detail regarding their application to waveguides in Sec. 4 and ring resonators in Sec. 5 for various SPDC and SFWM processes. In Sec. 6 we discussed the deep connection between spontaneous and stimulated processes, showing how the knowledge of the latter can be used to inform the former. Finally, in Sec. 7 we examined how tensor networks and MPSs can be used to move beyond linearized equations of motion and Gaussian state generation. Throughout this tutorial, we have also emphasized connections to relevant theoretical and experimental works.

It is our hope that the tools provided here will benefit the entire integrated quantum photonics community, leading to devices with still greater complexity and enhanced functionalities than their present-day counterparts, especially as state generation moves beyond the single or few photon pairs regime. Indeed, integrated quantum photonic technologies promise ultimate information security, enhanced measurement precision, and greater computing power than traditional technologies, and their realization is now within reach.

Appendix A: Field Normalization and Dispersion

As a first example we look at a uniform medium in the dispersionless limit, where $\varepsilon _{1}(\boldsymbol {r})=\varepsilon _{1}$. Then the equations (18) for $\boldsymbol {B}_{\alpha }(\boldsymbol {r})$ become simply

(A1)$$\begin{aligned} & \boldsymbol{\nabla}\times\left(\boldsymbol{\nabla}\times\boldsymbol{B}_{\alpha}(\boldsymbol{r})\right)=\frac{\omega_{\alpha}^{2}\varepsilon_{1}}{c^{2}}\boldsymbol{B}_{\alpha}(\boldsymbol{r}),\\ & \boldsymbol{\nabla}\cdot\boldsymbol{B}_{\alpha}(\boldsymbol{r})=0. \end{aligned}$$

Solutions satisfying periodic boundary conditions are of the form $\boldsymbol{b}\exp(i\boldsymbol{k}\cdot\boldsymbol{r})$, where the vector $\boldsymbol {b}$ is perpendicular to $\boldsymbol {k}$ and

(A2)$$ \omega_{\alpha}=c\left|\boldsymbol{k}\right|/n, $$

where $n=\sqrt {\varepsilon _{1}}$. If the normalization volume $V$ is a cube with each side of length $L$, $V=L^{3}$, the allowed $\boldsymbol {k}$ are of the form

(A3)$$\boldsymbol{k}=\frac{2\pi}{L}(m_{x}\boldsymbol{\hat{x}}+m_{y}\boldsymbol{\hat{y}}+m_{z}\boldsymbol{\hat{z}}), $$

where the $m_{i}$ are integers and $\boldsymbol {\hat {x}}$ is a unit vector pointing in the $x$ direction, etc. There are two independent vectors perpendicular to $\boldsymbol {k}$; call them $\boldsymbol {\hat {e}}_{1\boldsymbol {k}}$ and $\boldsymbol {\hat {e}}_{2\boldsymbol {k}}$, and define them such that

(A4)$$\boldsymbol{\hat{e}}_{1\boldsymbol{k}}\times\boldsymbol{\hat{e}}_{2\boldsymbol{k}}=\boldsymbol{\hat{k}}. $$

In a uniform medium it is often useful to work with left- and right-handed circularly polarized components of light, which can be identified by the polarization vectors

(A5)$$\begin{aligned} & \boldsymbol{\hat{e}}_{L\boldsymbol{k}}={-}\frac{1}{\sqrt{2}}\left(\boldsymbol{\hat{e}}_{1\boldsymbol{k}}+i\boldsymbol{\hat{e}}_{2\boldsymbol{k}}\right),\\ & \boldsymbol{\hat{e}}_{R\boldsymbol{k}}=\frac{1}{\sqrt{2}}\left(\boldsymbol{\hat{e}}_{1\boldsymbol{k}}-i\boldsymbol{\hat{e}}_{2\boldsymbol{k}}\right), \end{aligned}$$

respectively. Using $I$ to denote either $L$ or $R$ (positive or negative helicity respectively), and $\overline {I}$ to equal the opposite helicity, we have

(A6)$$\begin{aligned} & \boldsymbol{\hat{e}}_{I\boldsymbol{k}}\cdot\boldsymbol{\hat{e}}_{I\boldsymbol{k}}^{*}=1,\\ & \boldsymbol{\hat{e}}_{I\boldsymbol{k}}\cdot\boldsymbol{\hat{e}}_{\overline{I}\boldsymbol{k}}^{*}=0, \end{aligned}$$

and

(A7)$$\begin{aligned} & i\boldsymbol{\hat{k}}\times\boldsymbol{\hat{e}}_{I\boldsymbol{k}}=s_{I}\boldsymbol{\hat{e}}_{I\boldsymbol{k}},\\ & \hat{\boldsymbol{e}}_{I\boldsymbol{k}}\times\boldsymbol{\hat{e}}_{I\boldsymbol{k}}^{*}={-}is_{I}\boldsymbol{\hat{k}}, \end{aligned}$$

where $\boldsymbol{\hat{k}}\equiv{\boldsymbol{k}}/\left|\boldsymbol{k}\right|$ and $s_{L}=1$, $s_{R}=-1$ identifies the helicity. It is convenient to link the definitions of $\boldsymbol {\hat {e}}_{1\boldsymbol {k}}$ and $\boldsymbol {\hat {e}}_{2\boldsymbol {k}}$ to those of $\boldsymbol {\hat {e}}_{1(-\boldsymbol {k})}$ and $\boldsymbol {\hat {e}}_{2(-\boldsymbol {k})}$ according to

(A8)$$\begin{aligned} & \boldsymbol{\hat{e}}_{1(-\boldsymbol{k})}=\boldsymbol{\hat{e}}_{1\boldsymbol{k}},\\ & \boldsymbol{\hat{e}}_{2(-\boldsymbol{k})}={-}\boldsymbol{\hat{e}}_{2\boldsymbol{k}}, \end{aligned}$$

for then

(A9)$$\begin{aligned} & \boldsymbol{\hat{e}}_{I\boldsymbol{k}}^{*}=\boldsymbol{\hat{e}}_{I(-\boldsymbol{k})},\\ & \boldsymbol{\hat{e}}_{I\boldsymbol{k}}\cdot\boldsymbol{\hat{e}}_{I'(-\boldsymbol{k})}=\delta_{II'}. \end{aligned}$$

We can label the modes $\alpha$ by the pair of indices $I,\boldsymbol {k}$, and seeking $\boldsymbol {b}$ vectors proportional to $s_{I}\boldsymbol {\hat {e}}_{I\boldsymbol {k}}$ the normalization condition (23) means we can take

(A10)$$\boldsymbol{B}_{I\boldsymbol{k}}(\boldsymbol{r})=\sqrt{\frac{\mu_{0}}{V}}s_{I}\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}},$$

whereas from the last of (18) we then find

(A11)$$\boldsymbol{D}_{I\boldsymbol{k}}(\boldsymbol{r})=in\sqrt{\frac{\epsilon_{0}}{V}}\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}$$

and for the mode frequencies we write

(A12)$$\omega_{I\boldsymbol{k}}=\frac{c}{n}\left|\boldsymbol{k}\right|, $$

which, of course, only depend on $\left |\boldsymbol {k}\right |$. Note that the partner mode of a mode of polarization $I$ with propagation wave vector $\boldsymbol {k}$ is the mode of polarization $I$ with propagation wave vector $-\boldsymbol {k}$. Clearly other choices of the modes in addition to (A10) and (A11) are possible; one could multiply both by an overall phase factor or the factor $s_{I}$, or one could work with the linearly polarized basis $\boldsymbol {\hat {e}}_{1\boldsymbol {k}}$ and $\boldsymbol {\hat {e}}_{2\boldsymbol {k}}$ at each $\boldsymbol {k}$. However, choosing (A10) and (A11) in (38) we find

(A13)$$\begin{aligned} & \boldsymbol{D}(\boldsymbol{r},t)=i\sum_{I,\boldsymbol{k}}\sqrt{\frac{\hbar cn\left|\boldsymbol{k}\right|\epsilon_{0}}{2V}}a_{I\boldsymbol{k}}(t)\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}+\text{H.c.},\\ & \boldsymbol{B}(\boldsymbol{r},t)=\sum_{I,\boldsymbol{k}}\sqrt{\frac{\hbar c\left|\boldsymbol{k}\right|\mu_{0}}{2nV}}a_{I\boldsymbol{k}}(t)s_{I}\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}+\text{H.c.} \end{aligned}$$

Finally, we consider the passage to an infinite normalization volume. Here the analogs of (49) and (51) in three dimensions are

(A14)$$\sum_{\boldsymbol{k}}\rightarrow\frac{\int \textrm{d}\boldsymbol{k}}{\left(\frac{2\pi}{L}\right)^{3}}=\frac{V}{8\pi^{3}}\int {\rm d}\boldsymbol{k}$$

and

(A15)$$a_{I\boldsymbol{k}}(t)\rightarrow\sqrt{\frac{8\pi^{3}}{V}}a_{I}(\boldsymbol{k},t).$$

We then have

(A16)$$\begin{aligned} & \left[a_{I}(\boldsymbol{k},t),a_{I'}^{{\dagger}}(\boldsymbol{k'},t)\right]=\delta_{II'}\delta(\boldsymbol{k}-\boldsymbol{k'}),\\ & \left[a_{I}(\boldsymbol{k},t),a_{I'}(\boldsymbol{k'},t)\right]=0. \end{aligned}$$

In terms of these variables $a_{I}(\boldsymbol {k},t)$ (A13) become

(A17)$$\begin{aligned} & \boldsymbol{D}(\boldsymbol{r},t)=i\sum_{I}\int\sqrt{\frac{\hbar cn\left|\boldsymbol{k}\right|\epsilon_{0}}{16\pi^{3}}}a_{I}(\boldsymbol{k},t)\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}\textrm{d}\boldsymbol{k}+\text{H.c.},\\ & \boldsymbol{B}(\boldsymbol{r},t)=\sum_{I}\int\sqrt{\frac{\hbar c\left|\boldsymbol{k}\right|\mu_{0}}{16\pi^{3}n}}a_{I}(\boldsymbol{k},t)s_{I}\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}\textrm{d}\boldsymbol{k}+\text{H.c.}, \end{aligned}$$

where we keep the same notation for $\boldsymbol {\hat {e}}_{I\boldsymbol {k}}$ and $\omega _{I\boldsymbol {k}}$ as in (A13), although $\boldsymbol {k}$ now ranges continuously. The Hamiltonian (40) becomes

(A18)$$H=\sum_{I}\int\hbar\omega_{I\boldsymbol{k}}a_{I}(\boldsymbol{k},t)a_{I}(\boldsymbol{k,}t)\textrm{d}\boldsymbol{k}.$$

We now turn to the inclusion of dispersion. The consequences of this new factor $v_{p}(\boldsymbol {r};\omega _{\alpha })/v_{g}(\boldsymbol {r};\omega _{\alpha })$ in the normalization condition (67) can be easiest seen by first considering periodic boundary conditions. As in a uniform medium the relative dielectric constant and the velocities are independent of position, we can put $\varepsilon _{1}(\boldsymbol {r};\omega _{I\boldsymbol {k}})\rightarrow \varepsilon _{1}(\omega _{I\boldsymbol {k}})$, $n(\boldsymbol {r},\omega _{I\boldsymbol {k}})\rightarrow n(\omega _{I\boldsymbol {k}})$, $v_{p,g}(\boldsymbol {r};\omega _{I\boldsymbol {k}})\rightarrow v_{p,g}(\omega _{I\boldsymbol {k}})$, and the modes (A10) and (A11) in (38) in the dispersionless limit are replaced by

(A19)$$\boldsymbol{B}_{I\boldsymbol{k}}(\boldsymbol{r})=\sqrt{\frac{\mu_{0}v_{g}(\omega_{I\boldsymbol{k}})}{Vv_{p}(\omega_{I\boldsymbol{k}})}}s_{I} \boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}},$$

and

(A20)$$\boldsymbol{D}_{I\boldsymbol{k}}(\boldsymbol{r})=in(\omega_{I\boldsymbol{k}})\sqrt{\frac{\epsilon_{0}v_{g}(\omega_{I\boldsymbol{k}})}{Vv_{p} (\omega_{I\boldsymbol{k}})}}\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}},$$

where now the mode frequency for a given $\boldsymbol {k}$ is determined by self-consistently solving

(A21)$$n(\omega_{I\boldsymbol{k}})\omega_{I\boldsymbol{k}}=c\left|\boldsymbol{k}\right|$$

for $\omega _{I\boldsymbol {k}}$ (compare (A2)), and the field expansions (A17) become

(A22)$$\begin{aligned} & \boldsymbol{D}(\boldsymbol{r},t)=i\sum_{I,\boldsymbol{k}}\sqrt{\frac{\hbar\left|\boldsymbol{k}\right|\epsilon_{0}\varepsilon_{1} (\omega_{I\boldsymbol{k}})v_{g}(\omega_{I\boldsymbol{k}})}{2V}}a_{I\boldsymbol{k}}(t)\boldsymbol{\hat{e}}_{I\boldsymbol{k}} e^{i\boldsymbol{k}\cdot\boldsymbol{r}}+\text{H.c.},\\ & \boldsymbol{B}(\boldsymbol{r},t)=\sum_{I,\boldsymbol{k}}\sqrt{\frac{\hbar\left|\boldsymbol{k}\right|\mu_{0}v_{g}(\omega_{I\boldsymbol{k}})}{2V}} a_{I\boldsymbol{k}}(t)s_{I}\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}+\text{H.c.} \end{aligned}$$

Again we form the Poynting vector $\boldsymbol {S}(\boldsymbol {r},t)=\boldsymbol {E}(\boldsymbol {r},t)\times \boldsymbol {H}(\boldsymbol {r},t)$; using (63) and (64) we have

(A23)$$\begin{aligned} & \boldsymbol{E}(\boldsymbol{r},t)=i\sum_{I,\boldsymbol{k}}\sqrt{\frac{\hbar\left|\boldsymbol{k}\right|v_{g}(\omega_{I\boldsymbol{k}})}{ 2V\epsilon_{0}\varepsilon_{1}(\omega_{I\left|\boldsymbol{k}\right|})}}a_{I\boldsymbol{k}}(t)\boldsymbol{\hat{e}}_{I\boldsymbol{k}} e^{i\boldsymbol{k}\cdot\boldsymbol{r}}+\text{H.c.},\\ & \boldsymbol{H}(\boldsymbol{r},t)=\sum_{I,\boldsymbol{k}}\sqrt{\frac{\hbar\left|\boldsymbol{k}\right|v_{g}(\omega_{I\boldsymbol{k}})}{2V\mu_{0}}} a_{I\boldsymbol{k}}(t)s_{I}\boldsymbol{\hat{e}}_{I\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}+\text{H.c.}, \end{aligned}$$

and considering a state $\left |\Psi \right \rangle$ with one photon in mode $I,\boldsymbol {k}$,

$$\left|\Psi\right\rangle =a_{I\boldsymbol{k}}^{{\dagger}}(0)\left\vert\text{vac}\right\rangle$$

where $\left |vac\right \rangle$ is the vacuum state, we find

$$\int_{V}\left\langle \boldsymbol{S}(\boldsymbol{r},t)\right\rangle d\boldsymbol{r}=\hbar\omega_{I\boldsymbol{k}}v_{g}(\omega_{I\boldsymbol{k}}), $$

with the group velocity governing the energy flow, as expected.

For reference, we give the field expansion for a uniform infinite medium, the appropriate generalization of (A17):

(A24)$$\begin{aligned} & \boldsymbol{D}(\boldsymbol{r},t)=i\sum_{I}\int\sqrt{\frac{\hbar\left|\boldsymbol{k}\right|\epsilon_{0} \varepsilon_{1}(\omega_{I\boldsymbol{k}})v_{g}(\omega_{I\boldsymbol{k}})}{16\pi^{3}}}a_{I}(\boldsymbol{k},t) \boldsymbol{\hat{e}}_{I}(\boldsymbol{k})e^{i\boldsymbol{k}\cdot\boldsymbol{r}}d\boldsymbol{k}+\text{H.c.},\\ & \boldsymbol{B}(\boldsymbol{r},t)=\sum_{I}\int\sqrt{\frac{\hbar\left|\boldsymbol{k}\right|\mu_{0}v_{g}(\omega_{I\boldsymbol{k}})}{16\pi^{3}}} a_{I}(\boldsymbol{k},t)s_{I}\boldsymbol{\hat{e}}_{I}(\boldsymbol{k})e^{i\boldsymbol{k}\cdot\boldsymbol{r}}d\boldsymbol{k}+\text{H.c.}, \end{aligned}$$

where the dispersion relation is given by  (A21), although $\boldsymbol {k}$ now varies continuously, and the commutation relations  (A16) and the form of the Hamiltonian  (A18) are as they were in the limit of no dispersion.

Appendix B: Group Velocity

In this appendix we confirm the expression (73) for the channel mode group velocity. Temporarily introducing

(B1)$$\begin{aligned} & \mathcal{B}_{Ik}(\boldsymbol{r})=\boldsymbol{b}_{Ik}(y,z)e^{ikx},\\ & \mathcal{D}_{Ik}(\boldsymbol{r})=\boldsymbol{d}_{Ik}(y,z)e^{ikx}, \end{aligned}$$

we begin with the master equation (recall the first of (66)), which we can write as

(B2)$$\boldsymbol{\nabla}\times\left(\frac{\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right) =\frac{\omega_{Ik}^{2}}{c^{2}}\mathcal{B}_{Ik}(\boldsymbol{r}),$$

and take its derivative with respect to $k$. Using

(B3)$$\begin{aligned} & \frac{\partial\mathcal{B}_{Ik}(\boldsymbol{r})}{\partial k}=\frac{\partial\boldsymbol{b}_{Ik}(y,z)}{\partial k}e^{ikx}+ix\mathcal{B}_{Ik}(\boldsymbol{r}),\\ & \frac{\partial}{\partial k}\left(\frac{1}{\varepsilon_{1}(y,z;\omega_{Ik})}\right)=\frac{\partial\omega_{Ik}}{\partial k}\frac{\partial}{\partial\omega_{Ik}}\left(\frac{1}{\varepsilon_{1}(y,z;\omega_{Ik})}\right) \end{aligned}$$

and putting

(B4)$$ v_{Ik}\equiv\frac{\partial\omega_{Ik}}{\partial k}, $$

the group velocity of the mode at $k$ including both modal and material dispersion effects, we find

(B5)$$\begin{aligned} & \boldsymbol{\nabla}\times\left(\boldsymbol{C}_{1}(\boldsymbol{r})+\boldsymbol{C}_{2}(\boldsymbol{r})+v_{Ik}\boldsymbol{C}_{3}(\boldsymbol{r)}\right)\\ & =\frac{2\omega_{Ik}}{c^{2}}v_{Ik}\mathcal{B}_{Ik}(\boldsymbol{r})+\frac{\omega_{Ik}^{2}}{c^{2}}\left(\frac{\partial\boldsymbol{b}_{Ik}(y,z)}{\partial k}e^{ikx}+ix\mathcal{B}_{Ik}(\boldsymbol{r})\right),\end{aligned}$$

where

(B6)$$\begin{aligned} & \boldsymbol{C}_{1}(\boldsymbol{r})=\frac{\boldsymbol{\nabla}\times\left(\frac{\partial\boldsymbol{b}_{Ik}(y,z)}{\partial k}e^{ikx}\right)}{\varepsilon_{1}(y,z;\omega_{Ik})},\\ & \boldsymbol{C}_{2}(\boldsymbol{r})=\frac{\boldsymbol{\nabla}\times\left(ix\mathcal{B}_{Ik}(\boldsymbol{r})\right)}{\varepsilon_{1} (y,z;\omega_{Ik})}=\frac{i\boldsymbol{\hat{x}}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})} +\frac{ix\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})},\\ & \boldsymbol{C}_{3}(\boldsymbol{r})=\frac{\partial}{\partial\omega_{Ik}}\left(\frac{1}{\varepsilon_{1}(y,z;\omega_{Ik})}\right) \boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r}). \end{aligned}$$

Now write the curl of $\boldsymbol {C}_{2}(\boldsymbol {r})$ as

(B7)$$\boldsymbol{\nabla}\times\boldsymbol{C}_{2}(\boldsymbol{r})=\boldsymbol{\nabla}\times\boldsymbol{C}_{4}(\boldsymbol{r}) +ix\boldsymbol{\nabla}\times\left(\frac{\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right) +\boldsymbol{C}_{6}(\boldsymbol{r}),$$

where

(B8)$$\begin{aligned} & \boldsymbol{C}_{4}(\boldsymbol{r})=\frac{i\boldsymbol{\hat{x}}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})},\\ & \boldsymbol{C}_{6}(\boldsymbol{r})=i\hat{\boldsymbol{x}}\times\left(\frac{\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1} (y,z;\omega_{Ik})}\right). \end{aligned}$$

Using the expression  (B7) for $\boldsymbol {\nabla }\times \boldsymbol {C}_{2}(\boldsymbol {r})$ in  (B5), we see that the terms proportional to $x$ cancel out by virtue of the master equation  (B2). Further, we can combine the $\boldsymbol {\nabla }\times \boldsymbol {C}_{1}(\boldsymbol {r})$ term on the left-hand side of  (B5) with the term involving $\partial \boldsymbol {b}_{Ik}(y,z)/\partial k$ on the right-hand side of that equation to write  (B5) as

(B9)$$\begin{aligned} & \left(\mathcal{\overline{M}}-\frac{\omega_{Ik}^{2}}{c^{2}}\right)\left(\frac{\partial\boldsymbol{b}_{Ik}(y,z)}{\partial k}e^{ikx}\right)+\boldsymbol{\nabla}\times\boldsymbol{C}_{4}(\boldsymbol{r})+\boldsymbol{C}_{6}(\boldsymbol{r})\\ & =\frac{2\omega_{Ik}}{c^{2}}v_{Ik}\mathcal{B}_{Ik}(\boldsymbol{r})-v_{Ik}\boldsymbol{\nabla}\times\boldsymbol{C}_{3}(\boldsymbol{r}), \end{aligned}$$

where we have introduced the Hermitian operator $\mathcal {\overline {M}}$,

(B10)$$\mathcal{\overline{M}}\boldsymbol{g}(\boldsymbol{r})\equiv\boldsymbol{\nabla}\times\left(\frac{\boldsymbol{\nabla}\times\boldsymbol{g} (\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right)$$

for any vector function $\boldsymbol {g}(\boldsymbol {r}).$

Now choose a particular $k$, and introduce any length $L$. The derivation presented here can be easily extended to a photonic crystal structure as well, where instead of $\boldsymbol {b}_{Ik}(y,z)$ we would have $\boldsymbol {b}_{Ik}(x;y,z),$ with $\boldsymbol {b}_{Ik}(x;y,z)=\boldsymbol {b}_{Ik}(x+a;y,z)$, where $a$ is the length of the unit cell; in such a case, we would take $L=a$, but in a channel waveguide the choice of $L$ is irrelevant. Then take the scalar product of each side of  (B9) into $\mathcal {B}_{Ik}^{*}(\boldsymbol {r})$ and integrate over all $y$ and $z$, and over $x$ from $-L/2$ to $L/2.$ From the first term on the left-hand side of  (B9), we find

(B11)$$\begin{aligned} & \int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\left(\mathcal{\overline{M}}-\frac{\omega_{Ik}^{2}}{c^{2}}\right)\left(\frac{\partial\boldsymbol{b}_{Ik}(y,z)}{\partial k}e^{ikx}\right)\\ & =\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\left(\left(\mathcal{\overline{M}}-\frac{\omega_{Ik}^{2}}{c^{2}}\right)\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\right)\cdot\left(\frac{\partial\boldsymbol{b}_{Ik}(y,z)}{\partial k}e^{ikx}\right)\\ & =0, \end{aligned}$$

where in the second line we have partially integrated and used the fact that the terms in the integrand are the same at $z=-L/2$ and  $z=L/2$, and in the third line we have used the master Eq. (B2). From the second term on the left-hand side of  (B9) we find

(B12)$$\begin{aligned} & \int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\boldsymbol{\nabla}\times\boldsymbol{C}_{4}(\boldsymbol{r})\\ & =i\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\left(\boldsymbol{\nabla}\times\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\right).\left(\frac{\boldsymbol{\hat{x}}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right)\\ & ={-}\mu_{0}\omega_{Ik}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{D}_{Ik}^{*}(\boldsymbol{r}).\left(\frac{\boldsymbol{\hat{x}}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right)\\ & =\mu_{0}\omega_{Ik}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{\hat{x}}\cdot\left(\frac{\mathcal{D}_{Ik}^{*}(\boldsymbol{r})\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right). \end{aligned}$$

In the second line we have used the result

(B13)$$\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{U}(\boldsymbol{r})\cdot\left(\boldsymbol{\nabla}\times\boldsymbol{V}(\boldsymbol{r})\right)=\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{V}(\boldsymbol{r})\cdot\left(\boldsymbol{\nabla}\times\boldsymbol{U}(\boldsymbol{r})\right),$$

which holds for functions $\boldsymbol {U}(\boldsymbol {r})$ and $\boldsymbol {V}(\boldsymbol {r})$ that vanish at the infinite limits of integration for $\textrm{d}y$ and $\textrm{d}z$, and such that they yield no net contribution from the remaining limit of integration (for $\textrm{d}x$); in the third line we have used the result

(B14)$$\mathcal{D}_{Ik}(\boldsymbol{r})=\frac{i}{\mu_{0}\omega_{Ik}}\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r)},$$

which holds because $\mathcal {D}_{Ik}(\boldsymbol {r})\exp (-i\omega _{Ik}t)$ and $\mathcal {B}_{Ik}(\boldsymbol {r})\exp (-i\omega _{Ik}t)$ satisfy the linear Maxwell equations; in the fourth line we have used

(B15)$$\boldsymbol{u}\cdot\left(\boldsymbol{v}\times\boldsymbol{w}\right)=\boldsymbol{v}\cdot\left(\boldsymbol{w}\times\boldsymbol{u}\right),$$

for any three vectors $\boldsymbol {u}$, $\boldsymbol {v}$, and $\boldsymbol {w}$. Finally, from the last term on the right-hand side of (B9) we find

(B16)$$\begin{aligned} & \int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\boldsymbol{C}_{6}(\boldsymbol{r})\\ & =i\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{\hat{x}}\cdot\left(\left(\frac{\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right) \times\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\right)\\ & =\mu_{0}\omega_{Ik}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{\hat{x}}\cdot\left(\frac{\mathcal{D}_{Ik}(\boldsymbol{r})\times\mathcal{B}_{Ik}^{*}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right), \end{aligned}$$

where in the second line we have used (B15) and in the third line we have used (B14). Thus, from (B9) in all we have

(B17)$$\begin{aligned} & \mu_{0}\omega_{Ik}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{\hat{x}}\cdot\left(\frac{\mathcal{D}_{Ik}(\boldsymbol{r})\times\mathcal{B}_{Ik}^{*}(\boldsymbol{r})+\mathcal{D}_{Ik}^{*} (\boldsymbol{r})\times\mathcal{B}_{Ik}(\boldsymbol{r})}{\varepsilon_{1}(y,z;\omega_{Ik})}\right)\\ & =v_{Ik}\mathcal{T}, \end{aligned}$$

where

(B18)$$\begin{aligned} & \mathcal{T}=\frac{2\omega_{Ik}}{c^{2}}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\mathcal{B}_{Ik}(\boldsymbol{r})\\ & -\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\nabla\times\boldsymbol{C}_{3}(\boldsymbol{r}). \end{aligned}$$

Look first at the first term on the right-hand side. As $\mathcal {D}_{Ik}(\boldsymbol {r})\exp (-i\omega _{Ik}t)$ and $\mathcal {B}_{Ik}(\boldsymbol {r})\exp (-i\omega _{Ik}t)$ satisfy the linear Maxwell equations, from Faraday’s law we have

(B19)$$\mathcal{B}_{Ik}(\boldsymbol{r})=\frac{1}{i\omega_{Ik}}\boldsymbol{\nabla}\times\left(\frac{\mathcal{D}_{Ik}(\boldsymbol{r})}{\epsilon_{0} \varepsilon_{1}(y,z;\omega_{Ik})}\right),$$

so

(B20)$$\begin{aligned} & \frac{2\omega_{Ik}}{c^{2}}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\mathcal{B}_{Ik}(\boldsymbol{r})\\ & =\frac{2}{ic^{2}}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\left(\boldsymbol{\nabla}\times\left(\frac{\mathcal{D}_{Ik} (\boldsymbol{r})}{\epsilon_{0}\varepsilon_{1}(y,z;\omega_{Ik})}\right)\right)\\ & =\frac{2}{ic^{2}}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\left(\boldsymbol{\nabla}\times\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\right)\cdot\left(\frac{\mathcal{D}_{Ik} (\boldsymbol{r})}{\epsilon_{0}\varepsilon_{1}(y,z;\omega_{Ik})}\right)\\ & =\frac{2\mu_{0}\omega_{Ik}}{c^{2}}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\left(\frac{\mathcal{D}_{Ik}^{*}(\boldsymbol{r})\cdot\mathcal{D}_{Ik}(\boldsymbol{r})}{\epsilon_{0}\varepsilon_{1}(y,z;\omega_{Ik})}\right), \end{aligned}$$

where in the second expression on the right-hand side we have used  (B13) and in the last expression we have used  (B14). Looking at the second term on the right-hand side of  (B18) we have

(B21)$$\begin{aligned} & -\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\nabla\times\boldsymbol{C}_{3}(\boldsymbol{r})\\ & ={-}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\nabla\times\left(\frac{\partial}{\partial\omega_{Ik}}\left(\frac{1}{\varepsilon_{1} (y,z;\omega_{Ik})}\right)\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})\right)\\ & ={-}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{\nabla}\times\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\left(\frac{\partial}{\partial\omega_{Ik}} \left(\frac{1}{\varepsilon_{1}(y,z;\omega_{Ik})}\right)\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})\right)\\ & ={-}\mu_{0}^{2}\omega_{Ik}^{2}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\left(\mathcal{D}_{Ik}^{*}(\boldsymbol{r})\cdot\mathcal{D}_{Ik}(\boldsymbol{r})\right)\frac{\partial}{\partial\omega_{Ik}} \left(\frac{1}{\varepsilon_{1}(y,z;\omega_{Ik})}\right)\\ & ={-}\frac{2\mu_{0}\omega_{Ik}}{c^{2}}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\left(\mathcal{D}_{Ik}^{*}(\boldsymbol{r})\cdot\mathcal{D}_{Ik}(\boldsymbol{r})\right)\frac{\omega_{Ik}}{2} \frac{\partial}{\partial\omega_{Ik}}\left(\frac{1}{\epsilon_{0}\varepsilon_{1}(y,z;\omega_{Ik})}\right) \end{aligned}$$

or

(B22)$$\begin{aligned} & -\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\nabla\times\boldsymbol{C}_{3}(\boldsymbol{r})\\ &\quad ={-}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\boldsymbol{\nabla}\times\mathcal{B}_{Ik}^{*}(\boldsymbol{r})\cdot\left(\frac{\partial}{\partial\omega_{Ik}}\left(\frac{1}{\varepsilon_{1} (y,z;\omega_{Ik})}\right)\boldsymbol{\nabla}\times\mathcal{B}_{Ik}(\boldsymbol{r})\right)\\ & ={-}\frac{2\mu_{0}\omega_{Ik}}{c^{2}}\int_{{-}L/2}^{L/2}\textrm{d}x\int \textrm{d}y\textrm{d}z\left(\mathcal{D}_{Ik}^{*}(\boldsymbol{r})\cdot\mathcal{D}_{Ik}(\boldsymbol{r})\right)\frac{\omega_{Ik}}{2}\frac{\partial}{\partial\omega_{Ik}} \left(\frac{1}{\epsilon_{0}\varepsilon_{1}(y,z;\omega_{Ik})}\right), \end{aligned}$$

where in the first term on the right-hand side we have used (B13), and in the second we have used (B14). Combining (B20) and (B22) in (B18) we find

$$\mathcal{T}=\frac{2\mu_{0}\omega_{Ik}}{c^{2}}L,$$

where we have used (B1) and (70). Using this in (B17), and using (B1) and (74) we indeed find (73).

Appendix C: Interaction Picture

In this appendix, we show that it is always possible to go to an interaction picture where the terms proportional to $a_i^\dagger a_j$ in (188) can be eliminated if the coefficients $\Delta _{i,j}$ are time independent. We split the Hamiltonian as

(C1)$$H(t)= \underbrace{ \hbar \sum_{i,j =1}^\ell \Delta_{i,j} a_i^\dagger a_j}_{{\equiv} H_0} + \underbrace{\hbar i \frac12 \sum_{k,l=1}^\ell \left[ \zeta_{k,l}(t) a_k^{{\dagger} } a_{l}^{{\dagger}} -\text{H.c.} \right] }_{{\equiv} V(t)},$$

To study the dynamics of the problem let us derive the equation of motion obeyed by a modified time evolution operator $\mathcal { U}_I$ related to the original time evolution operator $\mathcal { U}$ via

(C2)$$\mathcal{ U}_I(t,t_0)\equiv e^{i H_0 (t-t_0)/\hbar} \mathcal{ U}(t,t_0).$$

Using the Schrödinger equation for $\mathcal { U}$ in (175) we find that

(C3)$$i \hbar \frac{\textrm{d}}{\textrm{d}t} \mathcal{ U}_I(t,t_0)=\left( e^{i H_0(t-t_0)/\hbar} V(t) e^{{-}i H_0(t-t_0)/\hbar} \right) \mathcal{ U}_I(t,t_0)=V_I(t) \mathcal{ U}_I(t,t_0).$$

In the last line we introduced the interaction picture Hamiltonian $V_I(t)$. For the explicit form of our Hamiltonian we can first obtain

(C4)$$e^{i H_0(t-t_0)/\hbar} \boldsymbol{a} e^{{-}i H_0(t-t_0)/\hbar} = \exp\left( i \boldsymbol{\Delta} (t-t_0) \right) \boldsymbol{a} = \boldsymbol{U}(t)\boldsymbol{a},$$

where $\boldsymbol{U}(t)$ is a time-dependent unitary matrix. With this solution we can write

(C5)$$\begin{aligned}V_I(t) = e^{i H_0(t-t_0)/\hbar} V(t) e^{{-}i H_0(t-t_0)/\hbar} & =\hbar i \frac12 \sum_{k,l=1}^\ell \left[ \zeta_{k,l}(t) e^{i H_0(t-t_0)/\hbar} a_k^{{\dagger} } a_{l}^{{\dagger}} e^{{-}i H_0(t-t_0)/\hbar} -\text{H.c.} \right] \end{aligned}$$

(C6)$$\begin{aligned}& =\hbar i \frac12 \sum_{i,j,k,l=1}^\ell \left[ \zeta_{k,l}(t) U_{l,i}^*(t) a_i^\dagger U_{k,j}^*(t)a_j^\dagger{-}\text{H.c.} \right] \end{aligned}$$

(C7)$$\begin{aligned}& =\hbar i \frac12 \sum_{i,j=1}^\ell \left[ \tilde{\zeta}_{i,j}(t) a_i^\dagger a_j^\dagger{-}\text{H.c.} \right], \end{aligned}$$

where we introduced the symmetric matrix $\tilde {\boldsymbol{\zeta }}(t) = \boldsymbol{U}^\dagger (t) \boldsymbol{\zeta }(t) \boldsymbol{U}^*(t)$. This completes the proof, because now we need to solve for the interaction picture evolution operator which is solely driven by terms of the form $a_i^\dagger a_j^\dagger$.

Appendix D: Symplectic Group

In Sec. 3.4 we showed that the solution of the Heisenberg equations of motion had the form (209), and respected the equal time commutation relations. We can write the canonical commutation relation of the ladder operators more succinctly as

(D1)$$[z_i, z_j^\dagger] = Z_{i,j}, \quad \boldsymbol{z} = \begin{pmatrix} \boldsymbol{a} \\ \boldsymbol{a}^\dagger \end{pmatrix}, \quad \boldsymbol{Z} = \begin{pmatrix} \mathbb{I}_\ell & \boldsymbol{0} \\ \boldsymbol{0} & -\mathbb{I}_\ell \end{pmatrix},$$

and then the fact that the Heisenberg transformation induced by $\boldsymbol{K}$ respects these commutation relations is equivalent to

(D2)$$\boldsymbol{K} \boldsymbol{Z} \boldsymbol{K}^\dagger{=} \boldsymbol{Z}.$$

To connect this equation to the usual definition of the symplectic group one introduces quadratures

(D3)$$\begin{aligned} \boldsymbol{r} \equiv \begin{pmatrix} \boldsymbol{q} \\ \boldsymbol{p} \end{pmatrix} = \sqrt{\hbar}\frac{1}{\sqrt{2}} \begin{pmatrix} \mathbb{I}_\ell & \mathbb{I}_\ell \\ -i \mathbb{I}_\ell & i \mathbb{I}_\ell \end{pmatrix} \begin{pmatrix} \boldsymbol{a} \\ \boldsymbol{a}^\dagger \end{pmatrix} = \sqrt{\hbar}\boldsymbol{R}\boldsymbol{z}. \end{aligned}$$

The canonical commutation relation for the operators in $\boldsymbol{r}$ follows from that in $\boldsymbol{z}$ as

(D4)$$[r_i,r_j] = i \hbar \Omega_{i,j} \quad \boldsymbol{\Omega} = \begin{pmatrix} 0 & \mathbb{I}_\ell \\ -\mathbb{I}_\ell & \boldsymbol{0} \end{pmatrix}.$$

Now it is easy to state the transformation for quadratures

(D5)$$\boldsymbol{z} \to \boldsymbol{z'} = \boldsymbol{K} \boldsymbol{z} \Longleftrightarrow \boldsymbol{r} \to \boldsymbol{r'} = \boldsymbol{S} \boldsymbol{r},$$

where

(D6)$$\begin{aligned} \boldsymbol{S} = \boldsymbol{R} \boldsymbol{K} \boldsymbol{R}^\dagger{=} \begin{pmatrix} \Re(\boldsymbol{V}) + \Re(\boldsymbol{W}) & \Im(\boldsymbol{W}) -\Im(\boldsymbol{V}) \\ \Im(\boldsymbol{V}) + \Im(\boldsymbol{W}) & \Re(\boldsymbol{V}) - \Re(\boldsymbol{W}) \end{pmatrix}. \end{aligned}$$

The condition in  D2 implies that the matrix $\boldsymbol{S}$ needs to satisfy

(D7)$$\boldsymbol{S} \boldsymbol{\Omega} \boldsymbol{S}^T = \boldsymbol{\Omega},$$

which implies that $\boldsymbol{S}$ is a symplectic matrix [5,64–67].

Appendix E: Form of the Evolution Operator in the High-Gain Regime

At the end of Sec. 3.5 we argued that the evolution operator associated with a general quadratic Hamiltonian can be written as the product of two exponentials

(E1)$$\mathcal{U} = \mathcal{U}_1 \mathcal{U}_0, \quad \mathcal{U}_1 = \exp\left[ \frac12 \sum_{k,l=1}^\ell J_{k,l} a^\dagger_{k} a^\dagger_l - \text{H.c.}\right],$$

where $\mathcal {U}_0$ is a unitary operator that satisfies $\mathcal {U}_0 | {\text {vac}}\rangle = | {\text {vac}}\rangle$. In this appendix we determine the form of the generator of $\mathcal {U}_0$. To pin down $\mathcal {U}_0$ we calculate

(E2)$$\mathcal{U}_0^\dagger \mathcal{U}_1^\dagger \boldsymbol{a} \mathcal{U}_1 \mathcal{U}_0 = \boldsymbol{F} \left[ \oplus_{\lambda=1}^\ell \cosh r_\lambda \right] \boldsymbol{F}^* \left\{ \mathcal{U}_0^\dagger \boldsymbol{a} \mathcal{U}_0 \right\} + \boldsymbol{F} \left[ \oplus_{\lambda=1}^\ell \sinh r_\lambda \right] \boldsymbol{F} \left\{ \mathcal{U}_0^\dagger \boldsymbol{a}^\dagger \mathcal{U}_0 \right\},$$

and compare with the solution in (216) which tells us that we want

(E3)$$\boldsymbol{F}^* \left\{ \mathcal{U}_0^\dagger \boldsymbol{a} \mathcal{U}_0 \right\} = \boldsymbol{G} \boldsymbol{a} \Longleftrightarrow \boldsymbol{G}^\dagger\boldsymbol{F}^* \left\{ \mathcal{U}_0^\dagger \boldsymbol{a} \mathcal{U}_0 \right\} = \boldsymbol{a},$$

To satisfy this constraint we recall the results from Appendix C and write

(E4)$$\mathcal{U}_0 = \exp\left[{-}i \sum_{j,k} \phi_{j,k}a^\dagger_j a_k \right],$$

such that

(E5)$$\mathcal{U}_0^\dagger \boldsymbol{a} \mathcal{U}_0 = \exp\left[ - i\boldsymbol{\phi} \right] \boldsymbol{a},$$

and comparing the two equations above tells us that we want $\boldsymbol{G}^\dagger \boldsymbol{F}^* \exp \left [ - i\boldsymbol{\phi } \right ] = \mathbb {I}_{\ell }$ or, equivalently, $\boldsymbol{\phi } = -i \log _e \boldsymbol{G}^\dagger \boldsymbol{F}^*$. Note that since both $\boldsymbol{G}^\dagger$ and $\boldsymbol{F}$ are unitaries their logarithm is $i$ times a Hermitian matrix which upon multiplication by $(-i)$ gives that $\boldsymbol{\phi }$ is a purely Hermitian matrix as it should. Also note that as expected $\mathcal {U}_0$ in (E4) has the vacuum as an eigenket with eigenvalue 1. It is interesting to note that the Heisenberg transformation generated by a general quadratic Hamiltonian requires the exponentiation of two generators. This happens because the exponential map from the Lie algebra $\textbf {sp}(2n, \mathbb {R})$ to the symplectic group Sp$(2n,\mathbb {R})$ is not surjective. However, any element of the group may be generated by the group multiplication of two elements of the group.

Appendix F: Waveguide Hamiltonians

Concerning waveguides in which a second-order nonlinearity is dominant, the process with $2k_{S}=k_{P}$ leads to

(F1)$$H_{\text{DSV}}^{\text{SPDC}}={-}\frac{\zeta_{\text{chan}}^{SSP}}{2}\int\textrm{d}x\,\psi_{S}^{{\dagger}}\left(x\right)\psi_{S}^{{\dagger}}\left(x\right)\psi_{P}\left(x\right)+\textrm{H.c.},$$

and the process with $k_{S}+k_{I}=k_{P}$ leads to

(F2)$$H_{\text{NDSV}}^{\text{SPDC}}={-}\zeta_{\text{chan}}^{SIP}\int\textrm{d}x\,\psi_{S}^{{\dagger}}\left(x\right)\psi_{I}^{{\dagger}}\left(x\right)\psi_{P}\left(x\right)+\textrm{H.c.},$$

where

(F3)$$\begin{aligned}\zeta_{\text{chan}}^{J_{1}J_{2}J_{3}} & =\frac{2}{\epsilon_{0}}\sqrt{\frac{\hbar^{3}\omega_{J_{1}}\omega_{J_{2}}\omega_{J_{3}}}{2^{3}}}\int\textrm{d}y\textrm{d}z\,\Gamma_{2}^{ijk}\left(\boldsymbol{r}\right)\left[d_{J_{1}k_{J_{1}}}^{i}\left(y,z\right)d_{J_{2}k_{J_{2}}}^{j}\left(y,z\right)\right]^{*} d_{J_{3}k_{J_{3}}}^{k}\left(y,z\right)\\ & =\hbar\sqrt{\frac{\hbar\omega_{J_{1}}\omega_{J_{2}}\omega_{J_{3}}}{2\epsilon_{0}\overline{n}_{J_{1}}\overline{n}_{J_{2}}\overline{n}_{J_{3}} A_{\text{chan}}^{J_{1}J_{2}J_{3}}}}\overline{\chi}_{2}, \end{aligned}$$

with

(F4)$$\frac{1}{\sqrt{A_{\text{chan}}^{J_{1}J_{2}J_{3}}}}=\frac{\int\textrm{d}y\textrm{d}z\frac{\chi_{2}^{ijk}\left(y,z\right)}{\overline{\chi}_{2}}\left[e_{J_{1}}^{i}\left(y,z\right)e_{J_{2}}^{j}\left(y,z\right)\right]^{*}e_{J_{3}}^{k}\left(y,z\right)}{\mathcal{N}_{J_{1}}\mathcal{N}_{J_{2}}\mathcal{N}_{J_{3}}\sqrt{c^{3}/\left(v_{J_{1}}v_{J_{2}}v_{J_{3}}\right)}},$$

and we have used (94), (138), and (307).

Concerning waveguides in which a third-order nonlinearity is dominant, the SP DSV process leads to

(F5)$$H_{\text{DSV}}^{\text{SP-SFWM}}={-}\frac{\gamma_{\text{chan}}^{PPPP}}{2}\hbar^{2}\omega_{P}v_{P}^{2}\int\textrm{d}x\,\psi_{P}^{{\dagger}}\left(x\right)\psi_{P}^{{\dagger}}\left(x\right)\psi_{P}\left(x\right)\psi_{P}\left(x\right),$$

the process with $k_{S}+k_{I}=2k_{P}$ leads to (303), the process with $2k_{S}=k_{P_{1}}+k_{P_{2}}$ leads to

(F6)$$\begin{aligned} H_{\text{DSV}}^{\text{DP-SFWM}} & ={-}\gamma_{\text{chan}}^{SSP_{1}P_{2}}\hbar^{2}\overline{\omega}_{SSP_{1}P_{2}}\overline{v}_{SSP_{1}P_{2}}^{2}\int\textrm{d}x\,\psi_{S}^{{\dagger}}\left(x\right)\psi_{S}^{{\dagger}}\left(x\right)\psi_{P_{1}}\left(x\right)\psi_{P_{2}}\left(x\right)+\text{H.c.}\\ & \quad-\frac{\hbar^{2}}{2}\sum_{J}\gamma_{\text{chan}}^{JJJJ}\omega_{J}v_{J}^{2}\int\textrm{d}x\,\psi_{J}^{{\dagger}}\left(x\right)\psi_{J}^{{\dagger}}\left(x\right)\psi_{J}\left(x\right)\psi_{J}\left(x\right)\\ & \quad-2\hbar^{2}\sum_{J\neq S}\gamma_{\text{chan}}^{JSJS}\sqrt{\omega_{J}\omega_{S}}v_{J}v_{S}\int\textrm{d}x\,\psi_{J}^{{\dagger}}\left(x\right)\psi_{S}^{{\dagger}}\left(x\right)\psi_{J}\left(x\right)\psi_{S}\left(x\right)\\ & \quad-2\gamma_{\text{chan}}^{P_{1}P_{2}P_{1}P_{2}}\hbar^{2}\sqrt{\omega_{P_{1}}\omega_{P_{2}}}v_{P_{1}}v_{P_{2}}\int\textrm{d}x\,\psi_{P_{1}}^{{\dagger}}\left(x\right)\psi_{P_{2}}^{{\dagger}}\left(x\right)\psi_{P_{1}}\left(x\right)\psi_{P_{2}}\left(x\right), \end{aligned}$$

and the process with $k_{S}+k_{I}=k_{P_{1}}+k_{P_{2}}$ leads to

(F7)$$\begin{aligned} H_{\text{NDSV}}^{\text{DP-SFWM}} & ={-}2\gamma_{\text{chan}}^{SIP_{1}P_{2}}\hbar^{2}\overline{\omega}_{SIP_{1}P_{2}}\overline{v}_{SIP_{1}P_{2}}^{2}\int\textrm{d}x\,\psi_{S}^{{\dagger}}\left(x\right)\psi_{I}^{{\dagger}}\left(x\right)\psi_{P_{1}}\left(x\right)\psi_{P_{2}}\left(x\right)+\textrm{H.c.}\\ & \quad-\frac{\hbar^{2}}{2}\sum_{J}\gamma_{\text{chan}}^{JJJJ}\omega_{J}v_{J}^{2}\int\textrm{d}x\,\psi_{J}^{{\dagger}}\left(x\right)\psi_{J}^{{\dagger}}\left(x\right)\psi_{J}\left(x\right)\psi_{J}\left(x\right)\\ & \quad-2\hbar^{2}\sum_{J\neq S}\gamma_{\text{chan}}^{JSJS}\sqrt{\omega_{J}\omega_{S}}v_{J}v_{S}\int\textrm{d}x\,\psi_{J}^{{\dagger}}\left(x\right)\psi_{S}^{{\dagger}}\left(x\right)\psi_{J}\left(x\right)\psi_{S}\left(x\right)\\ & \quad-2\hbar^{2}\sum_{J\neq I}\gamma_{\text{chan}}^{JIJI}\sqrt{\omega_{J}\omega_{I}}v_{J}v_{I}\int\textrm{d}x\,\psi_{J}^{{\dagger}}\left(x\right)\psi_{I}^{{\dagger}}\left(x\right)\psi_{J}\left(x\right)\psi_{I}\left(x\right)\\ & \quad-2\gamma_{\text{chan}}^{P_{1}P_{2}P_{1}P_{2}}\hbar^{2}\sqrt{\omega_{P_{1}}\omega_{P_{2}}}v_{P_{1}}v_{P_{2}}\int\textrm{d}x\,\psi_{P_{1}}^{{\dagger}}\left(x\right)\psi_{P_{2}}^{{\dagger}}\left(x\right)\psi_{P_{1}}\left(x\right)\psi_{P_{2}}\left(x\right), \end{aligned}$$

where we have used (304) and (306).

Appendix G: Heisenberg Equations of Motion for Waveguides

Concerning waveguides in which a second-order nonlinearity is dominant, and using (310) along with $2\omega _{S}=\omega _{P}$, the Heisenberg equations of motion yield

(G1)$$\begin{aligned} \left(\frac{\partial}{\partial t}+v_{P}\frac{\partial}{\partial x}-i\frac{v_{P}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle & =0,\\ \left(\frac{\partial}{\partial t}+v_{S}\frac{\partial}{\partial x}-i\frac{v_{S}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\overline{\psi}_{S}\left(x,t\right) & =\frac{i}{\hbar}\zeta_{\text{chan}}^{SSP}\left(x\right)\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle \overline{\psi}_{S}^{{\dagger}}\left(x,t\right), \end{aligned}$$

for the process with $2k_{S}=k_{P}$, and, using $\omega _{S}+\omega _{I}=\omega _{P}$

(G2)$$\begin{aligned} \left(\frac{\partial}{\partial t}+v_{P}\frac{\partial}{\partial x}-i\frac{v_{P}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle & =0,\\ \left(\frac{\partial}{\partial t}+v_{S}\frac{\partial}{\partial x}-i\frac{v_{S}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\overline{\psi}_{S}\left(x,t\right) & =\frac{i}{\hbar}\zeta_{\text{chan}}^{SIP}\left(x\right)\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle \overline{\psi}_{I}^{{\dagger}}\left(x,t\right),\\ \left(\frac{\partial}{\partial t}+v_{I}\frac{\partial}{\partial x}-i\frac{v_{I}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\overline{\psi}_{I}\left(x,t\right) & =\frac{i}{\hbar}\zeta_{\text{chan}}^{SIP}\left(x\right)\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle \overline{\psi}_{S}^{{\dagger}}\left(x,t\right), \end{aligned}$$

for the process with $k_{S}+k_{I}=k_{P}$. Note how the second of (G1) is of the form of (315), and the second and third of (G2) are of the form of (314).

Concerning waveguides in which a third-order nonlinearity is dominant, and again using (310), the Heisenberg equations of motion yield

(G3)$$\begin{aligned} \left(\frac{\partial}{\partial t}+v_{P}\frac{\partial}{\partial x}-i\frac{v_{P}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle & =i\gamma_{\text{chan}}^{PPPP}\hbar\omega_{P}v_{P}^{2}\left|\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle,\\ \left(\frac{\partial}{\partial t}+v_{P}\frac{\partial}{\partial x}-i\frac{v_{P}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\delta\overline{\psi}_{P}\left(x,t\right) & =i\gamma_{\text{chan}}^{PPPP}\hbar\omega_{P}v_{P}^{2}\left[\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle ^{2}\delta\overline{\psi}_{P}^{{\dagger}}\left(x,t\right)\right.\\&\left.+2\left|\left\langle \overline{\psi}_{P}\left(x,t\right)\right\rangle \right|^{2}\delta\overline{\psi}_{P}\left(x,t\right)\right], \end{aligned}$$

where we have put $\overline {\psi }_{P}\left (x,t\right )\rightarrow \left \langle \overline {\psi }_{P}\left (x,t\right )\right \rangle +\delta \overline {\psi }_{P}\left (x,t\right )$, for the SP DSV process. Similarly, using $\omega _{S}+\omega _{I}=2\omega _{P}$, they yield (311) and (312) for the process with $k_{S}+k_{I}=2k_{P}$, using $2\omega _{S}=\omega _{P_{1}}+\omega _{P_{2}}$ they yield

(G4)$$\begin{aligned} &\left(\frac{\partial}{\partial t}+v_{P_{1}}\frac{\partial}{\partial x}-i\frac{v_{P_{1}}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \\&\quad =i\gamma_{\text{chan}}^{P_{1}P_{1}P_{1}P_{1}}\hbar\omega_{P_{1}}v_{P_{1}}^{2}\left|\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle\\ & \qquad+2i\gamma_{\text{chan}}^{P_{1}P_{2}P_{1}P_{2}}\hbar\sqrt{\omega_{P_{1}}\omega_{P_{2}}}v_{P_{1}}v_{P_{2}}\left|\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle,\\ &\left(\frac{\partial}{\partial t}+v_{P_{2}}\frac{\partial}{\partial x}-i\frac{v_{P_{2}}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \\& \quad=i\gamma_{\text{chan}}^{P_{2}P_{2}P_{2}P_{2}}\hbar\omega_{P_{2}}v_{P_{2}}^{2}\left|\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle\\ & \qquad+2i\gamma_{\text{chan}}^{P_{1}P_{2}P_{1}P_{2}}\hbar\sqrt{\omega_{P_{1}}\omega_{P_{2}}}v_{P_{1}}v_{P_{2}}\left|\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle,\\ &\left(\frac{\partial}{\partial t}+v_{S}\frac{\partial}{\partial x}-i\frac{v_{S}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\overline{\psi}_{S}\left(x,t\right)\\ &\quad =2i\gamma_{\text{chan}}^{SSP_{1}P_{2}}\hbar\overline{\omega}_{SSP_{1}P_{2}}\bar{v}_{SSP_{1}P_{2}}^{2}\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \overline{\psi}_{S}^{{\dagger}}\left(x,t\right)\\ & \qquad+2i\gamma_{\text{chan}}^{P_{1}SP_{1}S}\hbar\sqrt{\omega_{P_{1}}\omega_{S}}v_{S}v_{P_{1}}\left|\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \right|^{2}\overline{\psi}_{S}\left(x,t\right)\\ & \qquad+2i\gamma_{\text{chan}}^{P_{2}SP_{2}S}\hbar\sqrt{\omega_{P_{2}}\omega_{S}}v_{S}v_{P_{2}}\left|\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \right|^{2}\overline{\psi}_{S}\left(x,t\right), \end{aligned}$$

for the process with $2k_{S}=k_{P_{1}}+k_{P_{2}}$, and using $\omega _{S}+\omega _{I}=\omega _{P_{1}}+\omega _{P_{2}}$ they yield

(G5)$$\begin{aligned} &\left(\frac{\partial}{\partial t}+v_{P_{1}}\frac{\partial}{\partial x}-i\frac{v_{P_{1}}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \\&\quad =i\gamma_{\text{chan}}^{P_{1}P_{1}P_{1}P_{1}}\hbar\omega_{P_{1}}v_{P_{1}}^{2}\left|\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle\\ & \qquad+2i\gamma_{\text{chan}}^{P_{1}P_{2}P_{1}P_{2}}\hbar\overline{\omega}_{P_{1}P_{2}P_{1}P_{2}}v_{P_{1}}v_{P_{2}}\left|\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle,\\ &\left(\frac{\partial}{\partial t}+v_{P_{2}}\frac{\partial}{\partial x}-i\frac{v_{P_{2}}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \\&\quad =i\gamma_{\text{chan}}^{P_{2}P_{2}P_{2}P_{2}}\hbar\omega_{P_{2}}v_{P_{2}}^{2}\left|\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle\\ & \qquad+2i\gamma_{\text{chan}}^{P_{1}P_{2}P_{1}P_{2}}\hbar\overline{\omega}_{P_{1}P_{2}P_{1}P_{2}}v_{P_{1}}v_{P_{2}}\left|\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \right|^{2}\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle,\\ &\left(\frac{\partial}{\partial t}+v_{S}\frac{\partial}{\partial x}-i\frac{v_{S}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\overline{\psi}_{S}\left(x,t\right)\\ &\quad =2i\gamma_{\text{chan}}^{SIP_{1}P_{2}}\hbar\overline{\omega}_{SIP_{1}P_{2}}\bar{v}_{SIP_{1}P_{2}}^{2}\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \overline{\psi}_{I}^{{\dagger}}\left(x,t\right)\\ & \qquad+2i\gamma_{\text{chan}}^{SP_{1}SP_{1}}\hbar\overline{\omega}_{SP_{1}SP_{1}}v_{S}v_{P_{1}}\left|\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \right|^{2}\overline{\psi}_{S}\left(x,t\right)\\ & \qquad+2i\gamma_{\text{chan}}^{SP_{2}SP_{2}}\hbar\overline{\omega}_{SP_{2}SP_{2}}v_{S}v_{P_{2}}\left|\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \right|^{2}\overline{\psi}_{S}\left(x,t\right),\\ &\left(\frac{\partial}{\partial t}+v_{I}\frac{\partial}{\partial x}-i\frac{v_{I}^{\prime}}{2}\frac{\partial^{2}}{\partial x^{2}}\right)\overline{\psi}_{I}\left(x,t\right)\\ &\quad =2i\gamma_{\text{chan}}^{SIP_{1}P_{2}}\hbar\overline{\omega}_{SIP_{1}P_{2}}\bar{v}_{SIP_{1}P_{2}}^{2}\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \overline{\psi}_{S}^{{\dagger}}\left(x,t\right)\\ & \qquad+2i\gamma_{\text{chan}}^{P_{1}IP_{1}I}\hbar\overline{\omega}_{P_{1}IP_{1}I}v_{P_{1}}v_{I}\left|\left\langle \overline{\psi}_{P_{1}}\left(x,t\right)\right\rangle \right|^{2}\overline{\psi}_{I}\left(x,t\right)\\ & \qquad+2i\gamma_{\text{chan}}^{P_{2}IP_{2}I}\hbar\overline{\omega}_{SP_{2}SP_{2}}v_{P_{2}}v_{I}\left|\left\langle \overline{\psi}_{P_{2}}\left(x,t\right)\right\rangle \right|^{2}\overline{\psi}_{I}\left(x,t\right), \end{aligned}$$

for the process with $k_{S}+k_{I}=k_{P_{1}}+k_{P_{2}}$. Note how the second of  (G3) and the third of  (G4) are both of the form of (315), and the third and fourth of  (G5) are of the form of (314).

Appendix H: Obtaining the Matrix Product Form of an Arbitrary Tensor

Consider an arbitrary tensor with $\ell$ indices $R_{i_0\ldots i_{\ell -1}} \in \mathbb {C}^{c_0 \times c_1 \times \cdots \times c_{\ell -1}}$ where the index $i_k$ can take $c_k$ values. In this appendix we show how one can construct its matrix product representation by leveraging the SVD of matrices, which we now recall. For a $c_0 \times c_1$ rectangular matrix with entries $M_{i j}$ one can construct its SVD as

(H1)$$M_{ij} = \sum_{k = 0}^{\min\{c_0,c_1\}-1} U_{ik} \lambda_{kk} V_{kj}.$$

In the last equation one can take $\boldsymbol{U}$ to be a $c_0 \times c_0$ unitary matrix, $\boldsymbol{V}$ to be $c_1 \times c_1$ unitary matrix and $\boldsymbol{\lambda }$ to be a $c_0 \times c_1$ diagonal matrix with nonnegative entries. Equivalently, one can take $\boldsymbol{U}$ to be a rectangular $c_0 \times \min \{c_0,c_1\}$ semi-unitary matrix, $\boldsymbol{V}$ to be $\min \{c_0,c_1\} \times D_1$ semi-unitary matrix and $\boldsymbol{\lambda }$ to be $\min \{c_0,c_1\} \times \min \{c_0,c_1\}$ square nonnegative matrix.

A second useful property of a multidimensional array is that we are free to reshape it as a matrix by grouping indices together. Consider, for example, a four-index tensor $R_{i_0 i_1 i_2 i_3}$ where each index takes two possible values, $i_k \in \{0,1\}$. We can group the last three indices $i_1 i_2 i_3$ into a common index that takes $2^3 = 8$ values. This allows us to consider the tensor $R$ as a $2 \times 8$ matrix. We can emphasize this partition by putting a semicolon between the indices we want to separate writing

(H2)$$R_{i_0;i_1 i_2 i_3} \in \mathbb{C}^{2 \times 8}.$$

We can instead group the first and last two indices to now obtain a $2^2 \times 2^2 = 4 \times 4$ matrix that we write as

(H3)$$R_{i_0 i_1 ; i_2 i_3} \in \mathbb{C}^{4 \times 4}.$$

Having introduced this useful notation and recalled the SVD for matrices we are now ready to decompose an arbitrary tensor.

To this end, we split our arbitrary input tensor as $R_{i_0;i_1\ldots i_{\ell -1}}$ and then write the singular decomposition for this $c_0 \times c_{1:\ell -1}$ (with $c_{1:\ell -1} = c_1 \times \cdots \times c_{\ell -1}$) matrix

(H4)$$R_{i_0;i_1\ldots i_{\ell-1}} = \sum_{\alpha_1} U_{i_0; \alpha_1}^{[0]} \lambda^{[0]}_{\alpha_1 ; \alpha_1} V_{\alpha_1;i_1\ldots i_{\ell-1}}^{[0]},$$

where we do not write explicitly the possible values of the dummy index $\alpha _1$. This can be easily inferred by looking at the dimensions of the two index vectors separated by the semicolon on the right-hand side. For this first step in our derivation $\alpha _1$ can take $\min \{c_0, c_1 \times \cdots \times c_{\ell -1}\}$ values. We can now look at the tensor $V_{\alpha _1,i_1\ldots i_{\ell -1}}$ and apply once more the SVDs grouping the first two indices together:

(H5)$$V_{\alpha_1,i_1;i_2\ldots i_{\ell-1}}^{[0]} = \sum_{\alpha_2} U_{\alpha_1,i_1;\alpha_2}^{[1]} \lambda^{[1]}_{\alpha_2;\alpha_2} V_{\alpha_2;i_2\ldots i_{\ell-1}}^{[1]},$$

which we can use to write

(H6)$$R_{i_0 i_1\ldots i_{\ell-1}} = \sum_{\alpha_1,\alpha_2 } U_{i_0; \alpha_1}^{[0]} \lambda^{[0]}_{\alpha_1 ; \alpha_1} U_{\alpha_1,i_1;\alpha_2}^{[1]} \lambda^{[1]}_{\alpha_2;\alpha_2} V_{\alpha_2;i_2\ldots i_{\ell-1}}^{[1]}.$$

Now note that the leftover higher-rank tensor $V^{[1]}$ has $\ell -1$ indices, one less than our original tensor $R$. We can continue as before and write the SVD of $V^{[1]}$ separating $\alpha _2$ and $i_2$ from the rest of the indices

(H7)$$V_{\alpha_2,i_2;i_3\ldots i_{\ell-1}}^{[1]} = \sum_{\alpha_3} U_{\alpha_2,i_2;\alpha_3}^{[2]} \lambda^{[2]}_{\alpha_3;\alpha_3} V_{\alpha_3;i_3\ldots i_{\ell-1}}^{[2]}.$$

We can then plug this decomposition into H6 to obtain

(H8)$$R_{i_0 i_1\ldots i_{\ell-1}} = \sum_{\alpha_1,\alpha_2 } U_{i_0; \alpha_1}^{[0]} \lambda^{[0]}_{\alpha_1 ; \alpha_1} U_{\alpha_1,i_1;\alpha_2}^{[1]} \lambda^{[1]}_{\alpha_2;\alpha_2} V_{\alpha_2;i_2\ldots i_{\ell-1}}^{[1]},$$

where we note that the leftover higher-order tensor $V^{[2]}$ has $\ell -2$ indices. By this point it is easy to guess that in the $k$th step of this recursive decomposition, we will have a tensor $V^{[k]}$ with $\ell -k$ indices $\{\alpha _{k+1},i_{k+1}\ldots,i_{\ell -1}\}$ on which we will apply the SVD to obtain

(H9)$$V_{\alpha_{k+1},i_{k+1};i_{k+2}\ldots i_{\ell-1}}^{[k]} = \sum_{\alpha_{k+2}} U^{[k+1]}_{\alpha_{k+1},i_{k+1}; \alpha_{k+2}} \lambda^{[k+1]}_{\alpha_{k+2};\alpha_{k+2}} V^{[k+1]}_{\alpha_{k+2};i_{k+2}\ldots i_{\ell-1}}.$$

Based on this iterative decomposition we finally arrive at the first matrix product form of the tensor $R$,

(H10)$$R_{i_0 i_1\ldots i_{\ell-1}} = \sum_{\alpha_1,\alpha_2,\ldots, \alpha_{\ell-1}} U_{i_0; \alpha_1}^{[0]} \lambda^{[0]}_{\alpha_1 ; \alpha_1} U_{\alpha_1,i_1;\alpha_2}^{[1]} \lambda^{[1]}_{\alpha_2;\alpha_2} \ldots U_{\alpha_{\ell-2},i_{\ell-2};\alpha_{\ell-1}}^{[1]} \lambda^{[1]}_{\alpha_{\ell-1};\alpha_{\ell-1}} V_{\alpha_{\ell-1},i_{\ell-1}}$$

which is the form used by Yanagimoto et al. [30]. A second useful form can be obtained by introducing $(S^{[k]})^{\alpha _k}_{i_k,\alpha _{k+1}} \equiv U_{\alpha _k,i_k;\alpha _{k+1}}^{[k]} \lambda ^{[k]}_{\alpha _{k+1};\alpha _{k+1}}$ as used in the main text.

Acknowledgments

N. Q. thanks A. M. Brańczyk, M. D Vidrighin, G. Triginer, A. L. Grimsmo, and I. Dhand for useful discussions and acknowledges support from the Ministère de l’Économie et de l’Innovation du Québec. M. L. acknowledges the support from Ministero dell’Istruzione, dell’Università e della Ricerca (Dipartimenti di Eccellenza Program (2018–2022)). J. E. S. acknowledges support from the Natural Sciences and Engineering Research Council of Canada.

Disclosures

The authors declare no conflicts of interest.

Data availability

No data were generated or analyzed in the presented research.

References and Notes

1. S. E. Harris, M. K. Oshman, and R. L. Byer, “Observation of tunable optical parametric fluorescence,” Phys. Rev. Lett. 18, 732–734 (1967). [CrossRef]  

2. D. Magde and H. Mahr, “Study in ammonium dihydrogen phosphate of spontaneous parametric interaction tunable from 4400 to 16 000 Å,” Phys. Rev. Lett. 18, 905–907 (1967). [CrossRef]  

3. A. I. Lvovsky, “Squeezed light,” in Photonics: Scientific Foundations, Technology and Applications (2015), Vol. 1, pp. 121–163. [CrossRef]  

4. U. L. Andersen, T. Gehring, C. Marquardt, and G. Leuchs, “30 years of squeezed light generation,” Phys. Scr. 91, 053001 (2016). [CrossRef]  

5. C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” Rev. Mod. Phys. 84, 621–669 (2012). [CrossRef]  

6. S. L. Braunstein, “Squeezing as an irreducible resource,” Phys. Rev. A 71, 055801 (2005). [CrossRef]  

7. T. Rudolph, “Why I am optimistic about the silicon-photonic route to quantum computing,” APL Photonics 2, 030901 (2017). [CrossRef]  

8. J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, K. K. Sabapathy, N. C. Menicucci, and I. Dhand, “Blueprint for a scalable photonic fault-tolerant quantum computer,” Quantum 5, 392 (2021). [CrossRef]  

9. T. R. Bromley, J. M. Arrazola, S. Jahangiri, J. Izaac, N. Quesada, A. D. Gran, M. Schuld, J. Swinarton, Z. Zabaneh, and N. Killoran, “Applications of near-term photonic quantum computers: software and algorithms,” Quantum Sci. Technol. 5, 034010 (2020). [CrossRef]  

10. M. V. Larsen, C. Chamberland, K. Noh, J. S. Neergaard-Nielsen, and U. L. Andersen, “Fault-tolerant continuous-variable measurement-based quantum computation architecture,” PRX Quantum 2, 030325 (2021). [CrossRef]  

11. M. Kues, C. Reimer, P. Roztocki, L. R. Cortés, S. Sciara, B. Wetzel, Y. Zhang, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, J. Azaña, and R. Morandotti, “On-chip generation of high-dimensional entangled quantum states and their coherent control,” Nature 546, 622–626 (2017). [CrossRef]  

12. C. Reimer, S. Sciara, P. Roztocki, M. Islam, L. Romero Cortés, Y. Zhang, B. Fischer, S. Loranger, R. Kashyap, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, W. J. Munro, J. Azaña, M. Kues, and R. Morandotti, “High-dimensional one-way quantum processing implemented on d-level cluster states,” Nat. Phys. 15, 148–153 (2019). [CrossRef]  

13. J. Aasi, J. Abadie, B. Abbott, et al., “Enhanced sensitivity of the LIGO gravitational wave detector by using squeezed states of light,” Nat. Photonics 7, 613–619 (2013). [CrossRef]  

14. E. Polzik, J. Carri, and H. Kimble, “Spectroscopy with squeezed light,” Phys. Rev. Lett. 68, 3020–3023 (1992). [CrossRef]  

15. C. A. Casacio, L. S. Madsen, A. Terrasson, M. Waleed, K. Barnscheidt, B. Hage, M. A. Taylor, and W. P. Bowen, “Quantum-enhanced nonlinear microscopy,” Nature 594, 201–206 (2021). [CrossRef]  

16. P. J. Mosley, J. S. Lundeen, B. J. Smith, P. Wasylczyk, A. B. U’Ren, C. Silberhorn, and I. A. Walmsley, “Heralded generation of ultrafast single photons in pure quantum states,” Phys. Rev. Lett. 100, 133601 (2008). [CrossRef]  

17. H. Takesue and K. Inoue, “Generation of polarization-entangled photon pairs and violation of Bell’s inequality using spontaneous four-wave mixing in a fiber loop,” Phys. Rev. A 70, 031802 (2004). [CrossRef]  

18. G. Harder, V. Ansari, B. Brecht, T. Dirmeier, C. Marquardt, and C. Silberhorn, “An optimized photon pair source for quantum circuits,” Opt. Express 21, 13975–13985 (2013). [CrossRef]  

19. G. Harder, T. J. Bartley, A. E. Lita, S. W. Nam, T. Gerrits, and C. Silberhorn, “Single-mode parametric-down-conversion states with 50 photons as a source for mesoscopic quantum optics,” Phys. Rev. Lett. 116, 143601 (2016). [CrossRef]  

20. M. A. Finger, T. S. Iskhakov, N. Y. Joly, M. V. Chekhova, and P. S. J. Russell, “Raman-free, noble-gas-filled photonic-crystal fiber source for ultrafast, very bright twin-beam squeezed vacuum,” Phys. Rev. Lett. 115, 143602 (2015). [CrossRef]  

21. J.-Y. Chen, Z. Li, Z. Ma, C. Tang, H. Fan, Y. M. Sua, and Y.-P. Huang, “Photon conversion and interaction on chip,”preprint arXiv:2105.00275 (2021).

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

23. J. B. Surya, X. Guo, C.-L. Zou, and H. X. Tang, “Efficient third-harmonic generation in composite aluminum nitride/silicon nitride microrings,” Optica 5, 103–108 (2018). [CrossRef]  

24. J. Lu, M. Li, C.-L. Zou, A. Al Sayem, and H. X. Tang, “Toward 1% single-photon anharmonicity with periodically poled lithium niobate microring resonators,” Optica 7, 1654–1659 (2020). [CrossRef]  

25. A. Cavanna, J. Hammer, C. Okoth, E. Ortiz-Ricardo, H. Cruz-Ramirez, K. Garay-Palmett, A. B. U’Ren, M. H. Frosz, X. Jiang, N. Y. Joly, and M. V. Chekhova, “Progress toward third-order parametric down-conversion in optical fibers,” Phys. Rev. A 101, 033840 (2020). [CrossRef]  

26. J. Flórez, J. S. Lundeen, and M. V. Chekhova, “Pump depletion in parametric down-conversion with low pump energies,” Opt. Lett. 45, 4264–4267 (2020). [CrossRef]  

27. S. D. Bartlett and B. C. Sanders, “Universal continuous-variable quantum computation: requirement of optical nonlinearity for photon counting,” Phys. Rev. A 65, 042304 (2002). [CrossRef]  

28. S. D. Bartlett and B. C. Sanders, “Efficient classical simulation of optical quantum information circuits,” Phys. Rev. Lett. 89, 207903 (2002). [CrossRef]  

29. U. Chabaud, D. Markham, and F. Grosshans, “Stellar representation of non-Gaussian quantum states,” Phys. Rev. Lett. 124, 063605 (2020). [CrossRef]  

30. R. Yanagimoto, E. Ng, L. G. Wright, T. Onodera, and H. Mabuchi, “Efficient simulation of ultrafast quantum nonlinear optics with matrix product states,” Optica 8, 1306–1315 (2021). [CrossRef]  

31. J. E. Bourassa, N. Quesada, I. Tzitrin, A. Száva, T. Isacsson, J. Izaac, K. K. Sabapathy, G. Dauphinais, and I. Dhand, “Fast simulation of bosonic qubits via Gaussian functions in phase space,” PRX Quantum 2, 040315 (2021). [CrossRef]  

32. M. Walschaers, “Non-Gaussian quantum states and where to find them,” PRX Quantum 2, 030204 (2021). [CrossRef]  

33. Y.-X. Huang, M. Li, K. Lin, Y.-L. Zhang, G.-C. Guo, and C.-L. Zou, “Classical-to-quantum transition in multimode nonlinear systems with strong photon-photon coupling,” preprint arXiv:2111.09557 (2021).

34. R. Yanagimoto, E. Ng, A. Yamamura, T. Onodera, L. G. Wright, M. Jankowski, M. Fejer, P. L. McMahon, and H. Mabuchi, “Onset of non-Gaussian quantum physics in pulsed squeezing with mesoscopic fields,” preprint arXiv:2111.13799 (2021).

35. Y.-X. Huang, M. Li, Z.-J. Chen, X.-B. Zou, G.-C. Guo, and C.-L. Zou, “Quantum effects beyond mean-field treatment in quantum optics,” preprint arXiv:2111.14636 (2021).

36. W. Xing and T. Ralph, “Pump depletion in parametric amplification,”preprint arXiv:2201.01372 (2022).

37. N. Quesada and J. E. Sipe, “Why you should not use the electric field to quantize in nonlinear optics,” Opt. Lett. 42, 3443 (2017). [CrossRef]  

38. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995).

39. M. Born and L. Infeld, “On the quantization of the new field equations I,” Proc. R. Soc. Lond. A 147, 522–546 (1934). [CrossRef]  

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

41. R. Loudon, The Quantum Theory of Light (Oxford University Press, 2000).

42. D. N. Klyshko, Photons and Nonlinear Optics (Routledge, 2018).

43. C. C. Gerry and P. Knight, Introductory Quantum Optics (Cambridge University Press, Cambridge, UK; New York, 2005).

44. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed. (Cambridge University Press, Cambridge; New York, 1999).

45. M. Hillery and L. D. Mlodinow, “Semiclassical expansion for nonlinear dielectric media,” Phys. Rev. A 31, 797–806 (1985). [CrossRef]  

46. M. Hillery and L. Mlodinow, “Quantized fields in a nonlinear dielectric medium: a microscopic approach,” Phys. Rev. A 55, 678–689 (1997). [CrossRef]  

47. P. Drummond and M. Hillery, “Quantum theory of dispersive electromagnetic modes,” Phys. Rev. A 59, 691–707 (1999). [CrossRef]  

48. M. Hillery and L. D. Mlodinow, “Quantization of electrodynamics in nonlinear dielectric media,” Phys. Rev. A 30, 1860–1865 (1984). [CrossRef]  

49. I. Abram and E. Cohen, “Quantum theory for light propagation in a nonlinear effective medium,” Phys. Rev. A 44, 500–517 (1991). [CrossRef]  

50. J. Sipe, N. A. Bhat, P. Chak, and S. Pereira, “Effective field theory for the nonlinear optical properties of photonic crystals,” Phys. Rev. E 69, 016604 (2004). [CrossRef]  

51. M. G. Raymer, “Quantum theory of light in a dispersive structured linear dielectric: a macroscopic Hamiltonian tutorial treatment,” J. Mod. Opt. 67, 196–212 (2020). [CrossRef]  

52. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals (Princeton University Press, 2011).

53. J. Sipe and M. Steel, “A Hamiltonian treatment of stimulated Brillouin scattering in nanoscale integrated waveguides,” New J. Phys. 18, 045004 (2016). [CrossRef]  

54. W. Xie, L. Chang, H. Shu, J. C. Norman, J. D. Peters, X. Wang, and J. E. Bowers, “Ultrahigh-q algaas-on-insulator microresonators for integrated nonlinear photonics,” Opt. Express 28, 32894–32906 (2020). [CrossRef]  

55. J. D. Jackson, Classical Electrodynamics (Wiley, 1998).

56. N. A. Bhat and J. Sipe, “Hamiltonian treatment of the electromagnetic field in dispersive and absorptive structured media,” Phys. Rev. A 73, 063808 (2006). [CrossRef]  

57. A. Judge, M. Steel, J. Sipe, and C. De Sterke, “Canonical quantization of macroscopic electrodynamics in a linear, inhomogeneous magnetoelectric medium,” Phys. Rev. A 87, 033824 (2013). [CrossRef]  

58. J. Sipe, “Photons in dispersive dielectrics,” J. Opt. A: Pure Appl. Opt. 11, 114006 (2009). [CrossRef]  

59. J. Heebner, R. Grover, and T. Ibrahim, Optical Microresonators Theory, Fabrication and Applications (Springer, 2008).

60. Z. Vernon and J. E. Sipe, “Spontaneous four-wave mixing in lossy microring resonators,” Phys. Rev. A 91, 053802 (2015). [CrossRef]  

61. P. Chak, S. Pereira, and J. Sipe, “Coupled-mode theory for periodic side-coupled microcavity and photonic crystal structures,” Phys. Rev. B 73, 035105 (2006). [CrossRef]  

62. M. Liscidini, L. Helt, and J. Sipe, “Asymptotic fields for a Hamiltonian treatment of nonlinear electromagnetic phenomena,” Phys. Rev. A 85, 013833 (2012). [CrossRef]  

63. R. W. Boyd, Nonlinear Optics (Academic Press, 2008).

64. A. Serafini, Quantum Continuous Variables: A Primer of Theoretical Methods (CRC Press, 2017).

65. R. Simon, E. C. G. Sudarshan, and N. Mukunda, “Gaussian pure states in quantum mechanics and the symplectic group,” Phys. Rev. A 37, 3028–3038 (1988). [CrossRef]  

66. B. Dutta, N. Mukunda, and R. Simon, “The real symplectic groups in quantum mechanics and optics,” Pramana 45, 471–497 (1995). [CrossRef]  

67. G. Adesso, S. Ragy, and A. R. Lee, “Continuous variable quantum information: Gaussian states and beyond,” Open Syst. Inf. Dyn. 21, 1440001 (2014). [CrossRef]  

68. R. Shankar, Principles of Quantum Mechanics (Springer Science Business Media, 2012).

69. J. J. Sakurai, Modern Quantum Mechanics (Addison-Wesley, 2011).

70. S. Barnett and P. M. Radmore, Methods in Theoretical Quantum Optics, Vol. 15 (Oxford University Press, 2002).

71. A. B. Klimov and S. M. Chumakov, A Group-theoretical Approach to Quantum Optics: Models of Atom–Field Interactions (John Wiley & Sons, 2009).

72. F. J. Dyson, “The radiation theories of Tomonaga, Schwinger, and Feynman,” Phys. Rev. 75, 486–502 (1949). [CrossRef]  

73. S. Blanes, F. Casas, J. Oteo, and J. Ros, “The magnus expansion and some of its applications,” Phys. Rep. 470, 151–238 (2009). [CrossRef]  

74. W. Magnus, “On the exponential solution of differential equations for a linear operator,” Comm. Pure Appl. Math. 7, 649–673 (1954). [CrossRef]  

75. M. Suzuki, “Generalized Trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems,” Commun. Math. Phys. 51, 183–190 (1976). [CrossRef]  

76. H. F. Trotter, “On the product of semi-groups of operators,” Proc. Am. Math. Soc. 10, 545–551 (1959). [CrossRef]  

77. In the case where Δ_ij is not time-dependent, this can always be assumed without loss of generality as shown in Appendix C.

78. R. A. Horn and C. R. Johnson, Matrix Analysis (Cambridge University Press, 2012).

79. N. Killoran, J. Izaac, N. Quesada, V. Bergholm, M. Amy, and C. Weedbrook, “Strawberry Fields: asoftware platform for photonic quantum computing,” Quantum 3, 129129 (2019). [CrossRef]  

80. B. Gupt, J. Izaac, and N. Quesada, “The Walrus: a library for the calculation of Hafnians, Hermite polynomials and Gaussian boson sampling,” JOSS 4, 1705 (2019). [CrossRef]  

81. A. Lvovsky, W. Wasilewski, and K. Banaszek, “Decomposing a pulsed optical parametric amplifier into independent squeezers,” J. Mod. Opt. 54, 721–733 (2007). [CrossRef]  

82. W. Wasilewski, A. I. Lvovsky, K. Banaszek, and C. Radzewicz, “Pulsed squeezed light: simultaneous squeezing of multiple modes,” Phys. Rev. A 73, 063819 (2006). [CrossRef]  

83. T. Lipfert, D. B. Horoshko, G. Patera, and M. I. Kolobov, “Bloch-Messiah decomposition and magnus expansion for parametric down-conversion with monochromatic pump,” Phys. Rev. A 98, 013815 (2018). [CrossRef]  

84. A. Gatti, R. Zambrini, M. San Miguel, and L. A. Lugiato, “Multiphoton multimode polarization entanglement in parametric down-conversion,” Phys. Rev. A 68, 053807 (2003). [CrossRef]  

85. A. K. Ekert and P. L. Knight, “Relationship between semiclassical and quantum-mechanical input–output theories of optical response,” Phys. Rev. A 43, 3934–3938 (1991). [CrossRef]  

86. A. Christ, B. Brecht, W. Mauerer, and C. Silberhorn, “Theory of quantum frequency conversion and type-II parametric down-conversion in the high-gain regime,” New J. Phys. 15, 053038 (2013). [CrossRef]  

87. K. E. Cahill and R. J. Glauber, “Ordered expansions in boson amplitude operators,” Phys. Rev. 177, 1857–1881 (1969). [CrossRef]  

88. Note that this equation does not imply that $\mathcal {U} = \exp \left \{ \tfrac 12 \sum _{k, l=1}^\ell J_{kl} a^\dagger _{k} a^\dagger _l - \text {H.c.}\right \}$; in particular, we can postmultiply the exponential by any unitary operator that satisfies $\mathcal {U}_{0} | {\text {vac}}\rangle = | {\text {vac}}\rangle$ and still satisfy Eq. (237). The form of the operator $\mathcal {U}_{0}$ is described in Appendix E.

89. The diagonal and unitary matrices obtained from the Takagi decomposition can be related to the diagonal and orthogonal matrices obtained by performing a Bloch–Messiah decomposition of the symplectic matrix that diagonalizes the covariance matrix of the pure state in Eq. (237). This symplectic matrix can be obtained by doing a Williamson decomposition of the covariance matrix of the quantum state as described by, e.g., Serafini [64].

90. A. M. Brańczyk, T. M. Stace, and T. Ralph, “Time ordering in spontaneous parametric down-conversion,” in AIP Conference Proceedings, Vol. 1363 (2011), pp. 335–338.

91. W. P. Grice and I. A. Walmsley, “Spectral information and distinguishability in type-II down-conversion with a broadband pump,” Phys. Rev. A 56, 1627–1634 (1997). [CrossRef]  

92. G. Triginer, M. D. Vidrighin, N. Quesada, A. Eckstein, M. Moore, W. S. Kolthammer, J. E. Sipe, and I. A. Walmsley, “Understanding high-gain twin-beam sources using cascaded stimulated emission,” Phys. Rev. X 10, 031063 (2020). [CrossRef]  

93. V. Bargmann, “On a Hilbert space of analytic functions and an associated integral transform part I,” Comm. Pure Appl. Math. 14, 187–214 (1961). [CrossRef]  

94. I. E. Segal and G. W. Mackey, Mathematical Problems of Relativistic Physics, Vol. 2 (American Mathematical Soc., 1963).

95. R. García-Patrón, J. J. Renema, and V. Shchesnovich, “Simulating boson sampling in lossy architectures,” Quantum 3, 169 (2019). [CrossRef]  

96. D. J. Brod and M. Oszmaniec, “Classical simulation of linear optics subject to nonuniform losses,” Quantum 4, 267 (2020). [CrossRef]  

97. O. F. Thomas, W. McCutcheon, and D. P. McCutcheon, “A general framework for multimode Gaussian quantum optics and photo-detection: application to Hong–Ou–Mandel interference with filtered heralded single photon sources,” APL Photonics 6, 040801 (2021). [CrossRef]  

98. K. Laiho, T. Dirmeier, M. Schmidt, S. Reitzenstein, and C. Marquardt, “Measuring higher-order photon correlations of faint quantum light: a short review,” Phys. Lett. A 435, 128059 (2022). [CrossRef]  

99. A. Christ, K. Laiho, A. Eckstein, K. N. Cassemiro, and C. Silberhorn, “Probing multimode squeezing with correlation functions,” New J. Phys. 13, 033027 (2011). [CrossRef]  

100. C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, “Gaussian boson sampling,” Phys. Rev. Lett. 119, 170501 (2017). [CrossRef]  

101. R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, “Detailed study of Gaussian boson sampling,” Phys. Rev. A 100, 032326 (2019). [CrossRef]  

102. N. Quesada, J. M. Arrazola, and N. Killoran, “Gaussian boson sampling using threshold detectors,” Phys. Rev. A 98, 062322 (2018). [CrossRef]  

103. A. Deshpande, A. Mehta, T. Vincent, N. Quesada, M. Hinsche, M. Ioannou, L. Madsen, J. Lavoie, H. Qi, J. Eisert, D. Hangleiter, B. Fefferman, and I. Dhand, “Quantum computational advantage via high-dimensional Gaussian boson sampling,” Sci. Adv. 8, eabi7894 (2022). [CrossRef]  

104. N. Quesada and J. M. Arrazola, “Exact simulation of Gaussian boson sampling in polynomial space and exponential time,” Phys. Rev. Res. 2, 023005 (2020). [CrossRef]  

105. H.-S. Zhong, L.-C. Peng, Y. Li, Y. Hu, W. Li, J. Qin, D. Wu, W. Zhang, H. Li, L. Zhang, Z. Wang, L. You, X. Jiang, L. Li, N.-L. Liu, J. P. Dowling, C.-Y. Lu, and J.-W. Pan, “Experimental Gaussian boson sampling,” Sci. Bull. 64, 511–515 (2019). [CrossRef]  

106. H.-S. Zhong, H. Wang, Y.-H. Deng, et al., “Quantum computational advantage using photons,” Science 370, 1460–1463 (2020). [CrossRef]  

107. L. G. Helt and N. Quesada, “Degenerate squeezing in waveguides: a unified theoretical approach,” JPhys Photonics 2, 035001 (2020). [CrossRef]  

108. N. Quesada, G. Triginer, M. D. Vidrighin, and J. E. Sipe, “Theory of high-gain twin-beam generation in waveguides: from Maxwell’s equations to efficient simulation,” Phys. Rev. A 102, 033519 (2020). [CrossRef]  

109. G. P. Agrawal, Nonlinear Fiber Optics, 4th ed. (Academic Press, Burlington, MA, 2007).

110. M.-D. Vidrighin, “Quantum optical measurements for practical estimation and information thermodynamics,” Ph.D. thesis, Imperial College London (2017).

111. C. McKinstrie and M. Marhic, “Quadrature and number fluctuations produced by parametric devices driven by pulsed pumps,” Opt. Express 21, 19437–19466 (2013). [CrossRef]  

112. D. B. Horoshko, L. La Volpe, F. Arzani, N. Treps, C. Fabre, and M. I. Kolobov, “Bloch–Messiah reduction for twin beams of light,” Phys. Rev. A 100, 013837 (2019). [CrossRef]  

113. M. I. Kolobov, “The spatial behavior of nonclassical light,” Rev. Mod. Phys. 71, 1539–1589 (1999). [CrossRef]  

114. Z. Yang, M. Liscidini, and J. E. Sipe, “Spontaneous parametric down-conversion in waveguides: a backward heisenberg picture approach,” Phys. Rev. A 77, 033808 (2008). [CrossRef]  

115. N. Quesada and J. E. Sipe, “Effects of time ordering in quantum nonlinear optics,” Phys. Rev. A 90, 063840 (2014). [CrossRef]  

116. A. M. Brańczyk, A. Fedrizzi, T. M. Stace, T. C. Ralph, and A. G. White, “Engineered optical nonlinearity for quantum light sources,” Opt. Express 19, 55–65 (2011). [CrossRef]  

117. P. B. Dixon, J. H. Shapiro, and F. N. C. Wong, “Spectral engineering by Gaussian phase-matching for quantum photonics,” Opt. Express 21, 5879–5890 (2013). [CrossRef]  

118. J.-L. Tambasco, A. Boes, L. G. Helt, M. J. Steel, and A. Mitchell, “Domain engineering algorithm for practical and effective photon sources,” Opt. Express 24, 19616–19626 (2016). [CrossRef]  

119. A. Dosseva, L. Cincio, and A. M. Brańczyk, “Shaping the joint spectrum of down-converted photons through optimized custom poling,” Phys. Rev. A 93, 013801 (2016). [CrossRef]  

120. F. Graffitti, D. Kundys, D. T. Reid, A. M. Brańczyk, and A. Fedrizzi, “Pure down-conversion photons through sub-coherence-length domain engineering,” Quantum Sci. Technol. 2, 035001 (2017). [CrossRef]  

121. A. B. U’Ren, C. Silberhorn, R. Erdmann, K. Banaszek, W. P. Grice, I. A. Walmsley, and M. G. Raymer, “Generation of pure-state single-photon wavepackets by conditional preparation based on spontaneous parametric downconversion,” arXiv preprint quant-ph/0611019 (2006).

122. N. Quesada and A. M. Brańczyk, “Gaussian functions are optimal for waveguided nonlinear-quantum-optical processes,” Phys. Rev. A 98, 043813 (2018). [CrossRef]  

123. N. Quesada and J. E. Sipe, “Time-ordering effects in the generation of entangled photons using nonlinear optical processes,” Phys. Rev. Lett. 114, 093903 (2015). [CrossRef]  

124. N. Quesada, “Very nonlinear quantum optics,” Ph.D. thesis, University of Toronto (2015).

125. D. G. Rabus, Integrated Ring Resonators (Springer, 2007).

126. Z. Vernon, “Microresonators for nonlinear quantum optics,” Ph.D. thesis, University of Toronto (2017).

127. A. Yariv and P. Yeh, Photonics: Optical Electronics in Modern Communications (Oxford University Press, 2007).

128. Z. Vernon, M. Menotti, C. Tison, J. Steidle, M. Fanto, P. Thomas, S. Preble, A. Smith, P. Alsing, M. Liscidini, and J. E. Sipe, “Truly unentangled photon pairs without spectral filtering,” Opt. Lett. 42, 3638–3641 (2017). [CrossRef]  

129. J. B. Christensen, J. G. Koefoed, K. Rottwitt, and C. McKinstrie, “Engineering spectrally unentangled photon pairs from nonlinear microring resonators by pump manipulation,” Opt. Lett. 43, 859–862 (2018). [CrossRef]  

130. V. D. Vaidya, B. Morrison, L. Helt, R. Shahrokshahi, D. Mahler, M. Collins, K. Tan, J. Lavoie, A. Repingon, M. Menotti, N. Quesada, R. C. Pooser, A. E. Lita, T. Gerrits, S. W. Nam, and Z. Vernon, “Broadband quadrature-squeezed vacuum and nonclassical photon number correlations from a nanophotonic device,” Sci. Adv. 6, eaba9186 (2020). [CrossRef]  

131. Z. Vernon, N. Quesada, M. Liscidini, B. Morrison, M. Menotti, K. Tan, and J. Sipe, “Scalable squeezed-light source for continuous-variable quantum sampling,” Phys. Rev. Appl. 12, 064024 (2019). [CrossRef]  

132. Y. Zhao, Y. Okawachi, J. K. Jang, X. Ji, M. Lipson, and A. L. Gaeta, “Near-degenerate quadrature-squeezed vacuum generation on a silicon-nitride chip,” Phys. Rev. Lett. 124, 193601 (2020). [CrossRef]  

133. Y. Zhang, M. Menotti, K. Tan, V. Vaidya, D. Mahler, L. Helt, L. Zatti, M. Liscidini, B. Morrison, and Z. Vernon, “Squeezed light from a nanophotonic molecule,” Nat. Commun. 12, 1–6 (2021). [CrossRef]  

134. H. Seifoory, Z. Vernon, D. Mahler, M. Menotti, Y. Zhang, and J. E. Sipe, “Degenerate squeezing in a dual-pumped integrated microresonator: parasitic processes and their suppression,”preprint arXiv:2109.03298 (2021).

135. M. Chen, N. C. Menicucci, and O. Pfister, “Experimental realization of multipartite entanglement of 60 modes of a quantum optical frequency comb,” Phys. Rev. Lett. 112, 120505 (2014). [CrossRef]  

136. Y. K. Chembo, “Quantum dynamics of Kerr optical frequency combs below and above threshold: spontaneous four-wave mixing, entanglement, and squeezed states of light,” Phys. Rev. A 93, 033820 (2016). [CrossRef]  

137. C. Reimer, M. Kues, P. Roztocki, B. Wetzel, F. Grazioso, B. E. Little, S. T. Chu, T. Johnston, Y. Bromberg, L. Caspani, D. J. Moss, and R. Morandotti, “Generation of multiphoton entangled quantum states by means of integrated frequency combs,” Science 351, 1176–1180 (2016). [CrossRef]  

138. S. Hagness, D. Rafizadeh, S. Ho, and A. Taflove, “FDTD microcavity simulations: design and experimental realization of waveguide-coupled single-mode ring and whispering-gallery-mode disk resonators,” J. Lightwave Technol. 15, 2154–2165 (1997). [CrossRef]  

139. H. P. Yuen and V. W. Chan, “Noise in homodyne and heterodyne detection,” Opt. Lett. 8, 177–179 (1983). [CrossRef]  

140. R. Wolf, Y. Jia, S. Bonaus, C. S. Werner, S. J. Herr, I. Breunig, K. Buse, and H. Zappe, “Quasi-phase-matched nonlinear optical frequency conversion in on-chip whispering galleries,” Optica 5, 872–875 (2018). [CrossRef]  

141. J.-Y. Chen, Z.-H. Ma, Y. M. Sua, Z. Li, C. Tang, and Y.-P. Huang, “Ultra-efficient frequency conversion in quasi-phase-matched lithium niobate microrings,” Optica 6, 1244–1245 (2019). [CrossRef]  

142. J. Lu, J. B. Surya, X. Liu, A. W. Bruch, Z. Gong, Y. Xu, and H. X. Tang, “Periodically poled thin-film lithium niobate microring resonators with a second-harmonic generation efficiency of 250,000%/w,” Optica 6, 1455–1460 (2019). [CrossRef]  

143. R. Luo, Y. He, H. Liang, M. Li, J. Ling, and Q. Lin, “Optical parametric generation in a lithium niobate microring with modal phase matching,” Phys. Rev. Appl. 11, 034026 (2019). [CrossRef]  

144. R. Luo, H. Jiang, S. Rogers, H. Liang, Y. He, and Q. Lin, “On-chip second-harmonic generation and broadband parametric down-conversion in a lithium niobate microresonator,” Opt. Express 25, 24531–24539 (2017). [CrossRef]  

145. Z. Yang, P. Chak, A. D. Bristow, H. M. van Driel, R. Iyer, J. S. Aitchison, A. L. Smirl, and J. Sipe, “Enhanced second-harmonic generation in AlGaAs microring resonators,” Opt. Lett. 32, 826–828 (2007). [CrossRef]  

146. J. Lin, Y. Xu, J. Ni, M. Wang, Z. Fang, L. Qiao, W. Fang, and Y. Cheng, “Phase-matched second-harmonic generation in an on-chip LiNbO₃ microresonator,” Phys. Rev. Appl. 6, 014002 (2016). [CrossRef]  

147. G. Lin, J. U. Fürst, D. V. Strekalov, and N. Yu, “Wide-range cyclic phase matching and second harmonic generation in whispering gallery resonators,” Appl. Phys. Lett. 103, 181107 (2013). [CrossRef]  

148. L. G. Helt, M. Liscidini, and J. E. Sipe, “How does it scale? Comparing quantum and classical nonlinear optical processes in integrated devices,” J. Opt. Soc. Am. B 29, 2199–2212 (2012). [CrossRef]  

149. M. Liscidini and J. Sipe, “Stimulated emission tomography,” Phys. Rev. Lett. 111, 193602 (2013). [CrossRef]  

150. A. De Rossi, V. Berger, M. Calligaro, G. Leo, V. Ortiz, and X. Marcadet, “Parametric fluorescence in oxidized aluminum gallium arsenide waveguides,” Appl. Phys. Lett. 79, 3758–3760 (2001). [CrossRef]  

151. V. Introini, M. J. Steel, J. E. Sipe, L. G. Helt, and M. Liscidini, “Spontaneous parametric down conversion in a doubly resonant one-dimensional photonic crystal,” Opt. Lett. 45, 1244 (2020). [CrossRef]  

152. S. Azzini, D. Grassani, M. Galli, L. C. Andreani, M. Sorel, M. J. Strain, L. G. Helt, J. E. Sipe, M. Liscidini, and D. Bajoni, “From classical four-wave mixing to parametric fluorescence in silicon microring resonators,” Opt. Lett. 37, 3807 (2012). [CrossRef]  

153. L. A. Rozema, C. Wang, D. H. Mahler, A. Hayat, A. M. Steinberg, J. E. Sipe, and M. Liscidini, “Characterizing an entangled-photon source with classical detectors and measurements,” Optica 2, 430 (2015). [CrossRef]  

154. B. Fang, M. Liscidini, J. E. Sipe, and V. O. Lorenz, “Multidimensional characterization of an entangled photon-pair source via stimulated emission tomography,” Opt. Express 24, 10013 (2016). [CrossRef]  

155. M. A. Ciampini, A. Geraldi, V. Cimini, C. Macchiavello, J. E. Sipe, M. Liscidini, and P. Mataloni, “Stimulated emission tomography: beyond polarization,” Opt. Lett. 44, 41 (2019). [CrossRef]  

156. A. Eckstein, G. Boucher, A. Lemaître, P. Filloux, I. Favero, G. Leo, J. E. Sipe, M. Liscidini, and S. Ducci, “High-resolution spectral characterization of two photon states via classical measurements,” Laser Photonics Rev. 8, L76–L80 (2014). [CrossRef]  

157. B. Fang, O. Cohen, M. Liscidini, J. E. Sipe, and V. O. Lorenz, “Fast and highly resolved capture of the joint spectral density of photon pairs,” Optica 1, 281 (2014). [CrossRef]  

158. I. Jizan, L. G. Helt, C. Xiong, M. J. Collins, D.-Y. Choi, C. Joon Chae, M. Liscidini, M. J. Steel, B. J. Eggleton, and A. S. Clark, “Bi-photon spectral correlation measurements from a silicon nanowire in the quantum and classical regimes,” Sci. Rep. 5, 12557 (2015). [CrossRef]  

159. D. Grassani, A. Simbula, S. Pirotta, M. Galli, M. Menotti, N. C. Harris, T. Baehr-Jones, M. Hochberg, C. Galland, M. Liscidini, and D. Bajoni, “Energy correlations of photon pairs generated by a silicon microring resonator probed by Stimulated Four Wave Mixing,” Sci. Rep. 6, 23564 (2016). [CrossRef]  

160. I. Jizan, B. Bell, L. G. Helt, A. C. Bedoya, C. Xiong, and B. J. Eggleton, “Phase-sensitive tomography of the joint spectral amplitude of photon pair sources,” Opt. Lett. 41, 4803–4806 (2016). [CrossRef]  

161. M. Borghi, “Phase-resolved joint spectra tomography of a ring resonator photon pair source using a silicon photonic chip,” Opt. Express 28, 7442 (2020). [CrossRef]  

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

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

164. S. Montangero, Introduction to Tensor Network Methods (Springer, 2018).

165. C. Schön, E. Solano, F. Verstraete, J. I. Cirac, and M. M. Wolf, “Sequential generation of entangled multiqubit states,” Phys. Rev. Lett. 95, 110503 (2005). [CrossRef]  

166. C. Schön, K. Hammerer, M. M. Wolf, J. I. Cirac, and E. Solano, “Sequential generation of matrix-product states in cavity QED,” Phys. Rev. A 75, 032311 (2007). [CrossRef]  

167. L. G. Helt, Z. Yang, M. Liscidini, and J. Sipe, “Spontaneous four-wave mixing in microring resonators,” Opt. Lett. 35, 3006–3008 (2010). [CrossRef]  

168. I. Dhand, M. Engelkemeier, L. Sansoni, S. Barkhofen, C. Silberhorn, and M. B. Plenio, “Proposal for quantum simulation via all-optically-generated tensor network states,” Phys. Rev. Lett. 120, 130501 (2018). [CrossRef]  

169. C. Roberts, A. Milsted, M. Ganahl, A. Zalcman, B. Fontaine, Y. Zou, J. Hidary, G. Vidal, and S. Leichenauer, “Tensornetwork: a library for physics and machine learning,”preprint arXiv:1905.01330 (2019).

170. M. Fishman, S. R. White, and E. M. Stoudenmire, “The ITensor software library for tensor network calculations,”preprint arXiv:2007.14822 (2020).

171. 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]  

172. 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.: Theory Exp. 2011, P05001 (2011). [CrossRef]  

Nicolás Quesada is an Assistant Professor in the Engineering Physics department at Polytechnique Montréal and holds an MEI research chair in quantum photonics. His group develops the theory and computational tools underpinning the next generation of nonclassical bright light sources and devices needed for building fault-tolerant quantum computers, quantum communication networks, and quantum sensors. Prior to Polytechnique, he worked at Xanadu Quantum Technologies as a lead developer of Strawberry Fields and The Walrus software libraries and led their theoretical efforts to demonstrate photonic quantum advantage. Before joining industry, he was an NSERC postdoctoral fellow and, as a PhD student, a Vanier and Stoicheff scholar at the University of Toronto.

Lukas G. Helt is a quantum photonics engineer and software developer at Xanadu Quantum Technologies. He received his PhD degree in physics in 2013 from the University of Toronto (Canada). Prior to joining Xanadu he was project leader of the CUDOS Nonlinear Quantum Photonics project while holding a research fellowship at Macquarie University (Australia).

Matteo Menotti received his PhD degree in Physics from the University of Pavia (Italy) in 2018 with a dissertation on the generation of nonclassical states of light in integrated devices. He has been a quantum photonics engineer at Xanadu Quantum Technologies since 2018, where he works on the design and simulation of integrated sources of squeezed states for photonic quantum computation.

Marco Liscidini is Associate Professor at the Department of Physics of the University of Pavia (Italy), and he serves as technical advisor to Xanadu Quantum Technologies. He received the PhD degree in physics from the University of Pavia (Italy) in 2006. From 2007 to 2009, he was a Post-Doctoral Fellow at the Department of Physics of the University of Toronto (Canada). His research activity focuses on the theoretical study and modeling of the light–matter interaction in micro- and nanostructures. He works in several areas of photonics, including classical and quantum nonlinear optics, photonic crystal structures, and optical sensing. His theoretical research activity is in strong collaboration with experimental groups and in the framework of national, European, US, and Canadian research programs. He is a Fellow of Optica.

John E. Sipe is a Professor of Physics at the University of Toronto, where he has been on the academic staff since 1981; he is a theorist who interacts strongly with experimental groups. His research interests are in quantum optics and condensed matter physics, particularly on topics such as quantum interference processes in solids, quantum nonlinear optics in integrated photonic structures, and foundational issues in the description of the interaction of light with matter. He is a Fellow of Optica, of the American Physical Society, and of the Royal Society of Canada.