# Simple way to incorporate loss when modelling multimode entangled state
generation

Colin Vendromin<EMAIL_ADDRESS>Marc M. Dignam  Department of
Physics, Engineering Physics, and Astronomy, Queen’s University, Kingston,
Ontario K7L 3N6, Canada

###### Abstract

We show that the light generated via spontaneous four-wave mixing or
parametric down conversion in multiple, coupled, lossy cavities is a multimode
squeezed thermal state. Requiring this state to be the solution of the
Lindblad master equation results in a set of coupled first-order differential
equations for the time-dependent squeezing parameters and thermal photon
numbers of the state. The benefit of this semi-analytic approach is that the
number of coupled equations scales linearly with the number of modes and is
independent of the number of photons generated. With this analytic form of the
state, correlation variances are easily expressed as analytic functions of the
time-dependent mode parameters. Thus, our solution makes it computationally
tractable and relatively straight forward to calculate the generation and
evolution of multimode entangled states in multiple coupled, lossy cavities,
even when there are a large number of modes and/or photons.

††preprint: APS/123-QED

##  I Introduction

Multimode squeezed states can be generated via a nonlinear interaction in
resonant structures, such as ring resonators (see Fig. 1(a)) or coupled-
resonator optical waveguides (CROWs) in a photonic crystal (see Fig. 1(c)).
They are a source of continuous-variable (CV) entanglement, since the
quadratures of the photons in different modes in the state can be correlated.
CV entanglement has applications in boson sampling [1, 2], quantum computing
[3, 4, 5], and CV cluster states [6, 7, 8].

The theoretical generation of multimode squeezed states via spontaneous four-
wave mixing (SFWM) or spontaneous parametric down conversion (SPDC) has been
studied extensively for ring resonators [9], nonlinear waveguides [10, 11],
and CROWs [12, 13]. Photon loss is an important problem in these systems,
since it can reduce the squeezing and inseparability of the state.

Loss in ring resonators and waveguides due to photon scattering can be handled
by introducing reservoir modes that photons can couple in to [9, 14]. The
waveguide-reservoir coupling parameters can be estimated with phenomenological
values taken from experiment. In contrast, loss in CROWs or coupled-cavities
(see Fig. 1(d)) can be handled intrinsically by calculating the complex
frequencies of the CROW Bloch modes or cavity modes [12, 13]. In this
approach, the evolution in these lossy systems is expressed as the non-unitary
evolution of the reduced density operator of the generated light, obtained
from the solution of the Lindblad master equation. This is the approach that
we use in this work.

(a)

(b)

(c)

(d)

Figure 1: Schematics of structures of potential interest. The thick arrow
indicates where the pump field is injected. The first structure (a) is a ring
resonator coupled to a waveguide. The second (b) are a series of coupled ring
resonators coupled to an input and output waveguide. The thin black arrows in
(a) and (b) indicate the direction of the field components at the coupling
points. The third (c) is a CROW made by periodically removing air holes (shown
as white circles) along a row in the crystal. And the fourth (d) is three
high-Q cavities that are coupled to a waveguide.

Before discussing our approach, we present in Fig. 1, four different lossy
coupled-cavity structures in which our theoretical method can be used to
determine the quantum state of light generated via SFWM or SPDC. In Fig. 1(a),
a ring resonator coupled to a waveguide is shown. The pump pulse is injected
into the straight waveguide; it then couples into the ring, where it generates
multimode entangled light via SFWM or SPDC in the multiple, lossy modes of the
ring. In Fig. 1(b) a series of three coupled rings is shown. This structure
allows for the generation and propagation of the multimode state in the
coupled-rings. Fig. 1(c) shows a CROW in a photonic crystal. The pump pulse
propagates in a set of Bloch modes of the CROW and generates signal and idler
photons as it propagates. Lastly, Fig. 1(d) shows a three-mode cavity coupled
to a waveguide in a photonic crystal. The state is predominantly generated
inside the cavity and couples out into the waveguide.

In this work, we prove that the density operator for entangled light generated
in multiple lossy modes via SFWM or SPDC takes the form of a multimode
squeezed thermal state (MSTS), where the thermal part of the state captures
the effects of photon loss. The loss can be handled either phenomenologically
or intrinsically, depending on the structure. For $M$ modes, we derive a set
of $3M$ coupled first-order differential equations that provide the complete
evolution of the squeezing parameters and thermal photon numbers in the
density operator. The number of equations is independent of the number of
generated photons which is generally many fewer equations than must be solved
when using numerical methods that rely on Fock states [15], since in those
methods the number of equations can be as large as $(M+N)![M!N!]^{-1}$, where
$N$ is the total photon number.

Our solution of the Lindblad master equation applies to the discrete modes of
a structure; for example, the modes of the rings in Figs. 1(a) and 1(b), and
the Bloch modes of the CROW in Fig. 1(c). Our solution can also be applied to
the structure in Fig. 1(d) as long as the modes are discretized (e.g. by using
periodic boundary conditions).

In the limit that the modes are lossless, our formalism reproduces the results
for a multimode squeezed vacuum state given by N. Quesada et al. [16]. Also,
in the case of only one or two lossy modes, our multimode formalism reproduces
the results of previous work [17, 18, 19].

The rest of the paper is organized as follows. In Sec. II we define the lossy
modes of a structure, called quasimodes. In Sec. III we present the nonlinear
Hamiltonian for SFWM, and the Lindblad master equation in the basis of
quasimodes. In Sec. IV we show that the analytic solution to the Lindblad
master equation is the density operator for a multimode squeezed thermal
state, and obtain a set of coupled first-order differential equations for the
squeezing amplitudes, squeezing phases, and thermal photons numbers. In Sec. V
we give expressions for the correlation variance and the expectation value of
the photon number operator. In Sec. VI we show that in the limit when the
modes are lossless our coupled equations describe a squeezed vacuum state that
agrees with other work, and that they produce results that are consistent with
previous work on single- and two- mode squeezed states. In Sec. VII we present
our results from solving the coupled differential equations for a squeezing
process in the Bloch modes of a four-cavity CROW. Finally, in Sec. VIII we
conclude. Additionally, there are three appendices. In Appendix A we derive
the Takagi factorization, starting from the symmetric singular value
decomposition. In Appendix B we give the details on deriving the coupled
equations for the squeezing process, and in Appendix C we discuss how to
numerically solve them for a general system initially in the vacuum state.

##  II Quasimodes

The positive frequency part of the electric field in a lossy mode $m$ takes
the form

$\tilde{\bm{E}}^{(+)}_{m}(\bm{r},t)=\tilde{\bm{N}}_{m}(\bm{r}){\rm
e}^{-i\tilde{\omega}_{m}t},$ (1)

where $\tilde{\bm{N}}_{m}(\bm{r})$ is the spatial profile of the mode and we
define $\tilde{\omega}_{m}=\omega_{m}-i\gamma_{m}$ the complex frequency with
real part $\omega_{m}$, and the imaginary part $\gamma_{m}$ quantifies the
energy leakage. The spatial profile solves the Helmholtz equation

$\displaystyle\bm{\nabla}\times\bm{\nabla}\times\tilde{\bm{N}}_{m}(\bm{r})-\frac{\tilde{\omega}^{2}_{m}}{c^{2}}\epsilon(\bm{r})\tilde{\bm{N}}_{m}(\bm{r})$
$\displaystyle=0,$ (2)

where $\epsilon(\bm{r})$ is the relative dielectric function of the structure.
Throughout this paper we refer to the solutions of Eq. (2) as quasimodes [20],
and we assume the quasimodes are orthonormal in the following sense

$\displaystyle\int
d^{3}\bm{r}\epsilon(\bm{r})\tilde{\bm{N}}^{*}_{m}(\bm{r})\cdot\tilde{\bm{N}}_{l}(\bm{r})$
$\displaystyle=\delta_{ml}.$ (3)

The orthogonality of the modes is strictly obeyed if there is symmetry in the
structure (e.g. rotational or translational). For example, the Bloch modes of
the CROW in Fig. 1(c) are orthogonal due to the translational symmetry of the
CROW, and the modes of the ring resonator are orthogonal due to the rotational
symmetry of the ring. In structures that lack symmetry, such as the cavity
coupled to a waveguide in Fig. 1(d), the orthonormality in Eq. (3) is not
generally strictly obeyed. However, modes that are separated in frequency by
more than their linewidths will be very nearly orthogonal [21]. Therefore, in
what follows, we will assume mode orthogonality.

In the next section we present a quantized theory of spontaneous four-wave
mixing and the Lindblad master equation in the basis of these quasimodes.

##  III The Lindblad master equation for a multimode lossy structure

Before discussing the solution to the Lindblad master equation for multiple
lossy modes, we present the system Hamiltonian in the undepleted pump
approximation and the Lindblad master equation that we reference throughout
the paper.

In a lossy dielectric structure, a discrete basis of orthogonal quasimodes can
be used to construct the system Hamiltonian

$\displaystyle H=H_{\rm L}+H_{\rm NL},$ (4)

where linear Hamiltonian in the basis of the quasimodes is

$H_{\rm L}=\sum_{m}\hbar\tilde{\omega}_{m}b^{\dagger}_{m}b_{m},$ (5)

where $b^{\dagger}_{m}$ and $b_{m}$ are the creation and annihilation
operators for photons in the $m{\rm th}$ quasimode satisfying the standard
commutation relation

$\left[b_{m},b^{\dagger}_{l}\right]=\delta_{ml}.$ (6)

Now we include the SFWM nonlinear interaction where two pump photons $p_{1}$
and $p_{2}$ are annihilated to create a pair of signal and idler photons $m$
and $l$. The nonlinear Hamiltonian in the basis of the quasimodes is [22]

$H_{\rm
NL}=\sum_{m,l,p_{1},p_{2}}G_{mlp_{1}p_{2}}b^{\dagger}_{m}b^{\dagger}_{l}b_{p_{1}}b_{p_{2}}+{\rm
H.c.},$ (7)

where

$\displaystyle G_{mlp_{1}p_{2}}$
$\displaystyle\equiv\frac{9\hbar^{2}}{16\epsilon_{0}}\sqrt{\omega_{m}\omega_{l}\omega_{p_{1}}\omega_{p_{2}}}\sum_{i,j,k,h}\int
d^{3}\bm{r}\chi^{(3)}_{ijkh}(\bm{r})$
$\displaystyle\times\tilde{N}^{*}_{mi}(\bm{r})\tilde{N}^{*}_{lj}(\bm{r})\tilde{N}_{p_{1}k}(\bm{r})\tilde{N}_{p_{2}h}(\bm{r})$
(8)

is the effective nonlinear parameter that depends on the spatial overlap of
the quasimodes $\tilde{\bm{N}}_{m}(\bm{r})$. The subscripts $i,j,k$, and $h$
label the Cartesian components of the electric field and the medium nonlinear
tensor $\chi^{(3)}$. In deriving this result we have assumed that the
imaginary part of the frequencies are small, such that we can neglect them
under the square root, so that
$\sqrt{\tilde{\omega}_{m}}\approx\sqrt{\omega_{m}}$. This is valid for modes
that have a quality factor on the order of $10^{3}$ or higher. For a CROW
studied in previous work [23], it was shown that the quality factors of the
Bloch modes varied from approximately $10^{4}$ to $10^{3}$. The relative
dielectric function does not appear in $H_{\rm NL}$, because we have
constructed it using the electric fields. If we used the displacement fields,
then the relative dielectric function would appear in the nonlinear parameter
in Eq. (III).

The dynamics of the density operator $\rho$ for the pump, signal, and idler
light can be modeled using the Lindblad master equation [24]

$\displaystyle\frac{d\rho}{dt}=-\frac{i}{\hbar}\left[H,\rho\right]+\sum_{m}\gamma_{m}\left(2b_{m}\rho
b^{\dagger}_{m}-b^{\dagger}_{m}b_{m}\rho-\rho b^{\dagger}_{m}b_{m}\right),$
(9)

where $H$ is the Hermitian part of the Hamiltonian in Eq. (4), given by

$\displaystyle H$
$\displaystyle=\sum_{m}\hbar\omega_{m}b^{\dagger}_{m}b_{m}+\sum_{m,l,p_{1},p_{2}}G_{mlp_{1}p_{2}}b^{\dagger}_{m}b^{\dagger}_{l}b_{p_{1}}b_{p_{2}}+{\rm
H.c.}.$ (10)

General expressions for the Hamiltonian and Lindblad master equation, that are
applicable to any lossy structure, can be derived using a set of non-
orthogonal quasimodes [25]. The extension from the orthogonal to the non-
orthogonal case is not straightforward. All derivations, however, in this
paper are done assuming the quasimodes are orthogonal according to Eq. (3).

###  III.1 Restrictions on the pump field and nonlinear parameter

The Hamiltonian presented in Eq. (10) is valid for a general pump field with a
quantum description. However, since we are concerned with generating squeezed
states, we let the pump be a classical field and make the undepleted pump
approximation. Furthermore, we restrict our analysis to pump fields that can
be factored into a function of space multiplied by a function of time. These
two restrictions on the pump are shown in Sec. IV to be necessary in order to
derive our solution to the Lindblad master equation.

To make the undepleted pump approximation, we let the classical pump field be
represented by a lossy multimode coherent state $\ket{\alpha(t)}$, defined as
a product of single-mode coherent states $\ket{\alpha_{p}(t)}$

$\displaystyle\ket{\alpha(t)}$
$\displaystyle\equiv\prod_{p}\ket{\alpha_{p}(t)},$ (11)

where the single-mode coherent states are defined as

$\displaystyle\ket{\alpha_{p}(t)}$ $\displaystyle=\exp(\alpha_{p}{\rm
e}^{-i\tilde{\omega}_{p}t}b^{\dagger}_{p}-{\rm H.c.})\ket{{\rm vac}}.$ (12)

For a mode $p$ with a frequency in the pump bandwidth we have

$\displaystyle b_{p}\ket{\alpha(t)}$ $\displaystyle=\alpha_{p}{\rm
e}^{-i\tilde{\omega}_{p}t}\ket{\alpha(t)},$ (13)

where $\alpha_{p}$ is the pump amplitude in the $p$th mode. The average total
photon number for the pump is

$\displaystyle\bra{\alpha(t)}\sum_{p}b^{\dagger}_{p}b_{p}\ket{\alpha(t)}=\sum_{p}|\alpha_{p}|^{2}{\rm
e}^{-2\gamma_{p}t}.$ (14)

To obtain the dynamics of the generated light only, the coherent state of the
pump is traced-out from the Lindblad master equation in Eq. (9). This amounts
to replacing the total density operator with the reduced density operator
$\bra{\alpha(t)}\rho(t)\ket{\alpha(t)}$ and replacing the nonlinear
Hamiltonian in Eq. (7) with its form in the undepleted pump approximation
[26], given by

$\displaystyle H_{\rm NL}$
$\displaystyle=\sum_{m,l,p_{1},p_{2}}G_{mlp_{1}p_{2}}\alpha_{p_{1}}\alpha_{p_{2}}{\rm
e}^{-i(\tilde{\omega}_{p_{1}}+\tilde{\omega}_{p_{2}})t}b^{\dagger}_{m}b^{\dagger}_{l}+{\rm
H.c.},$ (15)

where the pump operators $b_{p_{1}}$ and $b_{p_{2}}$ were replaced with their
the expectation value using the coherent state. We now make the crucial
assumption that the pump is in a single mode $P$. Putting $p_{1}=p_{2}\equiv
P$ into Eq. (15), we obtain

$\displaystyle H_{\rm NL}$
$\displaystyle=\alpha^{2}(t)\sum_{m,l}G_{mlPP}b^{\dagger}_{m}b^{\dagger}_{l}+{\rm
H.c.},$ (16)

where

$\displaystyle\alpha(t)$ $\displaystyle=\alpha_{P}{\rm
e}^{-i\tilde{\omega}_{P}t},$ (17)

is the time-dependent pump amplitude. For example in the CROW in Fig. 1(c),
this assumption corresponds to the pump being in a single Bloch mode. For some
systems, however, this assumption can be relaxed. If, for example, one were to
take a Gaussian pump pulse that was normal to the surface of the CROW (as was
done in Ref. [12]), then the pump would be in a continuum of free-space modes.
However, if the transverse profile of the pump inside the crystal slab does
not depend on frequency, then only a single pump mode is important to the
nonlinear interaction. Another example is where one considers a pulsed pump in
the continuum of modes of the channel waveguide in Fig. 1(a). If the duration
of the pump pulse in the channel is much longer than the ring round-trip time,
then the pump fits mostly into a single mode of the ring (as we have shown in
Ref. [27]). Since we neglect the nonlinear interaction in the channel and only
consider generation in the ring, then only a single pump mode is crucial to
the nonlinear interaction in the ring. In both of the last two examples we can
replace the constant pump amplitude $\alpha_{P}$ in Eq. (17) with a slowly-
varying temporal envelope $\alpha_{P}(t)$. The crucial point is that the
nonlinear Hamiltonian must take the form in Eq. (16) in order to derive our
solution to the Lindblad master equation.

Note that, although we are using the undepleted pump approximation, this
formalism does include linear pump loss, through the imaginary parts of the
pump mode frequencies. Note also that one can easily adapt this formalism for
SPDC, and all the results below will follow, by just replacing $\alpha(t)^{2}$
with a single pump amplitude $\alpha(t)$ and using the nonlinear parameter
$\chi^{(2)}$ for a second-order nonlinear process [12] in Eq. (16) instead.

In the next section we derive a semi-analytic solution to the Lindblad master
equation for the Hamiltonian given by Eqs. (4) and (16) .

##  IV Analytic solution to the Lindblad master equation

In previous work we studied squeezed light generation in two-mode lossy
cavities. We showed that the density operator for the generated light in the
cavity is a two-mode squeezed thermal state for all time [19]. In this work we
show that for a structure with many lossy quasimodes, the density operator for
the generated light is the multimode squeezed thermal state (MSTS),

$\rho(t)=S\left(t\right)\rho_{\rm th}(t)S^{\dagger}\left(t\right).$ (18)

Here, $S(t)$ is the unitary multimode squeezing operator given by [28]

$\displaystyle
S(t)=\exp\left(\frac{1}{2}\sum_{m,l}z_{ml}(t)b^{\dagger}_{m}b^{\dagger}_{l}-{\rm
H.c.}\right),$ (19)

where $z_{ml}(t)$ are the elements of the complex symmetric squeezing matrix
$z(t)$. The multimode thermal state $\rho_{\rm th}(t)$ is a product of single-
mode thermal states in each mode [17], such that

$\displaystyle\rho_{\rm th}(t)$
$\displaystyle=\prod_{m}\frac{1}{1+n_{m}(t)}\left(\frac{n_{m}(t)}{1+n_{m}(t)}\right)^{b^{\dagger}_{m}b_{m}},$
(20)

where $n_{m}(t)$ is defined as the average thermal photon number for the $m$th
mode. We stress that the thermal photons are not related to actual thermal
effects, but rather capture the process of photon loss. The presence of
$n_{m}(t)$ is due to a scattering process that breaks the entanglement between
the generated signal and idler photon pairs. At this point the matrix $z(t)$
and functions $n_{m}(t)$ are unknown functions of time; in Sec. IV.1 and Sec.
IV.3 we will derive equations of motion for them.

To show that the MSTS in Eq. (18) is the solution to the Lindblad master
equation in Eq. (9), we require that the equality

$\displaystyle I$ $\displaystyle=\rho_{\rm
th}^{-1/2}(t)S^{\dagger}(t)\rho(t)S(t)\rho_{\rm th}^{-1/2}(t)$ (21)

is satisfied for all time, where $I$ is the identity operator. Taking the time
derivative of Eq. (21), it can be shown that

$\displaystyle 0$ $\displaystyle=-\rho_{\rm
th}^{-1/2}S^{\dagger}\frac{dS}{dt}\rho_{\rm th}^{1/2}+\rho_{\rm
th}^{1/2}S^{\dagger}\frac{dS}{dt}\rho_{\rm th}^{-1/2}$
$\displaystyle+\rho_{\rm th}^{-1/2}S^{\dagger}\frac{d\rho}{dt}S\rho_{\rm
th}^{-1/2}+2\frac{d\rho_{\rm th}^{-1/2}}{dt}\rho_{\rm th}^{1/2},$ (22)

where we drop the time-dependence for convenience.

The majority of the rest of this section is devoted to simplifying the four
terms on the right-hand side of Eq. (IV).

Let the argument of the exponential squeezing operator be

$\displaystyle\sigma\equiv\frac{1}{2}\sum_{m,l}z_{ml}b^{\dagger}_{m}b^{\dagger}_{l}-{\rm
H.c.},$ (23)

such that $S=\exp(\sigma)$ and $S^{\dagger}=\exp(-\sigma)$. In order to
simplify Eq. (IV) we need to know the time derivative of the squeezing
operator of Eq. (19). This is not straightforward since $\sigma$ does not
commute with its time derivative. The derivative of the squeezing operator can
be written as

$\displaystyle\frac{dS}{dt}$
$\displaystyle=\frac{d}{dt}\left(1+\sigma+\frac{\sigma^{2}}{2!}+\frac{\sigma^{3}}{3!}+\ldots\right)$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\frac{\sigma^{n}\dot{\sigma}\sigma^{k}}{(n+k+1)!},$
(24)

where $\dot{\sigma}\equiv d\sigma/dt$. The sum in the last line of Eq. (IV)
has the integral representation

$\displaystyle\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\frac{\sigma^{n}\dot{\sigma}\sigma^{k}}{(n+k+1)!}$
$\displaystyle=\int_{0}^{1}dy{\rm e}^{(1-y)\sigma}\dot{\sigma}{\rm
e}^{y\sigma}.$ (25)

To prove this, we expand the exponentials in the integral in a series

$\displaystyle\int_{0}^{1}dy{\rm e}^{(1-y)\sigma}\dot{\sigma}{\rm
e}^{y\sigma}$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\frac{\sigma^{n}\dot{\sigma}\sigma^{k}}{n!k!}\int_{0}^{1}dy(1-y)^{n}y^{k}$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\frac{\sigma^{n}\dot{\sigma}\sigma^{k}}{(n+k+1)!},$
(26)

where the last line follows from the Euler integral of the first kind,

$\displaystyle\int_{0}^{1}dy(1-y)^{n}y^{k}$
$\displaystyle=\frac{n!k!}{(n+k+1)!}.$ (27)

Thus, putting Eq. (25) into Eq. (IV) we obtain

$\displaystyle\frac{dS}{dt}$ $\displaystyle=\int_{0}^{1}dy{\rm
e}^{(1-y)\sigma}\dot{\sigma}{\rm e}^{y\sigma}.$ (28)

Using the well-known Baker-Campbell-Hausdorff formula [29] on the integrand in
Eq. (28) and multiplying the equation by $S^{\dagger}=\exp(-\sigma)$ from the
left, we obtain

$\displaystyle S^{\dagger}\frac{dS}{dt}$
$\displaystyle=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(k+1)!}\mathcal{S}_{k},$
(29)

where the terms $\mathcal{S}_{k}$ can be obtained recursively from

$\displaystyle\mathcal{S}_{k+1}=[\sigma,\mathcal{S}_{k}],$ (30)

for $k=0,1,\ldots$, where

$\displaystyle\mathcal{S}_{0}$ $\displaystyle\equiv\dot{\sigma}.$ (31)

Each $\sigma$ or $\dot{\sigma}$ that appears in Eqs. (30) or (31) involves a
double sum over the quasimodes (see Eq. (23)). Thus, calculating the $k>0$
terms of $\mathcal{S}_{k}$ becomes intractable for multiple modes since each
term has $2(k+1)$ sums. However, by introducing a new Schmidt basis via a
Takagi factorization of $H_{\rm NL}$, in which the squeezing operator is
diagonal, we are able to calculate the terms $\mathcal{S}_{k}$ and show that
the sum in Eq. (29) converges even for the case of a multimode squeezing
operator (see Eq. (IV.2)). In the next subsection we introduce this new
Schmidt basis.

###  IV.1 The Schmidt basis and the diagonal multimode squeezing operator

In this subsection we diagonalize the nonlinear Hamiltonian in Eq. (16) and
the multimode squeezing operator in Eq. (19) so that we can simplify the four
terms in Eq. (IV) and prove that the MSTS is the solution to the Lindblad
master equation.

We start by performing a Takagi factorization of the nonlinear parameter $G$
in $H_{\rm NL}$. The Takagi factorization decomposes the complex symmetric
square matrix $G$ into the form

$\displaystyle G=U\Lambda U^{\rm T},$ (32)

where $\Lambda={\rm diag}(\lambda_{1},\lambda_{2},\ldots)$ is a diagonal
matrix with complex entries, $U$ is a unitary matrix with $U^{\dagger}U=1$ and
$U^{\rm T}$ is the transpose of $U$. The Takagi factorization is a special
case of symmetric singular value decomposition (SVD) [30]. In Appendix A we
show that the diagonal matrix $\Lambda$ from the Takagi factorization is just
a scaled version of the matrix of singular values obtained from the symmetric
SVD.

Putting Eq. (32) into Eq. (16) we obtain the diagonalized nonlinear
Hamiltonian

$H_{\rm NL}=\alpha^{2}(t)\sum_{\mu}\lambda_{\mu}B^{\dagger 2}_{\mu}+{\rm
H.c.},$ (33)

where $\lambda_{\mu}$ ($\mu=1,2,\ldots$) are the diagonal entries of
$\Lambda$, and we define the new creation and annihilation Schmidt operators

$\displaystyle B^{\dagger}_{\mu}$
$\displaystyle\equiv\sum_{m}U_{m\mu}b^{\dagger}_{m},$ (34) $\displaystyle
B_{\mu}$ $\displaystyle\equiv\sum_{m}U^{*}_{\mu m}b_{m},$ (35)

that have the standard commutation relation

$\displaystyle\left[B_{\mu},B^{\dagger}_{\nu}\right]=\delta_{\mu\nu}$ (36)

due to the orthogonality of the basis $(\sum_{m}U^{*}_{\mu
m}U_{m\nu}=\delta_{\mu\nu})$. We call the $B_{\mu}$ Schmidt operators, since
when there is no loss they give the Schmidt decomposition of the multimode
squeezed state. The inverse transformation is

$\displaystyle b_{m}$ $\displaystyle=\sum_{\mu}U_{m\mu}B_{\mu}.$ (37)

It is clear now why we have to assume the pump is separable into a spatial and
temporal part (see $H_{\rm NL}$ in Eq. (16)). This makes the Schmidt mode
basis independent of time. Otherwise, a new Takagi factorization would have to
be done at each time and the operators $B_{\mu}$ would become time-dependent.
In practice one calculates the nonlinear parameters $G$ for a given structure
by first calculating the quasimodes of the structure using a method such as
finite-difference time-domain (FDTD) and then calculating the spatial overlap
given in Eq. (III). Then one performs the Takagi factorization of $G$ to
obtain the matrices $U$ and $\Lambda$.

We now make the ansatz for the squeezing parameter matrix $z(t)$,

$\displaystyle z(t)=Ur(t){\rm e}^{i\phi(t)}U^{\rm T},$ (38)

where $r(t)$ and $\phi(t)$ are real diagonal matrices containing the squeezing
amplitudes and squeezing phases of the Schmidt modes, that is
$r(t)=\text{diag}(r_{1}(t),r_{2}(t),\ldots)$ and
$\phi(t)=\text{diag}(\phi_{1}(t),\phi_{2}(t),\ldots)$. Putting Eq. (38) into
Eq. (19), we obtain the diagonal multimode squeezing operator

$\displaystyle S(t)$ $\displaystyle=\prod_{\mu}S_{\mu}(t),$ (39)

where $S_{\mu}(t)$ is the single mode squeezing operator for the $\mu$th
Schmidt mode, given by

$\displaystyle S_{\mu}(t)$
$\displaystyle\equiv\exp\left(\frac{1}{2}r_{\mu}(t){\rm
e}^{i\phi_{\mu}(t)}B^{\dagger 2}_{\mu}-{\rm H.c.}\right).$ (40)

We show in Sec. IV.3 that this is the correct form of the multimode squeezing
operator, because when we impose this form, we can derive equations for the
$r_{\mu}(t)$ and $\phi_{\mu}(t)$ (see Eqs. (71) and (IV.3)) that describe a
MSTS that is a solution to the Lindblad master equation. Note, one can obtain
the squeezing matrix in the quasimode basis $z(t)$, after solving the
equations for the $r_{\mu}(t)$ and $\phi_{\mu}(t)$ and putting them into Eq.
(38).

Before moving to the next subsection, where we take the time derivative of
$S(t)$, we present the squeezing transformation for the $B_{\mu}$ and $b_{m}$
operators. These transformations are necessary when squeezing the terms in Eq.
(IV). Using the squeezing operator in Eq. (39), we obtain

$\displaystyle S^{\dagger}B_{\mu}S$
$\displaystyle=S_{\mu}^{\dagger}B_{\mu}S_{\mu}$
$\displaystyle=\cosh(r_{\mu})B_{\mu}-{\rm
e}^{i\phi_{\mu}}\sinh(r_{\mu})B^{\dagger}_{\mu},$ (41)

where the time-dependence is dropped for convenience. All other squeezing
transformations involving products of Schmidt operators (such as
$B_{\mu}B_{\nu}$) can be derived from Eq. (IV.1) and the fact that different
Schmidt operators commute. The squeezing transformation of the $b_{m}$
operators is derived with the help of the inverse transformation in Eq. (37).
For the operator $b_{m}$, we obtain

$\displaystyle S^{\dagger}b_{m}S$
$\displaystyle=\sum_{\mu}U_{m\mu}S_{\mu}^{\dagger}B_{\mu}S_{\mu}.$ (42)

###  IV.2 The time derivative of the multimode squeezing operator

In this subsection we simplify the first two terms in Eq. (IV) that involve
the derivative of the squeezing operator. We show that using the diagonal
squeezing operator in Eq. (39), we can calculate the sum in Eq. (29).

Using Eq. (39), along with the fact that the $S_{\mu}$ commute and that
$S^{\dagger}_{\mu}S_{\mu}=I$, the left-hand side of Eq. (29) can be written as

$\displaystyle S^{\dagger}\frac{dS}{dt}$
$\displaystyle=\sum_{\mu}S_{\mu}^{\dagger}\frac{dS_{\mu}}{dt}.$ (43)

Each term in the sum in Eq. (43) can be expanded using the right-hand side of
Eq. (29)

$\displaystyle S_{\mu}^{\dagger}\frac{dS_{\mu}}{dt}$
$\displaystyle=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(k+1)!}\mathcal{S}_{k\mu}.$
(44)

Putting Eq. (44) into Eq. (43), we obtain

$\displaystyle S^{\dagger}\frac{dS}{dt}$
$\displaystyle=\sum_{\mu}\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(k+1)!}\mathcal{S}_{k\mu},$
(45)

where the terms $\mathcal{S}_{k\mu}$ in Eq. (45) can be obtained recursively
from

$\displaystyle\mathcal{S}_{k+1\mu}=[\sigma_{\mu},\mathcal{S}_{k\mu}]$ (46)

for $k=0,1,\ldots$, where the $k=0$ term is defined as

$\displaystyle\mathcal{S}_{0\mu}$ $\displaystyle\equiv\dot{\sigma}_{\mu},$
(47)

where

$\displaystyle\sigma_{\mu}\equiv\frac{1}{2}r_{\mu}\exp(i\phi_{\mu})B^{\dagger
2}_{\mu}-{\rm H.c.}.$ (48)

Since the definition of $\sigma_{\mu}$ does not contain sums over Schmidt
modes, the problem is reduced to just evaluating a series of commutators
involving only single mode operators. This is straightforward to do using the
commutation relation for the $B_{\mu}$ operators. The odd $k$ and even $k$
terms of the sum are

$\displaystyle\mathcal{S}_{k\mu}$
$\displaystyle=\begin{cases}-\frac{i}{4}(2r_{\mu})^{k+1}\dot{\phi}_{\mu}\left(2B^{\dagger}_{\mu}B_{\mu}+1\right),\,\,\,\,\,{\rm
odd}\,k\geq 1\\\ \\\ \frac{i}{4}(2r_{\mu})^{k+1}\dot{\phi}_{\mu}\left({\rm
e}^{i\phi_{\mu}}B^{\dagger 2}_{\mu}+{\rm H.c.}\right),\,\,\,\,\,{\rm
even}\,k\geq 2.\end{cases}$ (49)

Putting $\dot{\sigma}_{\mu}$ and the terms in Eq. (49) into Eq. (45) and
summing over $k$, we obtain

$\displaystyle S^{\dagger}\frac{dS}{dt}$
$\displaystyle=\sum_{\mu}\Bigg{(}\frac{i}{2}\sinh^{2}(r_{\mu})\frac{d\phi_{\mu}}{dt}\left(2B^{\dagger}_{\mu}B_{\mu}+1\right)$
$\displaystyle+\left(\frac{1}{2}\frac{dr_{\mu}}{dt}+\frac{i}{4}\sinh(2r_{\mu})\frac{d\phi_{\mu}}{dt}\right){\rm
e}^{i\phi_{\mu}}B^{\dagger 2}_{\mu}$
$\displaystyle-\left(\frac{1}{2}\frac{dr_{\mu}}{dt}-\frac{i}{4}\sinh(2r_{\mu})\frac{d\phi_{\mu}}{dt}\right){\rm
e}^{-i\phi_{\mu}}B^{2}_{\mu}\Bigg{)}.$ (50)

In order to derive Eq. (IV.2) it is necessary that the Schmidt operators are
time-independent and have the standard commutator. In the next subsection we
put Eq. (IV.2) into Eq. (IV) and derive the set of coupled first-order
differential equations for the squeezing amplitudes ($r_{\mu}$), squeezing
phases ($\phi_{\mu}$), and thermal photon numbers ($n_{m}$).

###  IV.3 Differential equations for the squeezing parameters and thermal
photon numbers

We are now in a position to simplify the four terms in Eq. (IV). To make the
derivation clearer we redefine the terms in Eq. (IV) to obtain

$\displaystyle 0=T1+T2+T3,$ (51)

where

$\displaystyle T1$ $\displaystyle\equiv-\rho_{\rm
th}^{-1/2}S^{\dagger}(z)\frac{dS(z)}{dt}\rho_{\rm th}^{1/2}+\rho_{\rm
th}^{1/2}S^{\dagger}(z)\frac{dS(z)}{dt}\rho_{\rm th}^{-1/2},$ (52)
$\displaystyle T2$ $\displaystyle\equiv\rho_{\rm
th}^{-1/2}S^{\dagger}(z)\frac{d\rho}{dt}S(z)\rho_{\rm th}^{-1/2},$ (53)
$\displaystyle T3$ $\displaystyle\equiv 2\frac{d\rho_{\rm
th}^{-1/2}}{dt}\rho_{\rm th}^{1/2}.$ (54)

The transformations involving $\rho_{\rm th}^{\pm 1/2}$ in $T1$ and $T2$ are
performed using the identity

$\displaystyle\rho_{\rm th}^{\pm 1/2}B_{\mu}\rho_{\rm th}^{\mp 1/2}$
$\displaystyle=\sum_{m}U^{*}_{m\mu}b_{m}x_{m}^{\mp 1/2},$ (55)

where

$\displaystyle x_{m}\equiv\frac{n_{m}}{1+n_{m}}.$ (56)

It is simple to generalize this identity to transformations involving the
product of operators, such as $\rho_{\rm th}^{\pm 1/2}B_{\mu}B_{\nu}\rho_{\rm
th}^{\mp 1/2}$.

Our strategy is to reduce $T1$, $T2$, and $T3$ to expressions that are a sum
of Schrödinger operators multiplied by time-dependent coefficients that depend
on the squeezing amplitudes $r_{\mu}(t)$, squeezing phases $\phi_{\mu}(t)$,
and thermal photon numbers $n_{m}(t)$, respectively as well as their first
derivatives. Then we force these coefficients to be zero, such that Eq. (51)
is satisfied. The result is a set of coupled first-order differential
equations for $\dot{r}_{\mu}(t)$, $\dot{\phi}_{\mu}(t)$, and $\dot{n}_{m}(t)$.

Putting Eq. (IV.2) in $T1$ and using Eq. (55), we obtain

$\displaystyle T1$
$\displaystyle=\sum_{m,l}\left(D_{ml}b^{\dagger}_{m}b^{\dagger}_{l}+D^{*}_{ml}b_{m}b_{l}+F_{ml}b^{\dagger}_{m}b_{l}\right),$
(57)

where the time-dependent coefficients $D_{ml}$ and $F_{ml}$ are

$\displaystyle D_{ml}$
$\displaystyle=\frac{x_{m}x_{l}-1}{2\sqrt{x_{m}x_{l}}}\sum_{\mu}U_{m\mu}U_{l\mu}\left(\dot{r}_{\mu}+\frac{i}{2}\sinh(2r_{\mu})\dot{\phi}_{\mu}\right){\rm
e}^{i\phi_{\mu}},$ (58)

and

$\displaystyle F_{ml}$
$\displaystyle=i\,\frac{x_{m}-x_{l}}{\sqrt{x_{m}x_{l}}}\sum_{\mu}U_{m\mu}U^{*}_{l\mu}\dot{\phi}_{\mu}\sinh^{2}(r_{\mu}).$
(59)

Obtaining these coefficients requires many lines of algebra but it is a
straightforward exercise.

In order to derive the equations for $\dot{r}_{\mu}$ and $\dot{\phi}_{\mu}$,
we define new Schrödinger operators

$\displaystyle\hat{V}_{ml}$ $\displaystyle\equiv
b_{m}b_{l}+b^{\dagger}_{m}b^{\dagger}_{l},$ (60) $\displaystyle\hat{W}_{ml}$
$\displaystyle\equiv-ib_{m}b_{l}+ib^{\dagger}_{m}b^{\dagger}_{l},$ (61)

such that Eq. (57) can be written as

$\displaystyle T1$
$\displaystyle=\sum_{m,l}\left(\Re{D_{ml}}\hat{V}_{ml}+\Im{D_{ml}}\hat{W}_{ml}+F_{ml}b^{\dagger}_{m}b_{l}\right),$
(62)

where $\Re{D_{ml}}$ and $\Im{D_{ml}}$ are the real and imaginary parts of
$D_{ml}$.

Next we write $T2$ in terms of the operators $\hat{V}_{ml}$ and
$\hat{W}_{ml}$. The main calculation in $T2$ is the squeezing transformation
of $d\rho/dt$. This is the Lindblad master equation in Eq. (9). The nonlinear
Hamiltonian in the master equation is replaced with its diagonal form in Eq.
(33). The squeezing transformations are performed using Eq. (42). Then by
using Eq. (55), it can be shown that Eq. (53) simplifies to

$\displaystyle T2$
$\displaystyle=\sum_{m,l}\left(\Re{E_{ml}}\hat{V}_{ml}+\Im{E_{ml}}\hat{W}_{ml}+K_{ml}b^{\dagger}_{m}b_{l}\right)$
$\displaystyle+\sum_{\mu,\nu}\Gamma_{\mu\nu}\cosh(r_{\mu})\cosh(r_{\nu})\sum_{m}U_{m\mu}U^{*}_{m\nu}x_{m}$
$\displaystyle-\sum_{\mu}\Gamma_{\mu\mu}\sinh^{2}(r_{\mu}),$ (63)

where

$\displaystyle\Gamma_{\mu\nu}\equiv 2\sum_{m}\gamma_{m}U^{*}_{m\mu}U_{m\nu},$
(64)

and the time-dependent coefficients $E_{ml}$ and $K_{ml}$ are

$\displaystyle E_{ml}$ $\displaystyle=\sum_{\mu,\nu}U_{m\mu}U_{l\nu}{\rm
e}^{i\phi_{\nu}}\cosh(r_{\mu})\sinh(r_{\nu})\left(-i\Omega_{\mu\nu}\frac{x_{m}x_{l}-1}{\sqrt{x_{m}x_{l}}}+\frac{1}{2}\Gamma_{\mu\nu}\frac{1+x_{m}x_{l}-2x_{m}}{\sqrt{x_{m}x_{l}}}\right)$
$\displaystyle+\frac{i}{\hbar}\,\frac{x_{m}x_{l}-1}{\sqrt{x_{m}x_{l}}}\sum_{\mu}U_{m\mu}U_{l\mu}{\rm
e}^{i\phi_{\mu}}\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}\cosh^{2}(r_{\mu})+\alpha^{*2}\lambda^{*}_{\mu}{\rm
e}^{i\phi_{\mu}}\sinh^{2}(r_{\mu})\right),$ (65) $\displaystyle K_{ml}$
$\displaystyle=\sum_{\mu,\nu}U_{m\mu}U^{*}_{l\nu}\cosh(r_{\mu})\cosh(r_{\nu})\left(i\Omega_{\mu\nu}\frac{x_{m}-x_{l}}{\sqrt{x_{m}x_{l}}}+\frac{1}{2}\Gamma_{\mu\nu}\frac{2x_{m}x_{l}-x_{m}-x_{l}}{\sqrt{x_{m}x_{l}}}\right)$
$\displaystyle+\sum_{\mu,\nu}U_{m\mu}U^{*}_{l\nu}{\rm
e}^{i(\phi_{\mu}-\phi_{\nu})}\sinh(r_{\mu})\sinh(r_{\nu})\left(i\Omega^{*}_{\mu\nu}\frac{x_{m}-x_{l}}{\sqrt{x_{m}x_{l}}}+\frac{1}{2}\Gamma^{*}_{\mu\nu}\frac{2-x_{m}-x_{l}}{\sqrt{x_{m}x_{l}}}\right)$
$\displaystyle-\frac{2}{\hbar}\frac{x_{m}-x_{l}}{\sqrt{x_{m}x_{l}}}\sum_{\mu}U_{m\mu}U^{*}_{l\mu}\cosh(r_{\mu})\sinh(r_{\mu})\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}+\alpha^{*2}\lambda^{*}_{\mu}{\rm e}^{i\phi_{\mu}}\right),$
(66)

where

$\displaystyle\Omega_{\mu\nu}\equiv\sum_{m}\omega_{m}U^{*}_{m\mu}U_{m\nu}.$
(67)

Since the derivative of the thermal state $\rho_{\rm th}$ in Eq. (20) is easy,
$T3$ requires no special treatment in order to simplify. Putting Eq. (20) into
Eq. (54) we obtain

$\displaystyle T3$
$\displaystyle=\sum_{m}\left(-\frac{\dot{x}_{m}}{x_{m}}b^{\dagger}_{m}b_{m}+\frac{\dot{x}_{m}}{1-x_{m}}\right).$
(68)

To form the differential equations for $\dot{r}_{\mu}$ and $\dot{\phi}_{\mu}$,
we set the sum of the coefficients in front of $\hat{V}_{ml}$ and
$\hat{W}_{ml}$ in Eqs. (62) and (IV.3) to zero. This gives

$\displaystyle\Re{D_{ml}}$ $\displaystyle=-\Re{E_{ml}},$ (69)
$\displaystyle\Im{D_{ml}}$ $\displaystyle=-\Im{E_{ml}}.$ (70)

In Appendix B.1, we show that Eq. (69) leads to an equation for
$\dot{r}_{\mu}$ and Eq. (70) for $\dot{\phi}_{\mu}$. The resulting
differential equations are

$\displaystyle\dot{r}_{\mu}$
$\displaystyle=\frac{1}{i\hbar}\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}-{\rm
c.c.}\right)-\frac{1}{2}\sum_{\nu,\sigma}\cosh(r_{\nu})\sinh(r_{\sigma})\left(\Gamma_{\nu\sigma}{\rm
e}^{i(\phi_{\sigma}-\phi_{\mu})}\sum_{m,l}\frac{-n_{m}+n_{l}+1}{n_{m}+n_{l}+1}U_{m\nu}U_{l\sigma}U^{*}_{m\mu}U^{*}_{l\mu}+{\rm
c.c.}\right),$ (71) $\displaystyle\dot{\phi}_{\mu}$
$\displaystyle=2\Omega_{\mu\mu}-\frac{2}{\hbar\tanh(2r_{\mu})}\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}+{\rm c.c.}\right)$
$\displaystyle+\frac{i}{2}\sum_{\nu,\sigma}\frac{\cosh(r_{\nu})\sinh(r_{\sigma})}{\cosh(r_{\mu})\sinh(r_{\mu})}\left(\Gamma_{\nu\sigma}{\rm
e}^{i(\phi_{\sigma}-\phi_{\mu})}\sum_{m,l}\frac{-n_{m}+n_{l}+1}{n_{m}+n_{l}+1}U_{m\nu}U_{l\sigma}U^{*}_{m\mu}U^{*}_{l\mu}-{\rm
c.c.}\right).$ (72)

In order to obtain an equation for $\dot{n}_{m}$, the operators
$b^{\dagger}_{m}b_{l}$ in $T1$, $T2$, and $b^{\dagger}_{m}b_{m}$ in $T3$ are
expanded in terms of the Schmidt operator $B^{\dagger}_{\mu}B_{\nu}$ using Eq.
(37). This is necessary since it creates a common operator
$B^{\dagger}_{\mu}B_{\nu}$ that is shared between $T1$, $T2$, and $T3$. Then
setting the sum of the coefficients in front $B^{\dagger}_{\mu}B_{\nu}$ equal
to zero, we show in Appendix B.2 that this leads to the following equation for
$\dot{n}_{m}$,

$\displaystyle\dot{n}_{m}$
$\displaystyle=(1+n_{m})\sum_{\nu,\sigma}U_{m\nu}U^{*}_{m\sigma}\Gamma_{\nu\sigma}^{*}{\rm
e}^{i(\phi_{\nu}-\phi_{\sigma})}\sinh(r_{\nu})\sinh(r_{\sigma})-n_{m}\sum_{\nu,\sigma}U_{m\nu}U^{*}_{m\sigma}\Gamma_{\nu\sigma}\cosh(r_{\nu})\cosh(r_{\sigma}).$
(73)

In order to complete the derivation, we have to show that the sum of the
coefficients in front of the identity operator in $T2$ (Eq. (IV.3)) and $T3$
(Eq. (68)) is zero. Collecting the appropriate terms from Eq. (IV.3) and Eq.
(68), we obtain

$\displaystyle 0$
$\displaystyle=\sum_{m}\frac{\dot{n}_{m}}{1+n_{m}}-\sum_{\mu}\Gamma_{\mu\mu}\sinh^{2}(r_{\mu})$
$\displaystyle+\sum_{\mu,\nu}\Gamma_{\mu\nu}\cosh(r_{\mu})\cosh(r_{\nu})\sum_{m}U_{m\mu}U^{*}_{m\nu}\frac{n_{m}}{1+n_{m}}.$
(74)

Putting Eq. (73) into Eq. (IV.3) and summing over $m$, it is easily shown that
the right-hand side is equal to zero, where we use the fact that
$\sum_{m}U_{m\nu}U^{*}_{m\sigma}=\delta_{\nu\sigma}$, and
$\Gamma^{*}_{\mu\mu}=\Gamma_{\mu\mu}$. This completes the derivation since we
have shown that the coefficients in front of the operators $\hat{V}_{ml}$,
$\hat{W}_{ml}$, $B^{\dagger}_{\mu}B_{\nu}$, and the identity operator sum to
zero.

The key result of this work are Eqs. (71) - (73), which describe the time-
dependence of the squeezing amplitudes, squeezing phases, and thermal photon
numbers of the MSTS. For a system of $M$ quasimodes, they form a set of $3M$
coupled first-order differential equations, which can be easily solved on a
standard PC. The benefit of having these equations is that once they are
solved, we have the time-dependent density operator of the MSTS and with it we
can calculate any observables such as the correlation variance and photon
number. Note that the number of equations does not depend on the number of
photons at all. In contrast, numerical techniques for calculating the density
operator often require solving a large number of coupled equations that depend
on the photon number [15]. Therefore our approach can greatly reduce the
number of coupled equations, making it more feasible to study large,
multimode, lossy structures.

We note that the coupled equations as formulated apply only to a set of
discrete modes. If one has a continuum of modes, such as for the Bloch modes
of a CROW, one needs to simply discritize the modes by applying the
appropriate (e.g. periodic) boundary conditions. Otherwise, the equations
become a set of coupled integro-differential equations, which become difficult
to solve computationally without discretization. In Sec. VII we solve these
equations for a CROW with four cavities, where we use periodic boundary
conditions to quantize the allowed Bloch vectors. In Appendix C we discuss how
to solve Eqs. (71) - (73) when the initial state is the vacuum state.

In the next section we derive expressions for the expectation value of the
photon number operator and the same-time correlation variance between
different quasimodes using the MSTS.

##  V Expectation values and the correlation variance

In this section we use the MSTS to evaluate the expectation value of the
quasimode number operator and to calculate the same-time correlation variance
between pairs of quasimodes. Let $X_{m}(\varphi_{m})$ and $Y_{m}(\varphi_{m})$
be orthogonal quadrature operators defined as

$\displaystyle X_{m}(\varphi_{m})$ $\displaystyle=\frac{b_{m}{\rm
e}^{-i\varphi_{m}}+b^{\dagger}_{m}{\rm e}^{i\varphi_{m}}}{2},$ (75)
$\displaystyle Y_{m}(\varphi_{m})$ $\displaystyle=\frac{b_{m}{\rm
e}^{-i\varphi_{m}}-b^{\dagger}_{m}{\rm e}^{i\varphi_{m}}}{2i},$ (76)

where $\varphi_{m}$ is an angle in phase space. In order to quantify the
inseparability of the MSTS, the correlation variance [31] between two
quasimodes $m$ and $l$ is defined as

$\displaystyle\Delta_{ml}^{2}(t)$
$\displaystyle=\left<\left[\Delta\left(X_{m}\pm
X_{l}\right)\right]^{2}\right>+\left<\left[\Delta\left(Y_{m}\mp
Y_{l}\right)\right]^{2}\right>,$ (77)

where

$\displaystyle\left<\left[\Delta\left(X_{m}\pm X_{l}\right)\right]^{2}\right>$
$\displaystyle=\left<\left(X_{m}\pm X_{l}\right)^{2}\right>$
$\displaystyle\equiv{\rm Tr}\left[\rho(t)\left(X_{m}\pm
X_{l}\right)^{2}\right],$ (78) $\displaystyle\left<\left[\Delta\left(Y_{m}\mp
Y_{l}\right)\right]^{2}\right>$ $\displaystyle=\left<\left(Y_{m}\mp
Y_{l}\right)^{2}\right>$ $\displaystyle\equiv{\rm
Tr}\left[\rho(t)\left(Y_{m}\mp Y_{l}\right)^{2}\right]$ (79)

is defined as the variance of the operator and $\rho(t)$ is the density
operator for the MSTS. Here we have used the fact that
$\left<X_{m}\right>=\left<Y_{m}\right>=0$. It can be proved [32, 33] that the
MSTS is inseparable, and thus contains entanglement, iff $\Delta_{ml}^{2}<1$.

In order to calculate the correlation variances one needs to know the
expectation values $\left<b^{\dagger}_{m}b_{l}\right>$ and
$\left<b_{m}b_{l}\right>$. Using the Schmidt operators, we obtain

$\displaystyle\left<b^{\dagger}_{m}b_{l}\right>$
$\displaystyle=\sum_{\mu,\nu}U^{*}_{m\mu}U_{l\nu}\left<B^{\dagger}_{\mu}B_{\nu}\right>,$
(80) $\displaystyle\left<b_{m}b_{l}\right>$
$\displaystyle=\sum_{\mu,\nu}U_{m\mu}U_{l\nu}\left<B_{\mu}B_{\nu}\right>,$
(81)

where the expectation values of the Schmidt operators are

$\displaystyle\left<B^{\dagger}_{\mu}B_{\nu}\right>$ $\displaystyle={\rm
Tr}\left[\rho_{\rm th}(t)S^{\dagger}(t)B^{\dagger}_{\mu}B_{\nu}S(t)\right],$
(82) $\displaystyle\left<B_{\mu}B_{\nu}\right>$ $\displaystyle={\rm
Tr}\left[\rho_{\rm th}(t)S^{\dagger}(t)B_{\mu}B_{\nu}S(t)\right].$ (83)

These can be simplified using the squeezing transformation in Eq. (IV.1)
followed by the thermal state transformation in Eq. (55). Doing this, Eq. (80)
and Eq. (81) become

$\displaystyle\left<b^{\dagger}_{m}b_{l}\right>$
$\displaystyle=\sum_{\mu,\nu}\left(U^{*}_{m\mu}U_{l\nu}\cosh(r_{\mu})\cosh(r_{\nu})+U^{*}_{m\nu}U_{l\mu}\sinh(r_{\mu})\sinh(r_{\nu}){\rm
e}^{i(\phi_{\mu}-\phi_{\nu})}\right)\eta_{\mu\nu}+\sum_{\mu}U^{*}_{m\mu}U_{l\mu}\sinh^{2}(r_{\mu}),$
(84) $\displaystyle\left<b_{m}b_{l}\right>$
$\displaystyle=-\sum_{\mu,\nu}\left(U_{m\mu}U_{l\nu}+U_{m\nu}U_{l\mu}\right)\sinh(r_{\mu})\cosh(r_{\nu}){\rm
e}^{i\phi_{\mu}}\eta_{\mu\nu}-\sum_{\mu}U_{m\mu}U_{l\mu}\cosh(r_{\mu})\sinh(r_{\mu}){\rm
e}^{i\phi_{\mu}},$ (85)

where

$\displaystyle\eta_{\mu\nu}(t)$
$\displaystyle\equiv\sum_{m}U_{m\mu}U^{*}_{m\nu}n_{m}(t).$ (86)

The time-dependent expectation value of the photon number operator for the
$m$th quasimode $\left<b^{\dagger}_{m}b_{m}\right>$ can be obtained by letting
$m=l$ in Eq. (84). Also, putting Eq. (84) and Eq. (85) into Eq. (77) one can
calculate the time-dependent correlation variance between any two modes and
quantify the inseparability of the MSTS.

##  VI Limiting cases and comparison to other work

In this section we discuss three limiting cases to our multimode theory
presented above and compare the results to other work. The first limit is when
the modes are all lossless, such that the imaginary part of every quasimode is
equal to zero ($\gamma_{m}=0$). In this limit we show that our theory gives a
multimode squeezed vacuum state, and show that the squeezing parameter and
squeezing phase we obtain agree with other work [16]. The second limit we
consider is when there is only a single lossy mode that holds the squeezed
light. We show that in this limit our multimode theory gives a single-mode
squeezed state in agreement with previous work [18]. In the third limit we
allow two lossy modes for the squeezed light and show that our theory gives a
two-mode squeezed thermal state in agreement with our previous work [19].

###  VI.1 Lossless modes

If the modes are lossless, then the thermal photon number $n_{m}$ of each mode
is equal to zero and the thermal state in Eq. (18) gets replaced with the
vacuum state $\ket{0}\bra{0}$. Furthermore, the terms proportional to
$\gamma_{m}$ are dropped in the Lindblad master equation, so it becomes the
quantum Liouville equation. Therefore the density operator becomes

$\displaystyle\rho(t)$ $\displaystyle=S(t)\ket{0}\bra{0}S^{\dagger}(t),$ (87)

which is the density operator for a multimode squeezed vacuum state. One can
determine the squeezing parameter $r_{\mu}$ and squeezing phase $\phi_{\mu}$
for this state by letting $\gamma_{m}=0$ in Eq. (71) and Eq. (IV.3). Doing
this we obtain

$\displaystyle\dot{r}_{\mu}$
$\displaystyle=\frac{1}{i\hbar}\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}-{\rm c.c.}\right),$ (88) $\displaystyle\dot{\phi}_{\mu}$
$\displaystyle=2\Omega_{\mu\mu}-\frac{2\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}+{\rm c.c.}\right)}{\hbar\tanh(2r_{\mu})}.$ (89)

Assuming that the system starts initially in vacuum at $t=t_{i}$ with
$r_{\mu}(t_{i})=0$, the second term in Eq. (89) goes to infinity. To prevent
this we force the numerator to be zero initially and also at all later times.
Thus to obtain our final solution we put

$\displaystyle\Re{\alpha^{2}\lambda_{\mu}{\rm e}^{-i\phi_{\mu}}}$
$\displaystyle=0.$ (90)

Writing $\alpha^{2}(t)=|\alpha(t)|^{2}\exp(i\beta(t))$ and
$\lambda_{\mu}=|\lambda_{\mu}|\exp(i\theta_{\mu})$, where $\beta(t)$ is a
time-dependent phase and $\theta_{\mu}$ is a time-independent phase, we obtain
from Eq. (90) the phase condition

$\displaystyle\beta(t)+\theta_{\mu}-\phi_{\mu}(t)$ $\displaystyle=\pi/2.$ (91)

Putting Eq. (91) into Eq. (88) and Eq. (89) and integrating gives

$\displaystyle r_{\mu}(t)$
$\displaystyle=\frac{2|\lambda_{\mu}|}{\hbar}\int_{t_{i}}^{t}dt^{\prime}|\alpha(t^{\prime})|^{2},$
(92) $\displaystyle\phi_{\mu}(t)$
$\displaystyle=2(t-t_{i})\Omega_{\mu\mu}+\phi_{\mu}(t_{i}).$ (93)

Therefore the multimode squeezing operator can be written as

$\displaystyle S(t)$
$\displaystyle=\prod_{\mu}\exp\left(\frac{|\lambda_{\mu}|}{\hbar}r_{\mu}(t){\rm
e}^{-2i\Omega_{\mu\mu}t}B^{\dagger 2}_{\mu}-{\rm H.c.}\right),$ (94)

where we let $\phi_{\mu}(t_{i})=2\Omega_{\mu\mu}t_{i}$. The squeezed vacuum
state we obtain is given by $S(t)\ket{0}$. This is the same state obtained by
Quesada N. et al., following a similar procedure that uses the Takagi
factorization (see Eq. (204) and Eqs. (237) - (238) in Ref. [16]). We note
that this state has the same form one obtains in the weak pump limit (i.e.
$\alpha\ll 1$) by keeping only the first-order terms in the Dyson or Magnus
expansion of the evolution operator, a result that is also noted by Quesada N.
et al. [16]. For lossless modes the expectation values of Eq. (84) and Eq.
(85) become

$\displaystyle\left<b^{\dagger}_{m}b_{l}\right>$
$\displaystyle=\sum_{\mu}U^{*}_{m\mu}U_{l\mu}\sinh^{2}(r_{\mu}),$ (95)
$\displaystyle\left<b_{m}b_{l}\right>$
$\displaystyle=-\sum_{\mu}U_{m\mu}U_{l\mu}\cosh(r_{\mu})\sinh(r_{\mu}){\rm
e}^{i\phi_{\mu}}.$ (96)

These are the same results that were derived by Quesada N. et al. (see Eqs.
(234) - (235) in Ref. [16]).

###  VI.2 Lossy single-mode squeezed states

In this subsection we show that for a single lossy mode the coupled
differential equations Eqs. (71) - (73) give a squeezing amplitude, squeezing
phase, and thermal photon number that agrees with previous work on single-mode
squeezed thermal states [18]. In this case, the Takagi factorization of the
nonlinear parameter is trivial since there is only one mode. The matrix $U$
has only a single entry $U_{11}=1$. The single-mode squeezing operator can be
written as

$\displaystyle S(z(t))$
$\displaystyle=\exp\left(\frac{1}{2}z_{11}(t)b_{1}^{\dagger 2}-{\rm
H.c.}\right),$ (97)

where the single squeezing parameter $z_{11}(t)$ is given simply by

$\displaystyle z_{11}(t)=r_{1}(t){\rm e}^{i\phi_{1}(t)}.$ (98)

Thus, using Eq. (IV.1) the quasimode and Schmidt operators are identical
$b=B$. Only keeping the $m=l=\mu=\nu=1$ terms in Eqs. (71) - (73) and putting
$U_{11}=1$ we obtain

$\displaystyle\dot{r}_{1}$
$\displaystyle=\frac{2|\alpha|^{2}|\lambda|}{\hbar}-\frac{2\gamma_{1}}{2n_{1}+1}\cosh(r_{1})\sinh(r_{1}),$
(99) $\displaystyle\dot{\phi}_{1}$ $\displaystyle=2\omega_{1},$ (100)
$\displaystyle\dot{n}_{1}$
$\displaystyle=2\gamma_{1}\left(\sinh^{2}(r_{1})-n_{1}\right),$ (101)

where we have used the phase condition in Eq. (91) again to simplify the
equations. These equations are derived for SFWM, but they can be easily
adapted for SPDC. The only things that need to change are in Eq. (99),
$|\lambda|$ has to be redefined for SPDC, and $|\alpha|^{2}$ becomes
$|\alpha|$. Therefore, the structure of the equations is unchanged for SPDC.
In light of this, these equations are seen to agree with previous work for
single-mode squeezed states in a lossy cavity [18].

For a single lossy mode, the expectation values in Eq. (84) and Eq. (85)
become

$\displaystyle\left<b_{1}^{\dagger}b_{1}\right>$
$\displaystyle=n_{1}\cosh(2r_{1})+\sinh^{2}(r_{1}),$ (102)
$\displaystyle\left<b_{1}b_{1}\right>$
$\displaystyle=-\left(2n_{1}+1\right)\cosh(r_{1})\sinh(r_{1}).$ (103)

These are the same equations that were derived in previous work [17].

###  VI.3 Lossy two-mode squeezed states

In this subsection, we show that for two lossy modes the coupled equations
give results that agree with our previous work on two-mode squeezed states.
For two modes, the nonlinear parameter is a $2\times 2$ symmetric matrix and
its Takagi factorization is

$\displaystyle\begin{bmatrix}0&G_{12}\\\ G_{12}&0\end{bmatrix}$
$\displaystyle=U\begin{bmatrix}G_{12}&0\\\ 0&G_{12}\end{bmatrix}U^{\rm{T}},$
(104)

where we neglect generation of squeezed light in single modes by letting
$G_{11}=G_{22}=0$. The unitary matrix $U$ is given by

$\displaystyle U$ $\displaystyle=\frac{1}{2}\begin{bmatrix}1-i&1+i\\\
1+i&1-i\end{bmatrix}.$ (105)

Using Eq. (38), the squeezing parameter is given by

$\displaystyle\begin{bmatrix}0&z_{12}(t)\\\ z_{12}(t)&0\end{bmatrix}$
$\displaystyle=U\begin{bmatrix}r(t){\rm e}^{i\phi(t)}&0\\\ 0&r(t){\rm
e}^{i\phi(t)}\end{bmatrix}U^{\rm{T}},$ (106)

where $z_{12}\equiv r\exp(i\phi)$. There is a single squeezing amplitude
($r(t)$) and squeezing phase ($\phi(t)$), but the thermal photon number for
the two modes, $n_{1}(t)$ and $n_{2}(t)$, are allowed to be different. This
means that the photon loss rates are different for each mode. Putting this $U$
into Eqs. (71) - (73) and letting $\lambda=G_{12}$, we obtain

$\displaystyle\dot{r}$
$\displaystyle=\frac{2|\alpha|^{2}|\lambda|}{\hbar}-\frac{\cosh(r)\sinh(r)}{1+n_{1}+n_{2}}$
$\displaystyle\times\left[\gamma_{1}+\gamma_{2}+\left(\gamma_{1}-\gamma_{2}\right)\left(n_{2}-n_{1}\right)\right],$
(107) $\displaystyle\dot{\phi}$
$\displaystyle=\left(\omega_{1}+\omega_{2}\right),$ (108)
$\displaystyle\dot{n}_{1}$
$\displaystyle=2n_{1}\left[\gamma_{2}\sinh^{2}(r)-\gamma_{1}\cosh^{2}(r)\right]+2\gamma_{2}\sinh^{2}(r),$
(109) $\displaystyle\dot{n}_{2}$
$\displaystyle=2n_{2}\left[\gamma_{1}\sinh^{2}(r)-\gamma_{2}\cosh^{2}(r)\right]+2\gamma_{1}\sinh^{2}(r).$
(110)

These equations have the same form as those that were derived in our previous
work on two-mode squeezed states (see Eqs. (19) - (22) in Ref. [19]) generated
via SPDC. As mentioned in the previous subsection, we can adapt Eqs. (VI.3) -
(110) for SPDC by replacing $|\alpha|^{2}$ and $\lambda$ with $|\alpha|$ and
the nonlinear parameter for SPDC.

##  VII Results for a four-cavity CROW

In this section we present our results from solving the coupled equations Eqs.
(71) - (73) for SFWM in a system of four coupled-cavities. The structure is
the CROW in Fig. 1(c), but we take the structure to have only four cavities
and we enforce periodic boundary conditions at the two ends of the four-cavity
system. Thus, the quasimodes are Bloch modes, but the allowed Bloch vectors,
$k$, are quantized, such that there are only four allowed values. Although
such boundary conditions are not physically realizable, we choose this system
because it is simple and gives analytic expressions for the complex mode
frequencies and the effective nonlinear parameters. This toy-model is intended
to demonstrate how to go about solving the system of equations when there are
more than the two modes that have been modeled in previous work [19].

To obtain our results we use parameters from a CROW studied in previous work.
The photonic crystal slab that contains the coupled-cavities has a refractive
index of $n=3.4$ and a second-order refractive index of $n_{2}=4.5\times
10^{-18}\,\text{m}^{2}/\text{W}$, corresponding to Si at telecom wavelengths
[34]. The photonic crystal is a square lattice with lattice constant, $d$. The
slab thickness is $0.8d$ and the air-hole radius is $0.4d$. Using these
parameters, with find from finite-difference time domain calculations that
each individual cavity contains a single mode at the frequency
$\tilde{\omega}_{0}=(0.305-i7.71\times 10^{-6})2\pi c/d$, resulting in a
quality factor $Q_{0}=19800$ [23]. The Bloch modes and their frequencies are
obtained using a nearest-neighbour tight-binding model, where the individual
cavity modes form a basis [13]. The computed frequencies of the 4 Bloch modes
at $kD=\\{-\pi/2,0,\pi/2,\pi\\}$ are shown in Tab. 1(a), where $D\equiv 2d$ is
the periodicity of the CROW. The expression for the CROW dispersion is
approximately given by [12]

$\displaystyle\tilde{\omega}_{k}$
$\displaystyle\approx\tilde{\omega}_{0}\left(1-\tilde{\beta}_{1}\cos(kD)\right),$
(111)

where $\tilde{\beta}_{1}=9.87\times 10^{-3}-i1.97\times 10^{-5}$ is the
complex coupling parameter between individual nearest-neighbour cavity modes
[23]. This expression gives the loss dispersion of the CROW as well, it can be
obtained by $-2\Im{\tilde{\omega}_{k}}$.

We let the pump be a continuous-wave with frequency $\omega_{P}$, modelled by
a coherent state with amplitude

$\displaystyle\alpha(t)=\left|\alpha\right|{\rm e}^{-i\omega_{P}t},$ (112)

where $|\alpha|^{2}$ is the average pump photon number. The pump is in the
single Bloch mode with wavevector $k_{P}=\pi/(2D)$ and frequency
$\omega_{P}=0.305(2\pi c/d)$. Signal and idler photons are generated via SFWM
into modes with wavevectors $k_{1}$ and $k_{2}$, respectively. The effective
nonlinear parameter for this process (see Eq. (III)) has been shown to be
approximately given by [13]

$\displaystyle G_{k_{1}k_{2}k_{P}k_{P}}\approx G_{0}{\rm e}^{-i\Delta
kD/2}\frac{{\rm sinc}\left(M\Delta kD/2\right)}{{\rm sinc}\left(\Delta
kD/2\right)},$ (113)

where $\Delta k\equiv k_{1}+k_{2}-2k_{P}$ is the phase mismatch between the
different Bloch modes, $M$ is the total number of cavities, and

$\displaystyle G_{0}\equiv\frac{\hbar^{2}\omega_{0}^{2}\chi^{(3)}_{\rm
eff}}{16\epsilon_{0}(2N+1)V_{\rm eff}},$ (114)

where $\chi^{(3)}_{\rm eff}$ is the effective nonlinear coefficient, and
$V_{\rm eff}$ is the effective mode volume of the individual cavity mode.

Table 1: (a) Bloch mode wavevectors $k$ and their complex frequencies
$\tilde{\omega}_{k}$ for a CROW with 4 cavities. Here $D\equiv 2d$ is twice
the lattice spacing $d$. (b) Diagonal values $\lambda_{\mu}$ and Schmidt
frequencies $\Omega_{\mu\mu}$ for the $\mu$th Schmidt mode. The four Schmidt
modes are obtained by taking the Takagi factorization of the nonlinear
parameter in Eq. (113)

(a) $kD$ | $\tilde{\omega}_{k}d/(2\pi c)$
---|---
$\pi$ | $0.308010-i1.38\times 10^{-5}$
$\pm\pi/2$ | $0.305000-i7.71\times 10^{-6}$
0 | $0.301990-i1.63\times 10^{-6}$

(b) $\mu$ | $\lambda_{\mu}$ | $\Omega_{\mu\mu}d/(2\pi c)$
---|---|---
1 | 1.21 | 0.304877
2 | 1.16 | 0.304999
3 | 0.742 | 0.304953
4 | 0.665 | 0.305170

Now we perform a Takagi factorization of the nonlinear parameter in Eq. (113).
The resulting diagonal values $\lambda_{\mu}$ and Schmidt mode frequencies
$\Omega_{\mu\mu}$ (see Eq. (67)) are shown in Tab. 1(b). The frequencies
$\Omega_{\mu\mu}$ are all within $0.06\%$ of the pump frequency $\omega_{P}$.
This is due to the dispersion of the Bloch mode frequencies having a small
bandwidth and choosing the pump to be at the center of the band. This will not
be the case for a general structure. As we show below, however, the Schmidt
modes with frequencies that are nearly on resonance with the pump have larger
squeezing amplitudes than those that are off resonance.

To obtain a particular solution, we initially let the system be in the vacuum
state by setting all the squeezing amplitudes and thermal photon numbers equal
to zero. In Appendix C we discuss how to solve the coupled equations starting
from the vacuum, particularly how to choose the initial squeezing phase. At
time $t=0$ the pump is put into the Bloch mode $k_{P}=\pi/(2D)$ and then it
generates signal and idler photons in the Bloch modes $k_{1}$ and $k_{2}$,
respectively. Using these initial conditions, we solve Eqs. (71) - (73)
numerically using a Runge-Kutta method (such as the ode45 function in MATLAB).
For 4 modes it takes on the order of 10 seconds to solve the equations using a
standard PC.

In what follows, we use $t_{\rm c}$ as our unit of time, where $t_{\rm c}$ is
defined as the time it takes a pulse of light with group velocity $v$ to
propagate the total distance, $L$, of the coupled-cavity structure:

$\displaystyle t_{\rm c}\equiv\frac{L}{v}.$ (115)

For 4 cavities with a lattice constant of $d=480$nm, the length of the
structure is $L=2.9\mu$m. If the pulse is centered at $kD=\pi/2$, and assuming
linear dispersion in Eq. (111), then the group velocity of the pulse is
approximately $v=c/26.6$. Therefore, we obtain $t_{\rm c}=0.25$ps.
Additionally, we define the pumping strength dimensionless parameter,

$\displaystyle g\equiv\frac{4G_{0}\left|\alpha\right|^{2}t_{\rm c}}{\hbar},$
(116)

which scales the squeezing amplitude in all the Schmidt modes. In all that
follows, we set $g=1/12$. This value of $g$ can be achieved with an average
pump photon number of $|\alpha|^{2}=4.6\times 10^{7}$ and an effective mode
volume of $V_{\rm eff}=3(\mu{\rm m})^{3}$ for the cavity mode with frequency
$\tilde{\omega}_{0}$. It gives squeezing amplitudes in the Schmidt modes that
are on the order of one, while keeping the thermal photon number many orders
of magnitude below the number of photons in the pump.

We start by considering the Schmidt mode squeezing amplitudes in Fig. 2(a). At
$t=0$ the pump is turned on and the $r_{\mu}$ initially increase linearly with
time, where it can be shown using Eq. (71) that the slope is approximately
given by $\dot{r}_{\mu}\approx g|\lambda_{\mu}|/2$ (see also Eq. (144)). This
is because for short times the thermal noise and detuning
$|\Omega_{\mu\mu}-\omega_{P}|$ can be neglected. The Schmidt modes that are
detuned from the pump will have oscillations in their $r_{\mu}$ and an
$r_{\mu}$ that is smaller than that of the modes that are on resonance with
the pump. The oscillations exist for Schmidt modes that have a detuning from
the pump frequency that satisfies $|\Omega_{\mu\mu}-\omega_{P}|t_{\rm
c}>g|\lambda_{\mu}|/2$. Thus, these oscillations will disappear if we increase
the pumping strength $g$. If, however, $g$ is small, then the period of the
oscillations are approximately given by
$\pi(|\Omega_{\mu\mu}-\omega_{P}|t_{\rm c})^{-1}$. For example, using this
expression, the squeezing amplitude in the $\mu=1$ and $\mu=4$ Schmidt modes
have periods of approximately $26t_{\rm c}$ and $19t_{\rm c}$, respectively.
As time increases, the amplitude of the oscillations are dampened due to the
intrinsic losses of the Bloch modes.

Next, we consider the derivative of the Schmidt mode squeezing phases,
$\dot{\phi}_{\mu}$, in Fig. 2(b). Initially the derivative is approximately
constant, corresponding to the phases increasing linearly with time. In
Appendix C we show that for short times
$\dot{\phi}_{\mu}(t)\approx\Omega_{\mu\mu}+\omega_{P}$. The modes $\mu=2$ and
$\mu=3$ that are close to resonance with the pump have an approximately
constant phase derivative approaching $2\omega_{P}$. The detuned modes $\mu=1$
and $\mu=4$ have peaks in their phase derivatives whenever the squeezing
amplitude is close to zero. This is because in Eq. (IV.3) the term
proportional to $1/\tanh(2r_{\mu})$ is large when $r_{\mu}\ll 1$, which causes
$\dot{\phi}_{\mu}$ to increase.

(a)

(b)

Figure 2: (a) Squeezing amplitudes $r_{\mu}$ and (b) derivatives of the
squeezing phases $\dot{\phi}_{\mu}$ in the four Schmidt modes of four coupled-
cavities. Here $t_{\rm c}=0.25$ps is the time for a light pulse to cross the
length of the structure.

(a)

(b)

Figure 3: (a) Thermal photon number $n_{k}$ and (b) average total photon
number in the Bloch modes of four coupled-cavities.

Now, we present the average thermal photon numbers in the Bloch modes,
$n_{k}$, in Fig. 3(a). At short times $t<10t_{\rm c}$ the $n_{k}$ are
negligible. The $n_{k}$ do not go above approximately $10$ for the pumping
strength $g=1/12$. The majority of the thermal noise is generated in the Bloch
modes with positive wavevectors. The $n_{k}$ in the Bloch mode with
$kD=-\pi/2$ is only approximately $0.01$.

In Fig. 3(b) we calculate the average total photon number using Eq. (81).
There are on the order of approximately $10^{3}$ total photons generated by
the nonlinear process. The majority are generated in the forward-propagating
Bloch mode, $kD=\pi/2$ and $kD=\pi$ and $k=0$. Only approximately $0.1$
photons, on average, are generated in the backward-propagating Bloch mode,
$kD=-\pi/2$.

##  VIII Conclusion

The key finding in this paper is that the solution to the Lindblad master
equation for a set of lossy quasimodes is the density operator for a multimode
squeezed thermal state (see Eq. (18)). In order to prove this we introduce an
orthogonal Schmidt basis via the Takagi factorization in Sec. IV.1 to
diagonalize the nonlinear Hamiltonian and squeezing operator. The main result
is a set of coupled first-order differential equations that the squeezing
amplitude, squeezing phase, and thermal photon numbers must obey in order for
the multimode squeezed thermal state to be a solution (see Eqs. (71) - (73)).

In order to derive the solution, we make the undepleted pump approximation on
the nonlinear Hamiltonian and assume that the nonlinear parameter is
essentially the same for all pump modes in the pump bandwidth and neglect its
dependence on the pump frequencies (see Sec. III.1). The latter assumption is
valid for a pump pulse in a waveguide or free space as long as the pulse is
long in duration.

Our theory is applicable to the orthogonal set of discrete lossy modes in
structures, such as ring resonators, coupled-ring resonators, CROWs, and
high-Q cavities coupled to a waveguide in a photonic crystal (see Fig. 1).

Our results are consistent with previous work done on single-mode and two-mode
squeezed thermal states (see Refs. [18] and [19]) as presented in Sec. VI.
Also in the limit where the quasimodes become lossless, our solution reduces
to a multimode squeezed vacuum state, and it agrees with other work [16] that
used a similar approach also using a Takagi factorization.

To our knowledge, this is the first time the analytic solution for the quantum
state of light generated via nonlinear processes in $M>2$ lossy modes has been
derived. Not only is the solution of theoretical importance, it can greatly
reduce the number of coupled differential equations required to solve for the
density operator. One alternative way to determine the density operator
numerically, is by calculating its matrix elements using Fock states. The
number of possible states for $N$ photons in $M$ quasimodes is
$\sum_{j=0}^{N}\binom{j+M-1}{j}=(M+N)!/(M!N!)$ [15], where
$\binom{n}{k}=n!/[k!(n-k)!]$ is the binomial coefficient. In order to find the
matrix elements, the number of coupled equations one has to solve is the
square of this. However, using our results one can determine the density
operator by solving only $3M$ coupled first-order differential equations that
are independent of $N$. For example if $N=20$ and $M=4$ there are $10626$
basis states and thus $(10626)^{2}$ coupled equations. But our theory would
only require solving $12$ coupled equations.

The expressions we derive for the correlation variances (see Sec. V) only
contain double sums over the Schmidt modes and do not require one to solve any
additional coupled equations. Thus, the differential equations only need to be
solved once and then any same-time correlation variances of interest can be
quickly calculated.

We believe that our new solution to this important problem will make it more
feasible to study large multimode lossy structures and to optimize them for a
wide variety of quantum information applications. In future work, we will
apply our theory to more physically realistic structures, such as many-cavity
CROWs and ring resonators. We then hope to show that the generated multimode
squeezed thermal state is an approximate Gaussian cluster state [35, 7]. To do
this, one needs to determine under what conditions on the nonlinear parameter
and/or the squeezing matrix, does the state satisfy nullifier equations. If it
is possible to determine these conditions for the multimode squeezed thermal
state we derived, then it can be used as a cluster state resource for quantum
computations [8].

## Acknowledgements

This work was supported by Queen’s University and the Natural Sciences and
Engineering Research Council of Canada (NSERC).

## References

  * Hamilton  _et al._ [2017] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian boson sampling, Phys. Rev. Lett. 119, 170501 (2017).
  * Kruse  _et al._ [2019] 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).
  * Braunstein and van Loock [2005] S. L. Braunstein and P. van Loock, Quantum information with continuous variables, Rev. Mod. Phys. 77, 513 (2005).
  * Takeda and Furusawa [2019] S. Takeda and A. Furusawa, Toward large-scale fault-tolerant universal photonic quantum computing, APL Photonics 4, 060902 (2019).
  * Pfister [2020] O. Pfister, Continuous-variable quantum computing in the quantum optical frequency comb, J. Phys. B: At. Mol. Opt. Phys. 53, 012001 (2020).
  * Zhang and Braunstein [2006] J. Zhang and S. L. Braunstein, Continuous-variable gaussian analog of cluster states, Phys. Rev. A 73, 032318 (2006).
  * Menicucci  _et al._ [2011] N. C. Menicucci, S. T. Flammia, and P. van Loock, Graphical calculus for gaussian pure states, Phys. Rev. A 83, 042335 (2011).
  * Wu  _et al._ [2020] B. H. Wu, R. N. Alexander, S. Liu, and Z. Zhang, Quantum computing with multidimensional continuous-variable cluster states in a scalable photonic platform, Phys. Rev. Res. 2, 023138 (2020).
  * Vernon and Sipe [2015] Z. Vernon and J. E. Sipe, Strongly driven nonlinear quantum optics in microring resonators, Phys. Rev. A 92, 033840 (2015).
  * Yang  _et al._ [2008] 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).
  * Sharapova  _et al._ [2020] P. R. Sharapova, G. Frascella, M. Riabinin, A. M. Pérez, O. V. Tikhonova, S. Lemieux, R. W. Boyd, G. Leuchs, and M. V. Chekhova, Properties of bright squeezed vacuum at increasing brightness, Phys. Rev. Res. 2, 013371 (2020).
  * Seifoory  _et al._ [2019] H. Seifoory, L. G. Helt, J. E. Sipe, and M. M. Dignam, Counterpropagating continuous-variable entangled states in lossy coupled-cavity optical waveguides, Phys. Rev. A 100, 033839 (2019).
  * Dezfouli and Dignam [2017] M. K. Dezfouli and M. M. Dignam, Photon-pair generation in lossy coupled-resonator optical waveguides via spontaneous four-wave mixing, Phys. Rev. A 95, 033815 (2017).
  * Helt and Quesada [2020] L. Helt and N. Quesada, Degenerate squeezing in waveguides: A unified theoretical approach, J. Phys. Photonics 2, 035001 (2020).
  * Gregg and Richter [2016] M. Gregg and M. Richter, Efficient and exact numerical approach for many multi-level systems in open system CQED, New J. Phys. 18, 043037 (2016).
  * [16] N. Quesada, L. G. Helt, M. Menotti, M. Liscidini, and J. E. Sipe, Beyond photon pairs: Nonlinear quantum photonics in the high-gain regime, arXiv:2110.04340 .
  * Kim  _et al._ [1989] M. S. Kim, F. A. M. de Oliveira, and P. L. Knight, Properties of squeezed number states and squeezed thermal states, Phys. Rev. A 40, 2494 (1989).
  * Seifoory  _et al._ [2017] H. Seifoory, S. Doutre, M. Dignam, and J. E. Sipe, Squeezed thermal states: the result of parametric down conversion in lossy cavities, J. Opt. Soc. Am. B 34, 1587 (2017).
  * Vendromin and Dignam [2021] C. Vendromin and M. M. Dignam, Continuous-variable entanglement in a two-mode lossy cavity: An analytic solution, Phys. Rev. A 103, 022418 (2021).
  * Fussell and Dignam [2008] D. P. Fussell and M. M. Dignam, Quasimode-projection approach to quantum-dot–photon interactions in photonic-crystal-slab coupled-cavity systems, Phys. Rev. A 77, 053805 (2008).
  * Dignam  _et al._ [2006] M. M. Dignam, D. P. Fussell, M. J. Steel, C. Martijn de Sterke, and R. C. McPhedran, Spontaneous emission suppression via quantum path interference in coupled microcavities, Phys. Rev. Lett. 96, 103902 (2006).
  * Quesada  _et al._ [2020] 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).
  * Fussell and Dignam [2007] D. P. Fussell and M. M. Dignam, Engineering the quality factors of coupled-cavity modes in photonic crystal slabs, Appl. Phys. Lett. 90, 183121 (2007).
  * Breuer and Petruccione [2002] H. P. Breuer and F. Petruccione,  _The Theoy of Open Quantum Systems_ (Oxford University Press, 2002).
  * Dignam and Dezfouli [2012] M. M. Dignam and M. K. Dezfouli, Photon-quantum-dot dynamics in coupled-cavity photonic crystal slabs, Phys. Rev. A 85, 013809 (2012).
  * Garrison and Chiao [2008] J. C. Garrison and R. Y. Chiao,  _Quantum Optics_ (Oxford University Press, 2008).
  * Vendromin and Dignam [2020] C. Vendromin and M. M. Dignam, Optimization of a lossy microring resonator system for the generation of quadrature-squeezed states, Phys. Rev. A 102, 023705 (2020).
  * Ma and Rhodes [1990] X. Ma and W. Rhodes, Multimode squeeze operators and squeezed states, Phys. Rev. A 41, 4625 (1990).
  * [29] The formula we use is ${\rm e}^{-A}B{\rm e}^{A}=B-[A,B]+\frac{1}{2!}[A,[A,B]]-\frac{1}{3!}[A,[A,[A,B]]]+\ldots$, where $A$ and $B$ are general operators.
  * Bunse-Gerstner and Gragg [1988] A. Bunse-Gerstner and W. B. Gragg, Singular value decompositions of complex symmetric matrices, JCAM 21, 41 (1988).
  * Masada  _et al._ [2015] G. Masada, K. Miyata, A. Politi, T. Hashimoto, J. L. O’Brien, and A. Furusawa, Continuous-variable entanglement on a chip, Nat. Photonics 9, 316 (2015).
  * Duan  _et al._ [2000] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Inseparability criterion for continuous variable systems, Phys. Rev. Lett. 84, 2722 (2000).
  * Simon [2000] R. Simon, Peres-Horodecki separability criterion for continuous variable systems, Phys. Rev. Lett. 84, 2726 (2000).
  * Dinu  _et al._ [2003] M. Dinu, F. Quochi, and H. Garcia, Third-order nonlinearities in silicon at telecom wavelengths, Appl. Phys. Lett. 82, 2954 (2003).
  * Zippilli and Vitali [2020] S. Zippilli and D. Vitali, Possibility to generate any gaussian cluster state by a multimode squeezing transformation, Phys. Rev. A 102, 052424 (2020).

##  Appendix A Connection between SVD and the Takagi factorization

As mentioned in the text, the Takagi factorization is a special case of the
symmetric SVD, where the diagonal matrix from the Takagi factorization is just
a scaled version of the singular values from the SVD. In this section we
derive the Takagi factorization from the SVD, and show how the diagonal values
are related to the singular values.

The SVD of the nonlinear parameter $G$ is

$\displaystyle G=U\Sigma W^{\dagger},$ (117)

where $\Sigma$ is a diagonal matrix of real and positive singular values, and
$U^{\dagger}U=1$ and $W^{\dagger}W=1$. To obtain the Takagi factorization we
define the diagonal complex matrix

$\displaystyle\Lambda=\Sigma W^{\dagger}U^{*},$ (118)

where $U^{*}$ denotes the complex conjugate of $U$. Therefore the Takagi
factorization is

$\displaystyle G$ $\displaystyle=U\Lambda U^{\rm T},$ (119)

where the SVD can be recovered by putting Eq. (118) into Eq. (119), since
$U^{*}U^{\rm T}=1$. Therefore, to obtain the Takagi factorization of $G$, one
can perform the SVD to get the matrix $U$, and multiply the singular value
matrix $\Sigma$ by the diagonal matrix $W^{\dagger}U^{*}$.

Now we prove that $W^{\dagger}U^{*}$ in Eq. (118) is diagonal. If it is
diagonal, then it must be equal to its transpose, such that
$W^{\dagger}U^{*}=U^{\dagger}W^{*}$. Taking the Hermitian conjugate of both
sides, we obtain

$\displaystyle W^{\rm T}U$ $\displaystyle=U^{\rm T}W.$ (120)

Moreover, since $G$ is a symmetric matrix we have $G^{\rm T}=G$. Taking the
transpose of Eq. (117) we obtain

$\displaystyle W^{*}\Sigma U^{\rm T}$ $\displaystyle=U\Sigma W^{\dagger},$
$\displaystyle\Sigma U^{\rm T}W$ $\displaystyle=W^{\rm T}U\Sigma,$
$\displaystyle\Sigma U^{\rm T}W$ $\displaystyle=U^{\rm T}W\Sigma,$
$\displaystyle\Sigma U^{\rm T}W\Sigma^{-1}$ $\displaystyle=U^{\rm T}W,$ (121)

where to go from the second line to the third line we used Eq. (120). The last
line of Eq. (A) is the definition of a diagonalizable matrix, where the
diagonal form of $U^{\rm T}W$ is the same as $U^{\rm T}W$. Thus, $U^{\rm T}W$
is diagonal and so is its Hermitian conjugate $W^{\dagger}U^{*}$. Therefore we
have proved that $\Lambda$ is diagonal.

##  Appendix B Details on the derivation of the coupled differential equations
for $\dot{r}_{\mu}$, $\dot{\phi}$, and $\dot{n}_{m}$

In this section we provide the details on the derivation of the coupled
differential equations Eqs. (71) -(73).

###  B.1 Equations for $\dot{r}_{\mu}$ and $\dot{\phi}_{\mu}$

To start, we add Eqs. (69) and (70), to obtain

$\displaystyle D_{ml}=-E_{ml}.$ (122)

Putting the expression for $D_{ml}$ in Eq. (58) into Eq. (122), we obtain

$\displaystyle\frac{x_{m}x_{l}-1}{2\sqrt{x_{m}x_{l}}}\sum_{\mu}U_{m\mu}U_{l\mu}\left(\dot{r}_{\mu}+\frac{i\dot{\phi}_{\mu}}{2}\sinh(2r_{\mu})\right){\rm
e}^{i\phi_{\mu}}$ $\displaystyle=-E_{ml}.$ (123)

Multiplying both sides by $\sum_{m,l}U^{*}_{\mu m}U^{*}_{\mu l}$ and using the
orthogonality relation $\sum_{m}U^{*}_{\nu m}U_{m\mu}=\delta_{\nu\mu}$, we
obtain

$\displaystyle\dot{r}_{\mu}+\frac{i\dot{\phi}_{\mu}}{2}\sinh(2r_{\mu})$
$\displaystyle=\sum_{m,l}\frac{2\sqrt{x_{m}x_{l}}E_{ml}U^{*}_{\mu m}U^{*}_{\mu
l}{\rm e}^{-i\phi_{\mu}}}{1-x_{m}x_{l}}.$ (124)

Equating the real and imaginary parts of both sides of Eq. (124) gives the
following equations for $\dot{r}_{\mu}$ and $\dot{\phi}_{\mu}$

$\displaystyle\dot{r}_{\mu}$
$\displaystyle=\sum_{m,l}\frac{\sqrt{x_{m}x_{l}}\left(E_{ml}U^{*}_{\mu
m}U^{*}_{\mu l}{\rm e}^{-i\phi_{\mu}}+E^{*}_{ml}U_{m\mu}U_{l\mu}{\rm
e}^{i\phi_{\mu}}\right)}{1-x_{m}x_{l}},$ (125) $\displaystyle\dot{\phi}_{\mu}$
$\displaystyle=\sum_{m,l}\frac{2\sqrt{x_{m}x_{l}}\left(E_{ml}U^{*}_{\mu
m}U^{*}_{\mu l}{\rm e}^{-i\phi_{\mu}}-E^{*}_{ml}U_{m\mu}U_{l\mu}{\rm
e}^{i\phi_{\mu}}\right)}{i(1-x_{m}x_{l})\sinh(2r_{\mu})}.$ (126)

Now, to simplify Eq. (125) and Eq. (126), we focus on the first term of the
sum. Using the expression for $E_{ml}$ in Eq. (IV.3), we obtain

$\displaystyle\sum_{m,l}\frac{\sqrt{x_{m}x_{l}}E_{ml}U^{*}_{\mu m}U^{*}_{\mu
l}{\rm e}^{-i\phi_{\mu}}}{1-x_{m}x_{l}}$
$\displaystyle=\sum_{\sigma,\nu}\sum_{m,l}U_{m\sigma}U_{l\nu}U^{*}_{\mu
m}U^{*}_{\mu l}{\rm
e}^{i(\phi_{\nu}-\phi_{\mu})}\cosh(r_{\sigma})\sinh(r_{\nu})\left(i\Omega_{\sigma\nu}+\frac{1}{2}\Gamma_{\sigma\nu}\frac{1+x_{m}x_{l}-2x_{m}}{1-x_{m}x_{l}}\right)$
$\displaystyle-\frac{i}{\hbar}\sum_{\nu}\sum_{m,l}U_{m\nu}U_{l\nu}U^{*}_{\mu
m}U^{*}_{\mu l}{\rm
e}^{i(\phi_{\nu}-\phi_{\mu})}\left(\alpha^{2}\lambda_{\nu}{\rm
e}^{-i\phi_{\nu}}\cosh^{2}(r_{\nu})+\alpha^{*2}\lambda^{*}_{\nu}{\rm
e}^{i\phi_{\nu}}\sinh^{2}(r_{\nu})\right)$
$\displaystyle=i\Omega_{\mu\mu}\cosh(r_{\mu})\sinh(r_{\mu})-\frac{i}{\hbar}\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}\cosh^{2}(r_{\mu})+\alpha^{*2}\lambda^{*}_{\mu}{\rm
e}^{i\phi_{\mu}}\sinh^{2}(r_{\mu})\right)$
$\displaystyle+\frac{1}{2}\sum_{\sigma,\nu}{\rm
e}^{i(\phi_{\nu}-\phi_{\mu})}\Gamma_{\sigma\nu}\cosh(r_{\sigma})\sinh(r_{\nu})\sum_{m,l}U_{m\sigma}U_{l\nu}U^{*}_{\mu
m}U^{*}_{\mu l}\frac{1+x_{m}x_{l}-2x_{m}}{1-x_{m}x_{l}}.$ (127)

Using Eq. (56) we can write Eq. (B.1) in terms of the thermal photon numbers
as

$\displaystyle\sum_{m,l}\frac{\sqrt{x_{m}x_{l}}E_{ml}U^{*}_{\mu m}U^{*}_{\mu
l}{\rm e}^{-i\phi_{\mu}}}{1-x_{m}x_{l}}$
$\displaystyle=i\Omega_{\mu\mu}\cosh(r_{\mu})\sinh(r_{\mu})-\frac{i}{\hbar}\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}\cosh^{2}(r_{\mu})+\alpha^{*2}\lambda^{*}_{\mu}{\rm
e}^{i\phi_{\mu}}\sinh^{2}(r_{\mu})\right)$
$\displaystyle-\frac{1}{2}\sum_{\sigma,\nu}{\rm
e}^{i(\phi_{\nu}-\phi_{\mu})}\Gamma_{\sigma\nu}\cosh(r_{\sigma})\sinh(r_{\nu})\sum_{m,l}U_{m\sigma}U_{l\nu}U^{*}_{\mu
m}U^{*}_{\mu l}\frac{-n_{m}+n_{l}+1}{n_{m}+n_{l}+1}.$ (128)

The complex conjugate of Eq. (B.1) is

$\displaystyle\sum_{m,l}\frac{\sqrt{x_{m}x_{l}}E^{*}_{ml}U_{\mu m}U_{\mu
l}{\rm e}^{i\phi_{\mu}}}{1-x_{m}x_{l}}$
$\displaystyle=-i\Omega_{\mu\mu}\cosh(r_{\mu})\sinh(r_{\mu})+\frac{i}{\hbar}\left(\alpha^{2}\lambda_{\mu}{\rm
e}^{-i\phi_{\mu}}\sinh^{2}(r_{\mu})+\alpha^{*2}\lambda^{*}_{\mu}{\rm
e}^{i\phi_{\mu}}\cosh^{2}(r_{\mu})\right)$
$\displaystyle-\frac{1}{2}\sum_{\sigma,\nu}{\rm
e}^{-i(\phi_{\nu}-\phi_{\mu})}\Gamma^{*}_{\sigma\nu}\cosh(r_{\sigma})\sinh(r_{\nu})\sum_{m,l}U^{*}_{m\sigma}U^{*}_{l\nu}U_{\mu
m}U_{\mu l}\frac{-n_{m}+n_{l}+1}{n_{m}+n_{l}+1},$ (129)

where we have used the fact that $\Omega^{*}_{\mu\mu}=\Omega_{\mu\mu}$. Adding
Eq. (B.1) and Eq. (B.1) together gives the equation for $\dot{r}_{\mu}$ in Eq.
(71) in the main text. Subtracting Eq. (B.1) from Eq. (B.1) and dividing by
$\frac{i}{2}\sinh(2r_{\mu})$ (see Eq. (126)) gives the equation for
$\dot{\phi}_{\mu}$ in Eq. (IV.3) in the main text.

###  B.2 Equation for $\dot{n}_{m}$

As mentioned in Sec. IV.3 we obtain an equation for $\dot{n}_{m}$ by writing
the operators $b^{\dagger}_{m}b_{m}$ and $b^{\dagger}_{m}b_{l}$ found in the
expressions $T1$, $T2$, and $T3$ in terms of the Schmidt operator
$B^{\dagger}_{\mu}B_{\nu}$, and then force the sum of the coefficients in
front of $B^{\dagger}_{\mu}B_{\nu}$ equal to zero. The appropriate term in
$T3$ ($i.e.$ Eq. (54)) becomes

$\displaystyle-\sum_{m}b^{\dagger}_{m}b_{m}\frac{\dot{x}_{m}}{x_{m}}=-\sum_{\mu,\nu}\sum_{m}U^{*}_{m\mu}U_{m\nu}\frac{\dot{x}_{m}}{x_{m}}B^{\dagger}_{\mu}B_{\nu},$
(130)

where we have used Eq. (37). The appropriate terms in $T1$ and $T2$ ($i.e.$
Eq. (62) and Eq. (IV.3)) become

$\displaystyle\sum_{m,l}F_{ml}b^{\dagger}_{m}b_{l}=\sum_{\mu,\nu}\sum_{m,l}F_{ml}U^{*}_{m\mu}U_{l\nu}B^{\dagger}_{\mu}B_{\nu},$
(131)

and

$\displaystyle\sum_{m,l}K_{ml}b^{\dagger}_{m}b_{l}=\sum_{\mu,\nu}\sum_{m,l}K_{ml}U^{*}_{m\mu}U_{l\nu}B^{\dagger}_{\mu}B_{\nu}.$
(132)

Since we require that $0=T1+T2+T3$ (see Eq. (51)), the sum of the coefficients
multiplying $B^{\dagger}_{\mu}B_{\nu}$ in Eqs. (130) - (132) must be equal to
zero

$\displaystyle\sum_{m,l}U^{*}_{m\mu}U_{l\nu}\left(F_{ml}+K_{ml}-\frac{\dot{x}_{m}}{x_{m}}\right)=0.$
(133)

Multiplying Eq. (133) by $\sum_{\mu,\nu}U^{*}_{m\nu}U_{m\mu}$ and using the
orthogonality relation $\sum_{\mu}U^{*}_{m\mu}U_{l\mu}=\delta_{ml}$, we obtain

$\displaystyle F_{mm}+K_{mm}-\frac{\dot{x}_{m}}{x_{m}}=0,$ (134)

but from Eq. (59) we have that $F_{mm}=0$. Therefore

$\displaystyle\frac{\dot{x}_{m}}{x_{m}}=K_{mm}.$ (135)

Using Eq. (IV.3) in Eq. (135), we obtain

$\displaystyle\frac{\dot{x}_{m}}{1-x_{m}}$
$\displaystyle=\sum_{\mu,\nu}U_{m\mu}U^{*}_{m\nu}\bigg{(}-x_{m}\Gamma_{\mu\nu}\cosh(r_{\mu})\cosh(r_{\nu})$
$\displaystyle+\Gamma^{*}_{\mu\nu}{\rm
e}^{i(\phi_{\mu}-\phi_{\nu})}\sinh(r_{\mu})\sinh(r_{\nu})\bigg{)}.$ (136)

Expressing Eq. (B.2) in terms of $n_{m}$ gives

$\displaystyle\frac{\dot{n}_{m}}{1+n_{m}}$
$\displaystyle=\sum_{\mu,\nu}U_{m\mu}U^{*}_{m\nu}\bigg{(}-\frac{n_{m}}{1+n_{m}}\Gamma_{\mu\nu}\cosh(r_{\mu})\cosh(r_{\nu})$
$\displaystyle+\Gamma^{*}_{\mu\nu}{\rm
e}^{i(\phi_{\mu}-\phi_{\nu})}\sinh(r_{\mu})\sinh(r_{\nu})\bigg{)},$ (137)

which is the same as Eq. (73) in the text.

##  Appendix C Solving the coupled equations, Eqs. (71) - (73)

In this section we discuss how to solve Eqs. (71) - (73) for a system that
initially starts in the vacuum state. At time $t=0$ the vacuum state is
defined by

$\displaystyle r_{\mu}(0)$ $\displaystyle=0,$ (138) $\displaystyle n_{m}(0)$
$\displaystyle=0,$ (139)

for all $\mu$ and $m$. Putting these initial conditions into Eq. (71) we
obtain

$\displaystyle\dot{r}_{\mu}(0)$
$\displaystyle=\frac{2\left|\alpha(0)\right|^{2}\left|\lambda_{\mu}\right|}{\hbar}\sin\left(-\phi_{\mu}(0)+\theta_{\mu}\right),$
(140)

where we define

$\displaystyle\alpha(t)^{2}$
$\displaystyle=|\alpha(t)|^{2}\exp(2i\omega_{P}t),$ (141)

and

$\displaystyle\lambda_{\mu}=|\lambda_{\mu}|\exp(i\theta_{\mu}),$ (142)

where $\theta_{\mu}$ is a real number. We choose the initial squeezing phase,
$\phi_{\mu}(0)$, to be

$\displaystyle\phi_{\mu}(0)=\theta_{\mu}-\frac{\pi}{2},$ (143)

such that is maximizes the squeezing amplitude at the next time-step,
$r_{\mu}(\Delta t)$:

$\displaystyle r_{\mu}(\Delta
t)=\frac{2\left|\alpha(0)\right|^{2}\left|\lambda_{\mu}\right|}{\hbar}\Delta
t+O\left(\left(\Delta t\right)^{2}\right).$ (144)

Now, let us move on to the equation for the squeezing phase, Eq. (IV.3). It is
easily shown that using Eqs. (138), (139), and (143) will result in the second
and third terms in Eq. (IV.3) being indeterminate (0/0). Thus, at $t=0$, we
write this equation as

$\displaystyle\dot{\phi}_{\mu}(0)$
$\displaystyle=2\Omega_{\mu\mu}-\zeta_{\mu},$ (145)

where we let $\zeta_{\mu}$ be the indeterminate form. We are unable to solve
these equations unless we define $\zeta_{\mu}$. We define $\zeta_{\mu}$ by
requiring that the derivatives of the squeezing phase at $t=\Delta t$ and
$t=0$ are the same,

$\displaystyle\dot{\phi}_{\mu}(\Delta t)=\dot{\phi}_{\mu}(0),$ (146)

such that initially the squeezing phase is a linear function of time. Putting
Eqs. (138), (139), (143), and (144) into Eq. (IV.3) and using the fact that
$n_{m}(\Delta t)=0$ (which can be proven by writing Eq. (73) as a difference
equation and using the initial conditions) it can be shown that

$\displaystyle\dot{\phi}_{\mu}(\Delta t)=2\omega_{P}+\zeta_{\mu}.$ (147)

Using Eq. (147) in Eq. (146), the indeterminate form $\zeta_{\mu}$ is defined
as

$\displaystyle\zeta_{\mu}=\Omega_{\mu\mu}-\omega_{P}.$ (148)

To solve Eqs. (71) - (73) we use MATLAB’s ode45 function, that is based on a
Runge-Kutta method. The initial conditions that we use for the squeezing
amplitudes, thermal photon numbers, and squeezing phases are in Eqs. (138),
(139), and (143). We have to write an additional condition in the code that
imposes the condition that at $t=0$ the derivatives of the squeezing phases
are equal to $\dot{\phi}_{\mu}(0)=\Omega_{\mu\mu}+\omega_{P}$, otherwise the
program will return a division-by-zero error (as discussed above). The
solution is sensitive to the initial squeezing phases $\phi_{\mu}(0)$. It is
crucial that they are set to precisely the values given in Eq. (143) in order
to obtain the results we present in Sec. VII. We find, however, that the
initial value of the derivative of the phase $\dot{\phi}_{\mu}(0)$ has little
impact on the final solution, since it quickly settles to the correct value,
given by $\Omega_{\mu\mu}+\omega_{P}$.