# Pump Depletion in Parametric Amplification

Wanli Xing  Centre for Quantum Computation and Communication Technology,
School of Mathematics and Physics, University of Queensland, St Lucia,
Queensland 4072, Australia  T.C.Ralph  Centre for Quantum Computation and
Communication Technology, School of Mathematics and Physics, University of
Queensland, St Lucia, Queensland 4072, Australia

###### Abstract

We derive analytic solutions for Heisenberg evolution under the trilinear
parametric Hamiltonian which are correct to second order in the interaction
strength but are valid for all pump amplitudes. The solutions allow pump
depletion effects to be incorporated in the description of parametric
amplification in experimentally relevant scenarios and the resulting new
phenomena to be rigorously described.

Introduction: Optical parametric amplification Walls_2008  is the work-horse
of quantum optics, being a source for single photons in the weak amplification
regime Hong_1986 ; JIN14  and a source of squeezed states in the strong
amplification regime Wu_1987 ; EBE13 . A simple description of this
interaction for the non-degenerate (two-mode squeezing) case is given by the
unitary Bachor_Ralph_2019

$U=\exp{-i\chi(a^{\dagger}b^{\dagger}\alpha+ab\alpha^{*})},$ (1)

where $a$ and $b$ are annihilation operators describing the squeezed modes,
$\chi$ is the interaction strength and $\alpha$ is the amplitude of the
coherent pump field. Applying this unitary to the vacuum produces the well
known two-mode squeezed state

$\ket{\lambda}_{ab}=U|0\rangle=\sqrt{1-\lambda^{2}}\sum\lambda^{n}\ket{n}_{a}\ket{n}_{b}$
(2)

where $\lambda=\tanh{\chi\alpha}$. Alternatively the Heisenberg evolution of
the annihilation operators is given by

$\displaystyle a_{0}$
$\displaystyle=U^{\dagger}aU=a\;\cosh{\chi\alpha}-ib^{\dagger}\;\sinh{\chi\alpha}$
$\displaystyle b_{0}$
$\displaystyle=U^{\dagger}bU=b\;\cosh{\chi\alpha}-ia^{\dagger}\;\sinh{\chi\alpha}.$
(3)

Being quadratic in the operators the squeezing unitary is Gaussian, i.e.
mapping Gaussian states to Gaussian states and so the first and second moments
of the Heisenberg operators and their Hermitian conjugates are sufficient to
completely characterise squeezed states Weedbrook_2012 .

Sophisticated models based on this interaction can be built that successfully
describe a large range of devices and protocols in quantum optics
Baumgartner_1979 ; Dodonov_2002 , quantum communication Ou_1992 , quantum
computing MEN06  and quantum metrology Clerk_2010 ;
hudelist_kong_liu_jing_ou_zhang_2014 . Yet at a fundamental level this
interaction is unphysical as it is not energy conserving. This is because the
pump laser is treated as a reservoir that is unaffected, i.e. undepleted, by
the interaction. Under typical experimental conditions this is a good
approximation as the efficiency of the interaction is very low, however
efficiencies are improving all the time and experiments are moving into the
regime where depletion effects cannot be neglected ALL14 ; FLO20 . Whilst full
numerical solutions have been known for many years Walls_1970 ; Drobny_1992
and have been used in theoretical studies NAT10 ; BIR20  in conjunction with
short time perturbative approaches, these rapidly become intractable when
treating realistic systems where the pump power is large.

In this work we derive Heisenberg equations of motion which include the lowest
order non-trivial corrections to the standard equations due to pump depletion
in a consistent manner that allows for large pump powers. Although non-linear
in the mode operators, our equations are straightforward to work with and
allow an exploration of the novel physics that arises and description of the
most accessible experimental signatures of pump depletion.

Heisenberg Evolution by the Trilinear Hamiltonian: The exact unitary
describing the parametric amplification process is given by

$U=\exp{-i\chi(a^{\dagger}b^{\dagger}c+abc^{\dagger})}.$ (4)

Notice that the approximation that leads us back to quadratic form in Eq.1 is
the replacement $c\to\langle c\rangle=\alpha$. We want to know the full form
of operators $a,b,c$ in the Heisenberg picture, however we no longer obtain
the simple closed form linear equations of Eq.3. Nevertheless they can be
evaluated to any desired order using the Baker-Campbell-Hausdorff formula
sakurai_napolitano_2017 ,

$\displaystyle a_{o}$ $\displaystyle=e^{G}ae^{-G}$
$\displaystyle=a+[G,a]+\tfrac{1}{2!}[G,[G,a]]+\tfrac{1}{3!}[G,[G,[G,a]]]+\dots$
(5)

The Heisenberg operators get very complicated as we include terms of higher
and higher orders. The evolved operators up to order $\chi^{8}$ are given
explicitly in Appendix A. However, we find that a brute force expansion in
orders of $\chi$ is not the most tractable approach in situations of
experimental interest. This is because, for typical experimental parameters,
not all terms with the same power of $\chi$ contribute equally when we
calculate expectation values. To see this we can write the pump operator as
$\hat{c}=\alpha+{\hat{\delta}c}$, where the expectation value of $\hat{c}$,
$\langle\hat{c}\rangle=\alpha$, is the coherent amplitude of the pump, which
we assume to be real here; $\hat{\delta c}$ is an operator representing the
noise/quantum part of the pump; and $\langle{\hat{\delta c}}\rangle=0$ by
definition. If we carry out this expansion, then, for example, the first few
terms of $a_{o}$ become:

$\displaystyle a_{o}=a-i\chi\alpha b^{\dagger}-i\chi b^{\dagger}\delta c$
$\displaystyle+\tfrac{\chi^{2}}{2!}\quantity(-ab^{\dagger}b+a\alpha^{2}+\alpha
a\delta c^{\dagger}+\alpha a\delta c+a\delta c^{\dagger}\delta c)+\dots$ (6)

$\alpha$ can be much larger than 1, whereas $\chi$ is a small number much less
than 1. We assume $\alpha\chi$ is of order unity. In this case we see that
terms like $-\tfrac{\chi^{2}}{2!}ab^{\dagger}b$, which is of order $\chi^{2}$,
will contribute much less than terms like $\tfrac{\chi^{2}}{2!}\alpha^{2}a$,
which is of order $(\alpha\chi)^{2}\sim 1$. The result is that we need to
consider both $\alpha$ and $\chi$ when doing the expansion, and the size of a
term is determined by the difference in powers between $\alpha$ and $\chi$. In
the end we wish to derive consistent Heisenberg operator equations which can
be used to evaluate first, second and third order expectation values that are
accurate to second order in $\chi$ and to all orders in $\chi\alpha$.

Hence we perform the $\hat{c}=\alpha+\hat{\delta c}$ expansion, and only keep
terms of the form $\alpha^{n}\chi^{n}$ and $\alpha^{n-1}\chi^{n}$, and ignore
any other terms (there are also no terms where the power of $\alpha$ is higher
than the power of $\chi$). Then the operator for the signal looks like

$\displaystyle a_{o}$ $\displaystyle=a-i\alpha\chi b^{\dagger}-i\chi
b^{\dagger}\delta c+\tfrac{\chi^{2}}{2!}\quantity(\alpha^{2}a+\alpha a\delta
c^{\dagger}+\alpha a\delta c)$
$\displaystyle+\tfrac{i\chi^{3}}{3!}\quantity(-\alpha^{3}b^{{\dagger}}-\alpha^{2}b^{{\dagger}}\delta
c^{{\dagger}}-2\alpha^{2}b^{{\dagger}}\delta c)$
$\displaystyle+\tfrac{\chi^{4}}{4!}\quantity(\alpha^{4}a+2\alpha^{3}a\delta
c^{{\dagger}}+2\alpha^{3}a\delta c)$
$\displaystyle+\tfrac{i\chi^{5}}{5!}\quantity(-\alpha^{5}b^{{\dagger}}-2\alpha^{4}b^{{\dagger}}\delta
c^{{\dagger}}-3\alpha^{4}b^{{\dagger}}\delta c)$
$\displaystyle+\tfrac{\chi^{6}}{6!}\quantity(\alpha^{6}a+3\alpha^{5}a\delta
c^{{\dagger}}+3\alpha^{5}a\delta c)$
$\displaystyle+\tfrac{i\chi^{7}}{7!}\quantity(-\alpha^{7}b^{{\dagger}}-3\alpha^{6}b^{{\dagger}}\delta
c^{{\dagger}}-4\alpha^{6}b^{{\dagger}}\delta c)$
$\displaystyle+\tfrac{\chi^{8}}{8!}\quantity(\alpha^{8}a+4\alpha^{7}a\delta
c^{{\dagger}}+4\alpha^{7}a\delta c)+\dots$ (7)

Collecting terms with the same operators, we obtain several expansions in
$\alpha\chi$:

$\displaystyle a_{o}$
$\displaystyle=a(1+\tfrac{\alpha^{2}\chi^{2}}{2!}+\tfrac{\alpha^{4}\chi^{4}}{4!}+\tfrac{\alpha^{6}\chi^{6}}{6!}+\tfrac{\alpha^{8}\chi^{8}}{8!}+\dots)$
$\displaystyle+a(\delta c+\delta
c^{\dagger})(\tfrac{\alpha\chi^{2}}{2!}+\tfrac{2\alpha^{3}\chi^{4}}{4!}+\tfrac{3\alpha^{5}\chi^{6}}{6!}+\tfrac{4\alpha^{7}\chi^{8}}{8!}+\dots)$
$\displaystyle-
ib^{\dagger}(\alpha\chi+\tfrac{\alpha^{3}\chi^{3}}{3!}+\tfrac{\alpha^{5}\chi^{5}}{5!}+\tfrac{\alpha^{7}\chi^{7}}{7!}+\dots)$
$\displaystyle-ib^{\dagger}\delta
c(\chi+\tfrac{2\alpha^{2}\chi^{3}}{3!}+\tfrac{3\alpha^{4}\chi^{5}}{5!}+\tfrac{4\alpha^{6}\chi^{7}}{7!}+\dots)$
$\displaystyle-ib^{\dagger}\delta
c^{\dagger}(\tfrac{\alpha^{2}\chi^{3}}{3!}+\tfrac{2\alpha^{4}\chi^{5}}{5!}+\tfrac{3\alpha^{6}\chi^{7}}{7!}+\dots).$
(8)

We see a clear pattern from each of the expansions. Assuming the pattern
persists (checked to order $\chi^{15}$), we can write the terms in each
bracket as an infinite sum, which are found to have closed form expressions.
The coefficients of each operator are:

$\displaystyle a:\qquad$
$\displaystyle\sum_{n=0}^{\infty}\tfrac{\alpha^{2n}\chi^{2n}}{(2n)!}=\cosh\alpha\chi,$
$\displaystyle a(\delta c+\delta c^{\dagger}):$
$\displaystyle\sum_{n=1}^{\infty}\tfrac{\alpha^{2n-1}\chi^{2n}n}{(2n)!}=\tfrac{\chi}{2}\sinh\alpha\chi,$
$\displaystyle-ib^{\dagger}:\qquad$
$\displaystyle\sum_{n=0}^{\infty}\tfrac{\alpha^{2n+1}\chi^{2n+1}}{(2n+1)!}=\sinh\alpha\chi,$
$\displaystyle-b^{\dagger}\delta c:\qquad$
$\displaystyle\sum_{n=0}^{\infty}\tfrac{\alpha^{2n}\chi^{2n+1}(n+1)}{(2n+1)!}$
$\displaystyle\quad=\tfrac{i\chi}{2}\cosh\alpha\chi+\tfrac{i}{2\alpha}\sinh\alpha\chi,$
$\displaystyle-b^{\dagger}\delta c^{\dagger}:\qquad$
$\displaystyle\sum_{n=1}^{\infty}\tfrac{\alpha^{2n}\chi^{2n+1}n}{(2n+1)!}$
$\displaystyle\quad=\tfrac{i\chi}{2}\cosh\alpha\chi-\tfrac{i}{2\alpha}\sinh\alpha\chi.$
(9)

These terms give us expressions that are valid up to order $\chi$. Now, we are
interested in the second order moments which will be to order $\chi^{2}$,
hence to be safe we should expand to order $\chi^{2}$ in a similar way to what
we have done when expanding to order $\chi$. This full expansion is performed
in Appendix. B However, the vast majority of the second order terms derived in
the appendix do not contribute to the expectation values at
$\mathcal{O}(\chi^{2})$. We find that for the purpose of calculating
expectation values, we may take

$\displaystyle a_{o}$ $\displaystyle=a\quantity[\cosh\chi^{\prime}+(\delta
c+\delta c^{\dagger})\tfrac{\chi}{2}\sinh\chi^{\prime}]$ $\displaystyle-
ib^{\dagger}\big{[}\sinh\chi^{\prime}+\tfrac{\chi}{2}\cosh\chi^{\prime}(\delta
c+\delta c^{\dagger})$
$\displaystyle\qquad+\tfrac{\chi}{2\chi^{\prime}}\sinh\chi^{\prime}(\delta
c-\delta c^{\dagger})\big{]}$
$\displaystyle\qquad+\chi^{2}(iAb^{\dagger}+Ba),$ (10) $\displaystyle b_{o}$
$\displaystyle=b\quantity[\cosh\chi^{\prime}+(\delta c+\delta
c^{\dagger})\tfrac{\chi}{2}\sinh\chi^{\prime}]$ $\displaystyle-
ia^{\dagger}\big{[}\sinh\chi^{\prime}+\tfrac{\chi}{2}\cosh\chi^{\prime}(\delta
c+\delta c^{\dagger})$
$\displaystyle\qquad+\tfrac{\chi}{2\chi^{\prime}}\sinh\chi^{\prime}(\delta
c-\delta c^{\dagger})\big{]}$
$\displaystyle\qquad+\chi^{2}(Bb+iAa^{\dagger}),$ (11) $\displaystyle c_{o}$
$\displaystyle=\alpha_{o}+\delta c_{o}$
$\displaystyle=\alpha-\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}+\chi^{3}C$
$\displaystyle\quad+\delta
c-(a^{\dagger}a+b^{\dagger}b)\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}$
$\displaystyle\quad-
ia^{\dagger}b^{\dagger}\tfrac{\chi}{2}(1-\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})$
$\displaystyle\quad-
iab\tfrac{\chi}{2}(1+\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})$
$\displaystyle\quad+\chi^{2}D\delta c^{\dagger}.$ (12)

where

$\displaystyle A$
$\displaystyle=\tfrac{-5\chi^{\prime}\cosh\chi^{\prime}+2\sinh\chi^{\prime}-\chi^{\prime
2}\sinh\chi^{\prime}+\sinh 3\chi^{\prime}}{8\chi^{\prime 2}},$ $\displaystyle
B$ $\displaystyle=-\tfrac{-\cosh\chi^{\prime}-\chi^{\prime
2}\cosh\chi^{\prime}+\cosh
3\chi^{\prime}-3\chi^{\prime}\sinh\chi^{\prime}}{8\chi^{\prime 2}},$
$\displaystyle C$ $\displaystyle=\tfrac{-3-4\chi^{\prime 2}+(2-4\chi^{\prime
2})\cosh 2\chi^{\prime}+\cosh 4\chi^{\prime}-2\chi^{\prime}\sinh
2\chi^{\prime}}{32\chi^{\prime 3}}$ $\displaystyle D$
$\displaystyle=-\tfrac{1-\cosh 2\chi^{\prime}+\chi^{\prime}\sinh
2\chi^{\prime}}{4\chi^{\prime 2}}.$

and $\chi^{\prime}=\alpha\chi$. These are ‘effective operators’ in the sense
that they give the correct results for all normally ordered second order
moments, i.e.
$\expectationvalue*{a_{o}^{\dagger}a_{o}},\expectationvalue*{a_{o}b_{o}},\expectationvalue*{\delta
c_{o}^{\dagger}\delta c_{o}},\expectationvalue*{\delta c_{o}\delta c_{o}}$ and
$\alpha_{o}^{2}$, and therefore all variances calculated from these operators
are correct. Terms that arise in any calculation that are not normally ordered
must first be reordered using the standard Boson commutator relations, e.g.
$[a_{o},a_{o}^{\dagger}]=1$, before proceeding with the calculation.

Eqs 10, 11 and 12 are the main results of this paper. They provide a tractable
and physically intuitive way to investigate the lowest order corrections to
the behaviour of two-mode squeezing when pump depletion becomes significant.

Expectation values: We are now in a position to investigate the physics of the
pump-depleted squeezer. Let us first consider the photon number in the pump
and the squeezed modes. Assuming the pump is initially in a coherent state and
the squeezed modes are initially in vacuum states, the photon number in the
pump after the interaction is given by

$\displaystyle\expectationvalue*{c^{\dagger}_{o}c_{o}}$
$\displaystyle=\alpha^{2}_{o}+\expectationvalue*{\delta c^{\dagger}_{o}\delta
c_{o}}$
$\displaystyle=\alpha^{2}-\sinh^{2}\chi^{\prime}+(\tfrac{\chi}{2\chi^{\prime}})^{2}\sinh^{4}\chi^{\prime}+2\chi^{2}\chi^{\prime}C$
$\displaystyle+\tfrac{\chi^{2}}{4}$
$\displaystyle(1-\tfrac{2}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime}+\tfrac{1}{\chi^{\prime
2}}\sinh^{2}\chi^{\prime}\cosh^{2}\chi^{\prime}).$ (13)

As expected the pump is now depleted by the interaction with
$\expectationvalue*{c^{\dagger}_{o}c_{o}}<\alpha^{2}$. In addition there is
now a coherent contribution to the photon number, $\alpha_{o}^{2}$ and an
incoherent contribution, $\expectationvalue*{\delta c^{\dagger}_{o}\delta
c_{o}}$. The photon numbers in the squeezed modes are given by

$\displaystyle\expectationvalue*{a^{\dagger}_{o}a_{o}}$
$\displaystyle=\expectationvalue*{b^{\dagger}_{o}b_{o}}=\sinh^{2}\chi^{\prime}-2\chi^{2}\sinh\chi^{\prime}A$
$\displaystyle+\tfrac{\chi^{2}}{4}$
$\displaystyle(\cosh^{2}\chi^{\prime}-\tfrac{2}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime}+\tfrac{1}{\chi^{\prime
2}}\sinh^{2}\chi^{\prime}).$ (14)

The photon number in the squeezed modes are also lower than that predicted by
the undepleted pump model. It is straightforward to confirm that energy
conservation now holds as

$\displaystyle\expectationvalue*{c^{\dagger}_{o}c_{o}}+{{1}\over{2}}(\expectationvalue*{a^{\dagger}_{o}a_{o}}+\expectationvalue*{b^{\dagger}_{o}b_{o}})=\alpha^{2},$
(15)

where we have taken into account that the energy of the squeezed mode photons
is half that of the pump photons.

The other non-zero expectation values up to third order can also be calculated
and give

$\displaystyle\expectationvalue*{a_{o}b_{o}}$
$\displaystyle=-\tfrac{i}{2}\sinh 2\chi^{\prime}$
$\displaystyle+\tfrac{i\chi^{2}}{16\chi^{\prime
2}}\Big{(}-4\chi^{\prime}-6\chi^{\prime}\cosh 2\chi^{\prime}$
$\displaystyle\qquad+(1-4\chi^{\prime 2})\sinh 2\chi^{\prime}+2\sinh
4\chi^{\prime}\Big{)},$ $\displaystyle\expectationvalue*{\delta c_{o}\delta
c_{o}}$
$\displaystyle=\tfrac{1}{32}\Big{(}-8\chi^{2}(1+\tfrac{1}{\chi^{\prime}}\sinh
2\chi^{\prime})$ $\displaystyle\qquad+\tfrac{\chi^{2}}{\chi^{\prime 2}}(8\cosh
2\chi^{\prime}+\cosh 4\chi^{\prime}-9)\Big{)},$
$\displaystyle\expectationvalue*{a_{o}b_{o}\delta c_{o}}$
$\displaystyle=\tfrac{i\chi}{2\chi^{\prime}}\sinh\chi^{\prime}\cosh^{3}\chi\prime-\tfrac{i\chi}{2}\cosh^{2}\chi^{\prime},$
$\displaystyle\expectationvalue*{a_{o}b_{o}\delta c^{\dagger}_{o}}$
$\displaystyle=\tfrac{i\chi}{2\chi^{\prime}}\sinh^{3}\chi^{\prime}\cosh\chi^{\prime}-\tfrac{i\chi}{2}\sinh^{2}\chi^{\prime}.$
(16)

From these we can calculate other interesting observables such as the
quadrature variances of the output pump beam. The amplitude variance is given
by

$\displaystyle V_{xc}$ $\displaystyle=\langle(\delta c_{o}+\delta
c_{o}^{\dagger})^{2}\rangle$ $\displaystyle=2\expectationvalue*{\delta
c^{\dagger}_{o}\delta c_{o}}+\expectationvalue*{\delta c_{o}\delta
c_{o}}+\expectationvalue*{\delta c_{o}^{\dagger}\delta c_{o}^{\dagger}}+1$
$\displaystyle=1-\tfrac{\chi^{2}}{\chi^{\prime}}\sinh 2\chi^{\prime}+$
$\displaystyle\;\;\;\;\tfrac{\chi^{2}}{8\chi^{\prime 2}}(-5+4\cosh
2\chi^{\prime}+\cosh 4\chi^{\prime}),$ (17)

whilst the phase variance is given by

$\displaystyle V_{pc}$ $\displaystyle=-\langle(\delta c_{o}-\delta
c_{o}^{\dagger})^{2}\rangle$ $\displaystyle=2\expectationvalue*{\delta
c^{\dagger}_{o}\delta c_{o}}-\expectationvalue*{\delta c_{o}\delta
c_{o}}-\expectationvalue*{\delta c_{o}^{\dagger}\delta c_{o}^{\dagger}}+1$
$\displaystyle=1-\chi^{2}(\tfrac{\sinh^{2}\chi^{\prime}}{\chi^{\prime 2}}-1).$
(18)

We notice that the output pump has $V_{p}<1<V_{x}$ indicating it has become
phase squeezed through the interaction. Because
$\expectationvalue*{a_{o}a_{o}}=\expectationvalue*{b_{o}b_{o}}=0$ the
quadrature variances of the signal and idler are isotropic and given by
$V_{xj}=V_{pj}=1+2\expectationvalue*{j^{\dagger}_{o}j_{o}}$ where
$j=\\{a,b\\}$. Given the phase convention we have adopted the correlations
between the signal and idler exist between orthogonal quadratures. Hence the
difference and sum squeezing between the signal and idler are given by

$\displaystyle V_{xp}^{\pm}$ $\displaystyle={{\langle(a_{o}+a_{o}^{\dagger}\pm
i(b_{o}^{\dagger}-b_{o}))^{2}\rangle}\over{2}}$
$\displaystyle=2\expectationvalue*{a^{\dagger}_{o}a_{o}}\mp
2i\expectationvalue*{a_{o}b_{o}}+1$
$\displaystyle=\frac{1}{16}(\cosh\chi^{\prime}\mp\sinh\chi^{\prime})\;\;\;\times$
$\displaystyle\quad\bigg{[}-\tfrac{3\chi^{2}}{\chi^{\prime 2}}\cosh
3\chi^{\prime}+$ $\displaystyle(\tfrac{3\chi^{2}}{\chi^{\prime
2}}\mp\tfrac{20\chi^{2}}{\chi^{\prime}}+8(2+\chi^{2}))\cosh\chi^{\prime}$
$\displaystyle\quad+2(\pm\tfrac{5\chi^{2}}{\chi^{\prime 2}}(1+\cosh
2\chi^{\prime})+$ $\displaystyle\tfrac{2\chi^{2}}{\chi^{\prime}}\mp
4(2+\chi^{2}))\sinh\chi^{\prime}\bigg{]}.$ (19)

We find $V_{xp}^{+}<1<V_{xp}^{-}$ indicating entanglement between the signal
and idler beams as expected.

In contrast to the undepleted case, neither the output pump or signal and
idler are in minimum uncertainty states. That is, we find $V_{xc}V_{pc}>1$ and
$V_{xp}^{+}V_{xp}^{-}>1$. Given the overall unitarity of the interaction this
indicates either non-Gaussianity or entanglement or both is emerging between
the pump and the signal and idler. Indeed it is both. There is no correlation
between the pump or signal or idler in the second order moments as would be
required if Gaussian entanglement was emerging. Instead we find correlation
between the pump and signal and idler in the third order moments indicating
non-Gaussian entanglement. In particular we can consider the quadrature
correlations $V_{abc}=\langle X_{ax}X_{bp}\delta X_{cx}\rangle$ and show

$\displaystyle V_{abc}=\langle
i(a_{o}+a_{o}^{\dagger})(b_{o}^{\dagger}-b_{o})(\delta c_{o}+\delta
c_{o}^{\dagger})\rangle$
$\displaystyle=-2i(\expectationvalue*{a_{o}b_{o}\delta
c_{o}}+\expectationvalue*{a_{o}b_{o}\delta c^{\dagger}_{o}})$
$\displaystyle=\tfrac{\chi}{2}(\tfrac{1}{\chi^{\prime}}\cosh\chi^{\prime}\sinh\chi^{\prime}-1)(\cosh^{2}\chi^{\prime}+\sinh^{2}\chi^{\prime}),$
(20)

The fact that this moment is non-zero (whilst all related first order moments
are zero) indicates a non-Gaussian quantum correlation, i.e. entanglement.

Figure 1: Amplitude and phase variances of the output pump as a function of
input pump amplitude $\alpha$: (a) the amplitude variance is shown in blue
(upper) trace and the phase variance is the orange (lower) trace. The red-
dashed line is the quantum noise limit. Here $\chi=0.02$. A small amount of
squeezing is seen for these parameters which are plotted using Eqs Pump
Depletion in Parametric Amplification and Pump Depletion in Parametric
Amplification; (b) the amplitude variance is shown in green (upper) trace and
the phase variance is the orange (lower) trace. Here $\chi=0.001$. The
stronger pump powers in this regime mean that now $V_{pc}=1$ and $V_{xc}$ is
given by Eqs 23.

Strong pump regime: One parameter regime which is expected to be relevant for
experimental tests of these effects is the strong pump regime. That is we take
$\alpha$ sufficiently large that $\chi^{\prime}>>1$, whilst still insisting
$\chi$ is sufficiently small that our second order expansion remains valid. We
note that although this regime is inaccessible to numerical approaches, it is
easily explored with our analytical expressions. By neglecting the negative
exponentials in our cosh and sinh terms and keeping only the largest of the
positive exponentials we can significantly simplify our expectation values.
The average photon numbers of the pump, signal and idler become:

$\displaystyle\expectationvalue*{c^{\dagger}_{o}c_{o}}$
$\displaystyle=\alpha^{2}-{\tfrac{e^{2\chi^{\prime}}}{4}}+\tfrac{\chi^{2}}{16\chi^{\prime
2}}\;\;e^{4\chi^{\prime}},$ (21)

and

$\displaystyle\expectationvalue*{a^{\dagger}_{o}a_{o}}$
$\displaystyle=\expectationvalue*{b^{\dagger}_{o}b_{o}}={\tfrac{e^{2\chi^{\prime}}}{4}}-\tfrac{\chi^{2}}{16\chi^{\prime
2}}\;\;e^{4\chi^{\prime}}.$ (22)

Figure 2: Third order correlation $V_{abc}$ as a function of input pump
amplitude $\alpha$. Here $\chi=0.001$. Significant affects are seen at
relatively low pump amplitude.

The pump amplitude quadrature variance becomes

$\displaystyle V_{xc}$ $\displaystyle=1+\tfrac{\chi^{2}}{16\chi^{\prime
2}}\;\;e^{4\chi^{\prime}},$ (23)

whilst the phase quadrature remains at the quantum noise level, $V_{pc}=1$,
given this approximation. The difference squeezing between the signal and
idler is given by

$\displaystyle V_{xp}^{+}$
$\displaystyle=e^{-2\chi^{\prime}}+\tfrac{\chi^{2}}{16\chi^{\prime
2}}\;\;e^{2\chi^{\prime}},$ (24)

whilst the sum squeezing is given by

$\displaystyle V_{xp}^{-}$
$\displaystyle=e^{2\chi^{\prime}}-\tfrac{\chi^{2}}{4\chi^{\prime
2}}\;\;e^{4\chi^{\prime}}.$ (25)

Notice this leads to the uncertainty product
$V_{xp}^{-}V_{xp}^{+}=1+{{\chi^{2}}\over{16\chi^{\prime
2}}}\;\;e^{{4\chi^{\prime}}}$ indicating the departure from a pure Gaussian
entangled state.

Perhaps surprisingly the strongest effect is seen in the third order
correlations. The quadrature correlation between phase quadrature of the idler
and the amplitude quadratures of the signal and pump becomes

$\displaystyle
V_{abc}=\tfrac{\chi}{16\chi^{\prime}}\;\;e^{4\chi^{\prime}}{-\tfrac{\chi}{4}e^{2\chi^{\prime}}}.$
(26)

As this moment is linear in $\chi$ it should be the first quantum effect to
become observable as we enter the pump depletion regime at high pump powers.

Conclusion: We have derived non-linear Heisenberg equations describing the
evolution of quantum fields through the trilinear Hamiltonian which models
parametric amplification with pump depletion. Unlike previous treatments we
perform our perturbative expansion in such a way as to allow the strong pump
regime to be explored. We expect our results to be immediately useful in
describing and motivating squeezing experiments in the strong pump regime.
Being Heisenberg picture equations they provide good intuition about the
physics and can be easily adapted to account for imperfections such as loss
and excess noise. We also expect our solutions to stimulate investigations
into novel quantum protocols and technologies which may be enabled by the non-
Gaussian correlations STR18  that emerge as we push further into the depleted
pump regime of squeezing.

Acknowledgements: This research was supported by the Australian Research
Council (ARC) under the Centre of Excellence for Quantum Computation and
Communication Technology (CE170100012).

Note added: After the completion of this work we became aware of a related,
but distinct approach to pump depletion in single-mode squeezing following a
Schrödinger picture approach YAN21

## References

  * (1) D.F.Walls and G.J.Milburn, Quantum Optics, (2nd Edition, Springer-Verlag 2008).
  * (2) C.K.Hong and L.Mandel, Experimental realization of a localized one-photon state, Phys.Rev.Lett. 56, 58 (1986).
  * (3) R.B Jin, R. Shimizu, K. Wakui, Mikio Fujiwara, T. Yamashita, S. Miki, H. Terai, Z. Wang, M. Sasaki, Pulsed Sagnac polarization-entangled photon source with a PPKTP crystal at telecom wavelength, Opt. Express 22, 11498 (2014)
  * (4) Ling-An Wu and Min Xiao and H. J. Kimble, Squeezed states of light from an optical parametric oscillator, J. Opt. Soc. Am. B, 4, 1465 (1987).
  * (5) T. Eberle, V. Händchen, R. Schnabel, Stable control of 10 dB two-mode squeezed vacuum states of light, Optics Express 21, 11546 (2013).
  * (6) H.-A. Bachor and T. C. Ralph, A guide to experiments in quantum optics (3rd Edition Wiley, New York, 2019).
  * (7) Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J. Cerf, Timothy C. Ralph, Jeffrey H. Shapiro, and Seth Lloyd, Gaussian quantum information, Rev. Mod. Phys. 82, 1155 (2012).
  * (8) R. Baumgartner and R. Byer, Optical parametric amplification, IEEE Journal of Quantum Electronics 15, 432 (1979).
  * (9) V. V. Dodonov, ‘Nonclassical’ states in quantum optics: a ‘squeezed’ review of the first 75 years, Journal of Optics B: Quantum and Semiclassical Optics 4, R1 (2002).
  * (10) Z. Y. Ou, S. F. Pereira, H. J. Kimble, and K. C. Peng, Realization of the Einstein-Podolsky-Rosen paradox for continuous variables, Phys. Rev. Lett. 68, 3663 (1992).
  * (11) Nicolas C. Menicucci, Peter van Loock, Mile Gu, Christian Weedbrook, Timothy C. Ralph, Michael A. Nielsen, Universal Quantum Computation with Continuous-Variable Cluster States, Phys. Rev. Lett. 97, 110501 (2006).
  * (12) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Introduction to quantum noise, measurement and amplification, Rev. Mod. Phys. 82, 1155 (2010).
  * (13) F. Hudelist, J. Kong, C. Liu, J. Jing, Z. Ou, and W. Zhang, Quantum metrology with parametric amplifier-based photon correlation interferometers, Nature Communications 5, 10 (2014).
  * (14) A. Allevi, O. Jedrkiewicz, E. Brambilla, A. Gatti, J. Peřina, Jr., O. Haderka, and M. Bondani, Coherence properties of high-gain twin beams, Phys. Rev. A 90, 063812 (2014).
  * (15) Jefferson Flórez, Jeff S. Lundeen, and Maria V. Chekhova, Pump depletion in parametric down-conversion with low pump energies, Optics Letters 45, 4264 (2020).
  * (16) D. F. Walls and R. Barakat, Quantum-mechanical amplification and frequency conversion with a trilinear hamiltonian, Phys. Rev. A 1, 446 (1970).
  * (17) G. Drobny and I. Jex, Quantum properties of field modes in trilinear optical processes, Phys. Rev. A 46, 499 (1992).
  * (18) P. D. Nation and M. P. Blencowe, The trilinear Hamiltonian: a zero-dimensional model of Hawking radiation from a quantized source, New Journal of Physics 12, 095013 (2010).
  * (19) S. Ding, G. Maslennikov, R. Hablutzel and D. Matsukevich, Quantum Simulation with a Trilinear Hamiltonian, Phys. Rev. Lett. 121 130502 (2018).
  * (20) Richard J. Birrittella, Paul M. Alsing, Christopher C. Gerry, Phase effects in coherently-stimulated down-conversion with a quantized pump field, Phys. Rev. A 101, 013813 (2020).
  * (21) J. J. Sakurai and J. Napolitano, Modern Quantum Mechanics, 2nd ed. (Cambridge University Press, 2017).
  * (22) I. Straka, L. Lachman, J. Hloušek, M. Miková, M. Mičuda, M. Ježek, and R. Filip, Quantum non-Gaussian multiphoton light, npj Quantum Inf. 4, 4 (2018).
  * (23) Ryotatsu Yanagimoto,et al, Onset of non-Gaussian quantum physics in pulsed squeezing with mesoscopic fields, arXiv:2111.13799

##  I Appendix

###  I.1  Full Mode Expansions to $\chi^{8}$

The exact unitary describing the parametric amplification process is given by

$\displaystyle U=\exp{-i\chi(a^{\dagger}b^{\dagger}c+abc^{\dagger})}.$

Notice that the approximation that leads us back to the quadratic form in Eq.1
of the main text is the replacement $c\to\langle c\rangle=\alpha$. We want to
know the full forms of operators $a,b,c$ in the Heisenberg picture, which we
denote $a_{o},b_{o},c_{o}$. We no longer obtain the simple closed form linear
equations of Eq.3, nevertheless these Heisenberg operators can be evaluated to
any desired order using the Baker-Campbell-Hausdorff formula. For example, the
signal mode is

$\displaystyle a_{o}=e^{G}ae^{-G}$
$\displaystyle=a+[G,a]+\tfrac{1}{2!}[G,[G,a]]+\tfrac{1}{3!}[G,[G,[G,a]]]+\dots$

In this case $G=i\chi(a^{\dagger}b^{\dagger}c+abc^{\dagger})$. For reference,
the Heisenberg operator for the signal mode $a_{o}$ to order $\chi^{8}$ is:

$\displaystyle a_{o}$ $\displaystyle=a-i\chi
b^{\dagger}c+\tfrac{\chi^{2}}{2!}\quantity(-ab^{{\dagger}}b+ac^{{\dagger}}c)$
$\displaystyle+\tfrac{i\chi^{3}}{3!}\quantity(2a^{{\dagger}}ab^{{\dagger}}c-b^{{\dagger}}c^{{\dagger}}c^{2}+b^{{\dagger}}c+b^{{\dagger}2}bc-2a^{2}bc^{{\dagger}})$
$\displaystyle+\tfrac{\chi^{4}}{4!}\quantity(4a^{{\dagger}}b^{{\dagger}2}c^{2}+4a^{{\dagger}}a^{2}b^{{\dagger}}b-4a^{{\dagger}}a^{2}c^{{\dagger}}c+ab^{{\dagger}}b-10ab^{{\dagger}}bc^{{\dagger}}c+ab^{{\dagger}2}b^{2}-7ac^{{\dagger}}c+ac^{{\dagger}2}c^{2})$
$\displaystyle+\tfrac{i\chi^{5}}{5!}\bigg{(}28a^{{\dagger}}ab^{{\dagger}}c^{{\dagger}}c^{2}-16a^{{\dagger}}ab^{{\dagger}}c-28a^{{\dagger}}ab^{{\dagger}2}bc+8a^{{\dagger}}a^{3}bc^{{\dagger}}-8a^{{\dagger}2}a^{2}b^{{\dagger}}c+25b^{{\dagger}}c^{{\dagger}}c^{2}-b^{{\dagger}}c^{{\dagger}2}c^{3}-b^{{\dagger}}c$
$\displaystyle\qquad+14b^{{\dagger}2}bc^{{\dagger}}c^{2}-3b^{{\dagger}2}bc-b^{{\dagger}3}b^{2}c+16a^{2}b^{{\dagger}}b^{2}c^{{\dagger}}+12a^{2}bc^{{\dagger}}-16a^{2}bc^{{\dagger}2}c\bigg{)}$
$\displaystyle+\tfrac{\chi^{6}}{6!}\bigg{(}44a^{{\dagger}}b^{{\dagger}2}c^{{\dagger}}c^{3}-72a^{{\dagger}}b^{{\dagger}2}c^{2}-44a^{{\dagger}}b^{{\dagger}3}bc^{2}-28a^{{\dagger}}a^{2}b^{{\dagger}}b+216a^{{\dagger}}a^{2}b^{{\dagger}}bc^{{\dagger}}c-44a^{{\dagger}}a^{2}b^{{\dagger}2}b^{2}$
$\displaystyle\qquad+68a^{{\dagger}}a^{2}c^{{\dagger}}c-44a^{{\dagger}}a^{2}c^{{\dagger}2}c^{2}-72a^{{\dagger}2}ab^{{\dagger}2}c^{2}-16a^{{\dagger}2}a^{3}b^{{\dagger}}b+16a^{{\dagger}2}a^{3}c^{{\dagger}}c-ab^{{\dagger}}b+216ab^{{\dagger}}bc^{{\dagger}}c$
$\displaystyle\qquad-91ab^{{\dagger}}bc^{{\dagger}2}c^{2}-3ab^{{\dagger}2}b^{2}+91ab^{{\dagger}2}b^{2}c^{{\dagger}}c-ab^{{\dagger}3}b^{3}+41ac^{{\dagger}}c-85ac^{{\dagger}2}c^{2}+ac^{{\dagger}3}c^{3}-40a^{3}b^{2}c^{{\dagger}2}\bigg{)}$
$\displaystyle+\tfrac{i\chi^{7}}{7!}\bigg{(}-1386a^{{\dagger}}ab^{{\dagger}}c^{{\dagger}}c^{2}+270a^{{\dagger}}ab^{{\dagger}}c^{{\dagger}2}c^{3}+98a^{{\dagger}}ab^{{\dagger}}c-1204a^{{\dagger}}ab^{{\dagger}2}bc^{{\dagger}}c^{2}+598a^{{\dagger}}ab^{{\dagger}2}bc$
$\displaystyle\qquad+270a^{{\dagger}}ab^{{\dagger}3}b^{2}c-416a^{{\dagger}}a^{3}b^{{\dagger}}b^{2}c^{{\dagger}}-128a^{{\dagger}}a^{3}bc^{{\dagger}}+416a^{{\dagger}}a^{3}bc^{{\dagger}2}c+160a^{{\dagger}2}b^{{\dagger}3}c^{3}$
$\displaystyle\qquad-496a^{{\dagger}2}a^{2}b^{{\dagger}}c^{{\dagger}}c^{2}+144a^{{\dagger}2}a^{2}b^{{\dagger}}c+496a^{{\dagger}2}a^{2}b^{{\dagger}2}bc-32a^{{\dagger}2}a^{4}bc^{{\dagger}}+32a^{{\dagger}3}a^{3}b^{{\dagger}}c$
$\displaystyle\qquad-401b^{{\dagger}}c^{{\dagger}}c^{2}+264b^{{\dagger}}c^{{\dagger}2}c^{3}-b^{{\dagger}}c^{{\dagger}3}c^{4}+b^{{\dagger}}c-602b^{{\dagger}2}bc^{{\dagger}}c^{2}+135b^{{\dagger}2}bc^{{\dagger}2}c^{3}+7b^{{\dagger}2}bc$
$\displaystyle\qquad-135b^{{\dagger}3}b^{2}c^{{\dagger}}c^{2}+6b^{{\dagger}3}b^{2}c+b^{{\dagger}4}b^{3}c-338a^{2}b^{{\dagger}}b^{2}c^{{\dagger}}+700a^{2}b^{{\dagger}}b^{2}c^{{\dagger}2}c-138a^{2}b^{{\dagger}2}b^{3}c^{{\dagger}}$
$\displaystyle\qquad-70a^{2}bc^{{\dagger}}+910a^{2}bc^{{\dagger}2}c-138a^{2}bc^{{\dagger}3}c^{2}\bigg{)}$
$\displaystyle+\tfrac{\chi^{8}}{8!}\bigg{(}-4680a^{{\dagger}}b^{{\dagger}2}c^{{\dagger}}c^{3}+408a^{{\dagger}}b^{{\dagger}2}c^{{\dagger}2}c^{4}+1104a^{{\dagger}}b^{{\dagger}2}c^{2}-2064a^{{\dagger}}b^{{\dagger}3}bc^{{\dagger}}c^{3}+1752a^{{\dagger}}b^{{\dagger}3}bc^{2}$
$\displaystyle\qquad+408a^{{\dagger}}b^{{\dagger}4}b^{2}c^{2}+168a^{{\dagger}}a^{2}b^{{\dagger}}b-11616a^{{\dagger}}a^{2}b^{{\dagger}}bc^{{\dagger}}c+7272a^{{\dagger}}a^{2}b^{{\dagger}}bc^{{\dagger}2}c^{2}+936a^{{\dagger}}a^{2}b^{{\dagger}2}b^{2}$
$\displaystyle\qquad-7272a^{{\dagger}}a^{2}b^{{\dagger}2}b^{2}c^{{\dagger}}c+408a^{{\dagger}}a^{2}b^{{\dagger}3}b^{3}-840a^{{\dagger}}a^{2}c^{{\dagger}}c+4536a^{{\dagger}}a^{2}c^{{\dagger}2}c^{2}-408a^{{\dagger}}a^{2}c^{{\dagger}3}c^{3}$
$\displaystyle\qquad+896a^{{\dagger}}a^{4}b^{2}c^{{\dagger}2}-3216a^{{\dagger}2}ab^{{\dagger}2}c^{{\dagger}}c^{3}+3264a^{{\dagger}2}ab^{{\dagger}2}c^{2}+3216a^{{\dagger}2}ab^{{\dagger}3}bc^{2}+272a^{{\dagger}2}a^{3}b^{{\dagger}}b$
$\displaystyle\qquad-3872a^{{\dagger}2}a^{3}b^{{\dagger}}bc^{{\dagger}}c+912a^{{\dagger}2}a^{3}b^{{\dagger}2}b^{2}-496a^{{\dagger}2}a^{3}c^{{\dagger}}c+912a^{{\dagger}2}a^{3}c^{{\dagger}2}c^{2}+1088a^{{\dagger}3}a^{2}b^{{\dagger}2}c^{2}$
$\displaystyle\qquad+64a^{{\dagger}3}a^{4}b^{{\dagger}}b-64a^{{\dagger}3}a^{4}c^{{\dagger}}c+ab^{{\dagger}}b-3602ab^{{\dagger}}bc^{{\dagger}}c+10410ab^{{\dagger}}bc^{{\dagger}2}c^{2}-820ab^{{\dagger}}bc^{{\dagger}3}c^{3}$
$\displaystyle\qquad+7ab^{{\dagger}2}b^{2}-4134ab^{{\dagger}2}b^{2}c^{{\dagger}}c+3414ab^{{\dagger}2}b^{2}c^{{\dagger}2}c^{2}+6ab^{{\dagger}3}b^{3}-820ab^{{\dagger}3}b^{3}c^{{\dagger}}c+ab^{{\dagger}4}b^{4}$
$\displaystyle\qquad-239ac^{{\dagger}}c+3607ac^{{\dagger}2}c^{2}-810ac^{{\dagger}3}c^{3}+ac^{{\dagger}4}c^{4}+1392a^{3}b^{{\dagger}}b^{3}c^{{\dagger}2}+1792a^{3}b^{2}c^{{\dagger}2}-1392a^{3}b^{2}c^{{\dagger}3}c\bigg{)}$
$\displaystyle+\dots$

And the time-evolved operator for the idler mode $b_{o}$ can be found by
swapping $a$ and $b$ in the above formula. The formula for the pump mode is

$\displaystyle c_{o}$ $\displaystyle=c-i\chi
ab+\tfrac{\chi^{2}}{2}\quantity(-a^{{\dagger}}ac-b^{{\dagger}}bc-c)+\tfrac{i\chi^{3}}{3!}\quantity(2a^{{\dagger}}b^{{\dagger}}c^{2}+a^{{\dagger}}a^{2}b+ab^{{\dagger}}b^{2}+ab-2abc^{{\dagger}}c)$
$\displaystyle+\tfrac{\chi^{4}}{4!}\quantity(10a^{{\dagger}}ab^{{\dagger}}bc-4a^{{\dagger}}ac^{{\dagger}}c^{2}+3a^{{\dagger}}ac+a^{{\dagger}2}a^{2}c-4b^{{\dagger}}bc^{{\dagger}}c^{2}+3b^{{\dagger}}bc+b^{{\dagger}2}b^{2}c-4c^{{\dagger}}c^{2}-4a^{2}b^{2}c^{{\dagger}}+c)$
$\displaystyle+\tfrac{i\chi^{5}}{5!}\bigg{(}8a^{{\dagger}}b^{{\dagger}}c^{{\dagger}}c^{3}-20a^{{\dagger}}b^{{\dagger}}c^{2}-16a^{{\dagger}}b^{{\dagger}2}bc^{2}-14a^{{\dagger}}a^{2}b^{{\dagger}}b^{2}-3a^{{\dagger}}a^{2}b+28a^{{\dagger}}a^{2}bc^{{\dagger}}c-16a^{{\dagger}2}ab^{{\dagger}}c^{2}$
$\displaystyle\qquad-a^{{\dagger}2}a^{3}b-3ab^{{\dagger}}b^{2}+28ab^{{\dagger}}b^{2}c^{{\dagger}}c-ab^{{\dagger}2}b^{3}-ab+40abc^{{\dagger}}c-8abc^{{\dagger}2}c^{2}\bigg{)}$
$\displaystyle+\tfrac{\chi^{6}}{6!}\bigg{(}216a^{{\dagger}}ab^{{\dagger}}bc^{{\dagger}}c^{2}-148a^{{\dagger}}ab^{{\dagger}}bc-91a^{{\dagger}}ab^{{\dagger}2}b^{2}c+148a^{{\dagger}}ac^{{\dagger}}c^{2}-16a^{{\dagger}}ac^{{\dagger}2}c^{3}-7a^{{\dagger}}ac+44a^{{\dagger}}a^{3}b^{2}c^{{\dagger}}$
$\displaystyle\qquad-40a^{{\dagger}2}b^{{\dagger}2}c^{3}-91a^{{\dagger}2}a^{2}b^{{\dagger}}bc+44a^{{\dagger}2}a^{2}c^{{\dagger}}c^{2}-6a^{{\dagger}2}a^{2}c-a^{{\dagger}3}a^{3}c+148b^{{\dagger}}bc^{{\dagger}}c^{2}-16b^{{\dagger}}bc^{{\dagger}2}c^{3}-7b^{{\dagger}}bc$
$\displaystyle\qquad+44b^{{\dagger}2}b^{2}c^{{\dagger}}c^{2}-6b^{{\dagger}2}b^{2}c-b^{{\dagger}3}b^{3}c+60c^{{\dagger}}c^{2}-16c^{{\dagger}2}c^{3}+44a^{2}b^{{\dagger}}b^{3}c^{{\dagger}}+60a^{2}b^{2}c^{{\dagger}}-72a^{2}b^{2}c^{{\dagger}2}c-c\bigg{)}$
$\displaystyle+\tfrac{i\chi^{7}}{7!}\bigg{(}-704a^{{\dagger}}b^{{\dagger}}c^{{\dagger}}c^{3}+32a^{{\dagger}}b^{{\dagger}}c^{{\dagger}2}c^{4}+222a^{{\dagger}}b^{{\dagger}}c^{2}-416a^{{\dagger}}b^{{\dagger}2}bc^{{\dagger}}c^{3}+490a^{{\dagger}}b^{{\dagger}2}bc^{2}+138a^{{\dagger}}b^{{\dagger}3}b^{2}c^{2}$
$\displaystyle\qquad+208a^{{\dagger}}a^{2}b^{{\dagger}}b^{2}-1204a^{{\dagger}}a^{2}b^{{\dagger}}b^{2}c^{{\dagger}}c+135a^{{\dagger}}a^{2}b^{{\dagger}2}b^{3}+7a^{{\dagger}}a^{2}b-1022a^{{\dagger}}a^{2}bc^{{\dagger}}c+496a^{{\dagger}}a^{2}bc^{{\dagger}2}c^{2}$
$\displaystyle\qquad-416a^{{\dagger}2}ab^{{\dagger}}c^{{\dagger}}c^{3}+490a^{{\dagger}2}ab^{{\dagger}}c^{2}+700a^{{\dagger}2}ab^{{\dagger}2}bc^{2}+135a^{{\dagger}2}a^{3}b^{{\dagger}}b^{2}+6a^{{\dagger}2}a^{3}b-270a^{{\dagger}2}a^{3}bc^{{\dagger}}c$
$\displaystyle\qquad+138a^{{\dagger}3}a^{2}b^{{\dagger}}c^{2}+a^{{\dagger}3}a^{4}b+7ab^{{\dagger}}b^{2}-1022ab^{{\dagger}}b^{2}c^{{\dagger}}c+496ab^{{\dagger}}b^{2}c^{{\dagger}2}c^{2}+6ab^{{\dagger}2}b^{3}-270ab^{{\dagger}2}b^{3}c^{{\dagger}}c$
$\displaystyle\qquad+ab^{{\dagger}3}b^{4}+ab-522abc^{{\dagger}}c+848abc^{{\dagger}2}c^{2}-32abc^{{\dagger}3}c^{3}+160a^{3}b^{3}c^{{\dagger}2}\bigg{)}$
$\displaystyle+\tfrac{\chi^{8}}{8!}\bigg{(}-17472a^{{\dagger}}ab^{{\dagger}}bc^{{\dagger}}c^{2}+3872a^{{\dagger}}ab^{{\dagger}}bc^{{\dagger}2}c^{3}+1826a^{{\dagger}}ab^{{\dagger}}bc-7272a^{{\dagger}}ab^{{\dagger}2}b^{2}c^{{\dagger}}c^{2}+3246a^{{\dagger}}ab^{{\dagger}2}b^{2}c$
$\displaystyle\qquad+820a^{{\dagger}}ab^{{\dagger}3}b^{3}c-3768a^{{\dagger}}ac^{{\dagger}}c^{2}+3376a^{{\dagger}}ac^{{\dagger}2}c^{3}-64a^{{\dagger}}ac^{{\dagger}3}c^{4}+15a^{{\dagger}}ac-2064a^{{\dagger}}a^{3}b^{{\dagger}}b^{3}c^{{\dagger}}$
$\displaystyle\qquad-1512a^{{\dagger}}a^{3}b^{2}c^{{\dagger}}+3216a^{{\dagger}}a^{3}b^{2}c^{{\dagger}2}c-896a^{{\dagger}2}b^{{\dagger}2}c^{{\dagger}}c^{4}+2384a^{{\dagger}2}b^{{\dagger}2}c^{3}+1392a^{{\dagger}2}b^{{\dagger}3}bc^{3}$
$\displaystyle\qquad-7272a^{{\dagger}2}a^{2}b^{{\dagger}}bc^{{\dagger}}c^{2}+3246a^{{\dagger}2}a^{2}b^{{\dagger}}bc+3414a^{{\dagger}2}a^{2}b^{{\dagger}2}b^{2}c-2736a^{{\dagger}2}a^{2}c^{{\dagger}}c^{2}+912a^{{\dagger}2}a^{2}c^{{\dagger}2}c^{3}$
$\displaystyle\qquad+25a^{{\dagger}2}a^{2}c-408a^{{\dagger}2}a^{4}b^{2}c^{{\dagger}}+1392a^{{\dagger}3}ab^{{\dagger}2}c^{3}+820a^{{\dagger}3}a^{3}b^{{\dagger}}bc-408a^{{\dagger}3}a^{3}c^{{\dagger}}c^{2}+10a^{{\dagger}3}a^{3}c$
$\displaystyle\qquad+a^{{\dagger}4}a^{4}c-3768b^{{\dagger}}bc^{{\dagger}}c^{2}+3376b^{{\dagger}}bc^{{\dagger}2}c^{3}-64b^{{\dagger}}bc^{{\dagger}3}c^{4}+15b^{{\dagger}}bc-2736b^{{\dagger}2}b^{2}c^{{\dagger}}c^{2}+912b^{{\dagger}2}b^{2}c^{{\dagger}2}c^{3}$
$\displaystyle\qquad+25b^{{\dagger}2}b^{2}c-408b^{{\dagger}3}b^{3}c^{{\dagger}}c^{2}+10b^{{\dagger}3}b^{3}c+b^{{\dagger}4}b^{4}c-744c^{{\dagger}}c^{2}+1552c^{{\dagger}2}c^{3}-64c^{{\dagger}3}c^{4}-1512a^{2}b^{{\dagger}}b^{3}c^{{\dagger}}$
$\displaystyle\qquad+3216a^{2}b^{{\dagger}}b^{3}c^{{\dagger}2}c-408a^{2}b^{{\dagger}2}b^{4}c^{{\dagger}}-744a^{2}b^{2}c^{{\dagger}}+6384a^{2}b^{2}c^{{\dagger}2}c-1088a^{2}b^{2}c^{{\dagger}3}c^{2}+c\bigg{)}$
$\displaystyle+\dots$

###  I.2 Full Output Operators to $\mathcal{O}(\chi^{2})$

We perform a $c=\alpha+\delta c$ expansion, and retain only the terms of the
form $\alpha^{n}\chi^{n}$ (assumed to be of order 1) and
$\alpha^{n-1}\chi^{n}$ (of $\mathcal{O}(\chi)$). For example, with the pump
mode $c_{o}$, only a small number of terms could potentially contribute at
$\mathcal{O}(\chi)$:

$\displaystyle c_{o}$ $\displaystyle=c-i\chi
ab+\tfrac{\chi^{2}}{2}(-c-a^{\dagger}ac-b^{\dagger}bc)$
$\displaystyle+\tfrac{i\chi^{3}}{3!}(2a^{\dagger}b^{\dagger}c^{2}-2abc^{\dagger}c)+\tfrac{\chi^{4}}{4!}(-4c^{\dagger}c^{2}-4a^{\dagger}ac^{\dagger}c^{2}-4b^{\dagger}bc^{\dagger}c^{2})$
$\displaystyle+\tfrac{i\chi^{5}}{5!}(8a^{\dagger}b^{\dagger}c^{\dagger}c^{3}-8abc^{{\dagger}2}c^{2})+\tfrac{\chi^{6}}{6!}(-16c^{{\dagger}2}c^{3}-16a^{\dagger}ac^{{\dagger}2}c^{3}-16b^{\dagger}bc^{{\dagger}2}c^{3})$
$\displaystyle+\tfrac{i\chi^{7}}{7!}(32a^{\dagger}b^{\dagger}c^{{\dagger}2}c^{4}-32abc^{{\dagger}3}c^{3})+\tfrac{\chi^{8}}{8!}(-64c^{{\dagger}3}c^{4}-64a^{\dagger}ac^{{\dagger}3}c^{4}-64b^{\dagger}bc^{{\dagger}3}c^{4})+\dots$

Now substitute $c=\alpha+\delta c$ (and therefore $c^{\dagger}=\alpha+\delta
c^{\dagger})$, keeping only $\mathcal{O}(\chi)$ terms, we have

$\displaystyle c_{o}$ $\displaystyle=\ \alpha+\delta c-i\chi
ab+\tfrac{\chi^{2}}{2!}\quantity(-\alpha-\alpha a^{{\dagger}}a-\alpha
b^{{\dagger}}b)$
$\displaystyle+\tfrac{i\chi^{3}}{3!}\quantity(2\alpha^{2}a^{{\dagger}}b^{{\dagger}}-2\alpha^{2}ab)+\tfrac{\chi^{4}}{4!}\quantity(-4\alpha^{3}-4\alpha^{3}a^{{\dagger}}a-4\alpha^{3}b^{{\dagger}}b)$
$\displaystyle+\tfrac{i\chi^{5}}{5!}\quantity(8\alpha^{4}a^{{\dagger}}b^{{\dagger}}-8\alpha^{4}ab)+\tfrac{\chi^{6}}{6!}\quantity(-16\alpha^{5}-16\alpha^{5}a^{{\dagger}}a-16\alpha^{5}b^{{\dagger}}b)$
$\displaystyle+\tfrac{i\chi^{7}}{7!}\quantity(32\alpha^{6}a^{{\dagger}}b^{{\dagger}}-32\alpha^{6}ab)+\tfrac{\chi^{8}}{8!}\quantity(-64\alpha^{7}-64\alpha^{7}a^{{\dagger}}a-64\alpha^{7}b^{{\dagger}}b)+\dots$

Group the terms according to the operators they are multiplied to.

$\displaystyle\Rightarrow c_{o}$
$\displaystyle=\alpha-(\tfrac{\alpha\chi^{2}}{2!}+\tfrac{4\alpha^{3}\chi^{4}}{4!}+\tfrac{16\alpha^{5}\chi^{6}}{6!}+\tfrac{64\alpha^{7}\chi^{8}}{8!}+\dots)$
$\displaystyle+\delta
c-(a^{\dagger}a+b^{\dagger}b)(\tfrac{\alpha\chi^{2}}{2!}+\tfrac{4\alpha^{3}\chi^{4}}{4!}+\tfrac{16\alpha^{5}\chi^{6}}{6!}+\tfrac{64\alpha^{7}\chi^{8}}{8!}+\dots)$
$\displaystyle+ia^{\dagger}b^{\dagger}(\tfrac{2\alpha^{2}\chi^{3}}{3!}+\tfrac{8\alpha^{4}\chi^{5}}{5!}+\tfrac{32\alpha^{6}\chi^{7}}{7!}+\dots)$
$\displaystyle-
iab(\chi+\tfrac{2\alpha^{2}\chi^{3}}{3!}+\tfrac{8\alpha^{4}\chi^{5}}{5!}+\tfrac{32\alpha^{6}\chi^{7}}{7!}+\dots).$

Note just as we decomposed $c=\alpha+\delta c$, we can decompose the
Heisenberg operator into an amplitude part and a noise part,
$c_{o}=\alpha_{o}+\delta c_{o}$, both parts are time-dependent. We see from
above that the first line of $c_{o}$ consists only of pure numbers and no
operators, and is therefore the amplitude $\alpha_{o}$; anything on the second
line and below are the noise part $\delta c_{o}$.

There are clear patterns to the first few terms of each infinite series
appeared above. Assuming the patterns persist indefinitely (checked to order
$\alpha^{14}\chi^{15}$), we may express each infinite series as a sum, and use
Mathematica to find the closed forms of these series:

$\displaystyle\tfrac{\alpha\chi^{2}}{2!}+\tfrac{4\alpha^{3}\chi^{4}}{4!}+\tfrac{16\alpha^{5}\chi^{6}}{6!}+\dots$
$\displaystyle=\sum_{n=1}^{\infty}\tfrac{\alpha^{2n-1}\chi^{2n}\cdot
4^{n-1}}{(2n)!}=\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime},$
$\displaystyle\tfrac{2\alpha^{2}\chi^{3}}{3!}+\tfrac{8\alpha^{4}\chi^{5}}{5!}+\tfrac{32\alpha^{6}\chi^{7}}{7!}+\dots$
$\displaystyle=\sum_{n=1}^{\infty}\tfrac{\alpha^{2n}\chi^{2n+1}\cdot 2\cdot
4^{n-1}}{(2n+1)!}=-\tfrac{\chi}{2}+\tfrac{\chi}{2\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime},$
$\displaystyle\chi+\tfrac{2\alpha^{2}\chi^{3}}{3!}+\tfrac{8\alpha^{4}\chi^{5}}{5!}+\tfrac{32\alpha^{6}\chi^{7}}{7!}+\dots$
$\displaystyle=\tfrac{\chi}{2}+\tfrac{\chi}{2\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime},$

where we again defined $\chi^{\prime}\equiv\alpha\chi$. So $c_{o}$ is

$\displaystyle\begin{split}c_{o}&=\alpha_{o}+\delta c_{o}\\\
&=\alpha-\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}+\delta
c-(a^{\dagger}a+b^{\dagger}b)\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}\\\
&-ia^{\dagger}b^{\dagger}\tfrac{\chi}{2}(1-\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})-iab\tfrac{\chi}{2}(1+\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime}),\end{split}$

where we see that

$\displaystyle\alpha_{o}$
$\displaystyle=\alpha-\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime},$
$\displaystyle\text{and }\ \delta c_{o}$ $\displaystyle=\delta
c-(a^{\dagger}a+b^{\dagger}b)\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}-ia^{\dagger}b^{\dagger}\tfrac{\chi}{2}(1-\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})-iab\tfrac{\chi}{2}(1+\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime}).$

For a fully self-consistent model capable of calculating non-trivial
expectation values, we need to include the $\mathcal{O}(\chi^{2})$ terms in
the output modes $a_{o},b_{o},c_{o}$ as well. That is, after the
$c=\alpha+\delta c$ expansion, on top of the $\alpha^{n}\chi^{n}$ and
$\alpha^{n-1}\chi^{n}$ terms, we now also retain terms of the form
$\alpha^{n-2}\chi^{n}$. This introduces many additional infinite series. For
illustrative purposes, let’s focus on a couple of them. In $c_{o}$, consider
terms proportional to $a^{\dagger}ac^{{\dagger}n}c^{m}$, the terms that could
contribute are

$\displaystyle c_{o}$
$\displaystyle=\dots-\tfrac{\chi^{2}}{2}a^{\dagger}ac-\tfrac{\chi^{4}}{4!}4a^{\dagger}ac^{\dagger}c^{2}-\tfrac{\chi^{6}}{6!}16a^{\dagger}ac^{{\dagger}2}c^{3}-\tfrac{\chi^{8}}{8!}64a^{\dagger}ac^{{\dagger}3}c^{4}+\dots,$

now perform the $c=\alpha+\delta c$ expansion, we have

$\displaystyle c_{o}$
$\displaystyle=\dots-a^{\dagger}a(\tfrac{\chi^{2}}{2}(\alpha+\delta
c)+\tfrac{\chi^{4}}{4!}4(\alpha+\delta c^{\dagger})(\alpha+\delta
c)^{2}+\tfrac{\chi^{6}}{6!}16(\alpha+\delta c^{\dagger})^{2}(\alpha+\delta
c)^{3}+\tfrac{\chi^{8}}{8!}64(\alpha+\delta c^{\dagger})^{3}(\alpha+\delta
c)^{4}+\dots)$
$\displaystyle=\dots-a^{\dagger}a(\tfrac{\chi^{2}}{2}\alpha+\tfrac{4\chi^{4}}{4!}\alpha^{3}+\tfrac{16\chi^{6}}{6!}\alpha^{5}+\tfrac{64\chi^{8}}{8!}\alpha^{7}+\dots)$
$\displaystyle\qquad\ \,-a^{\dagger}a(\tfrac{\chi^{2}}{2}\delta
c+\tfrac{8\chi^{4}}{4!}\alpha^{2}\delta
c+\tfrac{48\chi^{6}}{6!}\alpha^{4}\delta
c+\tfrac{256\chi^{8}}{8!}\alpha^{6}\delta c+\dots)$ $\displaystyle\qquad\
\,-a^{\dagger}a(\tfrac{4\chi^{4}}{4!}\alpha^{2}\delta
c^{\dagger}+\tfrac{32\chi^{6}}{6!}\alpha^{4}\delta
c^{\dagger}+\tfrac{192\chi^{8}}{8!}\alpha^{6}\delta c^{\dagger}+\dots),$

giving three infinite series. The first one we’ve already seen, it is of
$\mathcal{O}(\chi)$ and equals to
$-a^{\dagger}a\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}$. The last
two series are new, and can be put into closed form expressions as

$\displaystyle c_{o}$ $\displaystyle=\dots-a^{\dagger}a\delta
c\,\chi^{2}\sum_{n=0}^{\infty}2^{2n-1}(2n+2)\cdot\frac{\chi^{\prime
2n}}{(2n+2)!}-a^{\dagger}a\delta
c^{\dagger}\,\chi^{2}\sum_{n=1}^{\infty}2^{2n}n\frac{\chi^{\prime
2n}}{(2n+2)!}$ $\displaystyle=\dots-a^{\dagger}a\delta
c\,\chi^{2}\tfrac{1-\cosh
2\chi^{\prime}+2\sinh^{2}\chi^{\prime}+\chi^{\prime}\sinh
2\chi^{\prime}}{4\chi^{\prime 2}}-a^{\dagger}a\delta
c^{\dagger}\,\chi^{2}\tfrac{1-\cosh 2\chi^{\prime}+\chi^{\prime}\sinh
2\chi^{\prime}}{4\chi^{\prime 2}}.$

With some work, all second order terms can be grouped into series which can
then be expressed as closed form expressions like the ones above.

We now simply list the final results. We find that the signal mode $a_{o}$ to
order $\chi^{2}$ is

$\displaystyle a_{o}$ $\displaystyle=a\quantity[\cosh\chi^{\prime}+(\delta
c+\delta
c^{\dagger})\tfrac{\chi}{2}\sinh\chi^{\prime}]-ib^{\dagger}\big{[}\sinh\chi^{\prime}+\tfrac{\chi}{2}\cosh\chi^{\prime}(\delta
c+\delta c^{\dagger})+\tfrac{\chi}{2\chi^{\prime}}\sinh\chi^{\prime}(\delta
c-\delta c^{\dagger})\big{]}$
$\displaystyle+\chi^{2}\bigg{(}A_{a}a+A_{b^{\dagger}}b^{\dagger}+A_{a^{2}b}a^{2}b+A_{ab^{\dagger}b}ab^{\dagger}b+A_{a^{\dagger}ab^{\dagger}}a^{\dagger}ab^{\dagger}+A_{a^{\dagger}b^{{\dagger}2}}a^{\dagger}b^{{\dagger}2}+A_{a^{\dagger}a^{2}}a^{\dagger}a^{2}+A_{a\delta
c^{2}}a\delta c^{2}$ $\displaystyle+A_{a\delta c^{{\dagger}2}}a\delta
c^{{\dagger}2}+A_{a\delta c^{\dagger}\delta c}a\delta c^{\dagger}\delta
c+A_{b^{{\dagger}2}b}b^{{\dagger}2}b+A_{b^{\dagger}\delta
c^{2}}b^{\dagger}\delta c^{2}+A_{b^{\dagger}\delta c^{\dagger}\delta
c}b^{\dagger}\delta c^{\dagger}\delta c+A_{b^{\dagger}\delta
c^{{\dagger}2}}b^{\dagger}\delta c^{{\dagger}2}\bigg{)},$ (27)

where

$\displaystyle A_{a}$
$\displaystyle=-\sum_{n=2}^{\infty}\quantity(\tfrac{9^{n}-1}{8}-\tfrac{n(n+1)}{2})\cdot\tfrac{\chi^{\prime
2n-2}}{(2n)!}=-\tfrac{-\cosh\chi^{\prime}-\chi^{\prime
2}\cosh\chi^{\prime}+\cosh
3\chi^{\prime}-3\chi^{\prime}\sinh\chi^{\prime}}{8\chi^{\prime 2}},$
$\displaystyle A_{b^{\dagger}}$
$\displaystyle=i\sum_{n=1}^{\infty}Y(n)\tfrac{\chi^{\prime
2n-1}}{(2n+1)!}=i\cdot\tfrac{-5\chi^{\prime}\cosh\chi^{\prime}+2\sinh\chi^{\prime}-\chi^{\prime
2}\sinh\chi^{\prime}+\sinh 3\chi^{\prime}}{8\chi^{\prime 2}},$ $\displaystyle
A_{a^{2}b}$
$\displaystyle=-i\sum_{n=1}^{\infty}(a(n)+n)\cdot\tfrac{\chi^{\prime
2n-1}}{(2n+1)!}=-i\cdot\tfrac{4\chi^{\prime}\cosh\chi^{\prime}-7\sinh\chi^{\prime}+\sinh
3\chi^{\prime}}{16\chi^{\prime 2}},$ $\displaystyle A_{ab^{\dagger}b}$
$\displaystyle=-\sum_{n=1}^{\infty}\tfrac{9^{n}-1}{8}\cdot\tfrac{\chi^{\prime
2n-2}}{(2n)!}=\tfrac{\cosh\chi^{\prime}-\cosh 3\chi^{\prime}}{8\chi^{\prime
2}},$ $\displaystyle A_{a^{\dagger}ab^{\dagger}}$
$\displaystyle=i\sum_{n=1}^{\infty}(Y(n)+\tfrac{n(n+1)}{2})\cdot\tfrac{\chi^{\prime
2n-1}}{(2n+1)!}=i\cdot\tfrac{-4\chi^{\prime}\cosh\chi^{\prime}+\sinh\chi^{\prime}+\sinh
3\chi^{\prime}}{8\chi^{\prime 2}},$ $\displaystyle
A_{a^{\dagger}b^{{\dagger}2}}$
$\displaystyle=-A_{a^{\dagger}a^{2}}=\sum_{n=2}^{\infty}(\tfrac{9^{n}-1}{16}-\tfrac{n}{2})\cdot\tfrac{\chi^{\prime
2n-2}}{(2n)!}=\tfrac{-\cosh\chi^{\prime}+\cosh
3\chi^{\prime}-4\chi^{\prime}\sinh\chi^{\prime}}{16\chi^{\prime 2}},$
$\displaystyle A_{a\delta c^{2}}$ $\displaystyle=A_{a\delta
c^{{\dagger}2}}=\sum_{n=1}^{\infty}\tfrac{n(n+1)}{2}\cdot\tfrac{\chi^{\prime
2n}}{(2n+2)!}=\tfrac{\chi^{\prime}\cosh\chi^{\prime}-\sinh\chi^{\prime}}{8\chi^{\prime}},$
$\displaystyle A_{a\delta c^{\dagger}\delta c}$
$\displaystyle=\sum_{n=0}^{\infty}(n+1)^{2}\tfrac{\chi^{\prime
2n}}{(2n+2)!}=\tfrac{\chi^{\prime}\cosh\chi^{\prime}+\sinh\chi^{\prime}}{4\chi^{\prime}},$
$\displaystyle A_{b^{{\dagger}2}b}$
$\displaystyle=i\sum_{n=1}^{\infty}a(n)\cdot\tfrac{\chi^{\prime
2n-1}}{(2n+1)!}=i\cdot\tfrac{-4\chi^{\prime}\cosh\chi^{\prime}+\sinh\chi^{\prime}+\sinh
3\chi^{\prime}}{16\chi^{\prime 2}},$ $\displaystyle A_{b^{\dagger}\delta
c^{2}}$ $\displaystyle=\frac{1}{2}A_{b^{\dagger}\delta c^{\dagger}\delta
c}=-i\sum_{n=1}^{\infty}\tfrac{n(n+1)}{2}\cdot\tfrac{\chi^{\prime
2n-1}}{(2n+1)!}=-i\cdot\tfrac{\chi^{\prime}\cosh\chi^{\prime}-\sinh\chi^{\prime}+\chi^{\prime
2}\sinh\chi^{\prime}}{8\chi^{\prime}2},$ $\displaystyle A_{b^{\dagger}\delta
c^{{\dagger}2}}$
$\displaystyle=-i\sum_{n=1}^{\infty}(\tfrac{n(n+1)}{2})\cdot\tfrac{\chi^{\prime
2n+1}}{(2n+3)!}=-i\cdot\tfrac{-3\chi^{\prime}\cosh\chi^{\prime}+3\sinh\chi^{\prime}+\chi^{\prime
2}\sinh\chi^{\prime}}{8\chi^{\prime 2}}.$

And for the pump mode:

$\displaystyle c_{o}$
$\displaystyle=\alpha-\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}+\chi^{3}C_{\alpha}$
$\displaystyle+\delta
c-(a^{\dagger}a+b^{\dagger}b)\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}-ia^{\dagger}b^{\dagger}\tfrac{\chi}{2}(1-\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})-iab\tfrac{\chi}{2}(1+\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})$
$\displaystyle+\chi^{2}\bigg{(}C_{\delta c}\delta c+C_{a^{\dagger}a\delta
c}a^{\dagger}a\delta c+C_{b^{\dagger}b\delta c}b^{\dagger}b\delta
c+C_{ab\delta c}ab\delta c+C_{ab\delta c^{\dagger}}ab\delta
c^{\dagger}+C_{a^{\dagger}b^{\dagger}\delta c}a^{\dagger}b^{\dagger}\delta c$
$\displaystyle+C_{\delta c^{\dagger}}\delta c^{\dagger}+C_{a^{\dagger}a\delta
c^{\dagger}}a^{\dagger}a\delta c^{\dagger}+C_{b^{\dagger}b\delta
c^{\dagger}}b^{\dagger}b\delta c^{\dagger}+C_{a^{\dagger}b^{\dagger}\delta
c^{\dagger}}a^{\dagger}b^{\dagger}\delta c^{\dagger}\bigg{)},$ (28)

where

$\displaystyle C_{\alpha}$
$\displaystyle=\sum_{n=1}^{\infty}Z(n)\cdot\tfrac{\chi^{\prime
2n-1}}{(2n+2)!}=\tfrac{-3-4\chi^{\prime 2}+(2-4\chi^{\prime 2})\cosh
2\chi^{\prime}+\cosh 4\chi^{\prime}-2\chi^{\prime}\sinh
2\chi^{\prime}}{32\chi^{\prime 3}}$ $\displaystyle C_{\delta c}$
$\displaystyle=C_{a^{\dagger}a\delta c}=C_{b^{\dagger}b\delta
c}=-\sum_{n=0}^{\infty}2^{2n-1}(2n+2)\cdot\tfrac{\chi^{\prime
2n}}{(2n+2)!}=-\tfrac{1-\cosh
2\chi^{\prime}+2\sinh^{2}\chi^{\prime}+\chi^{\prime}\sinh
2\chi^{\prime}}{4\chi^{\prime 2}},$ $\displaystyle C_{ab\delta c}$
$\displaystyle=C_{ab\delta
c^{\dagger}}=-i\sum_{n=1}^{\infty}2^{2n-2}(2n)\cdot\tfrac{\chi^{\prime
2n-1}}{(2n+1)!}=-i\tfrac{2\chi^{\prime}\cosh 2\chi^{\prime}-\sinh
2\chi^{\prime}}{8\chi^{\prime 2}},$ $\displaystyle
C_{a^{\dagger}b^{\dagger}\delta c}$
$\displaystyle=i\sum_{n=2}^{\infty}\tfrac{2^{2n-2}n}{2}\cdot\tfrac{\chi^{\prime
2n-3}}{(2n-1)!}=i\tfrac{-4\chi^{\prime}+2\chi^{\prime}\cosh
2\chi^{\prime}+\sinh 2\chi^{\prime}}{8\chi^{\prime 2}},$ $\displaystyle
C_{\delta c^{\dagger}}$ $\displaystyle=C_{a^{\dagger}a\delta
c^{\dagger}}=C_{b^{\dagger}b\delta
c^{\dagger}}=-\sum_{n=1}^{\infty}2^{2n}n\tfrac{\chi^{\prime
2n}}{(2n+2)!}=-\tfrac{1-\cosh 2\chi^{\prime}+\chi^{\prime}\sinh
2\chi^{\prime}}{4\chi^{\prime 2}},$ $\displaystyle
C_{a^{\dagger}b^{\dagger}\delta c^{\dagger}}$
$\displaystyle=i\sum_{n=1}^{\infty}2^{2n}(2n)\cdot\tfrac{\chi^{\prime
2n+1}}{(2n+3)!}=i\tfrac{4\chi^{\prime}+2\chi^{\prime}\cosh
2\chi^{\prime}-3\sinh 2\chi^{\prime}}{8\chi^{\prime 2}}.$

In the above, the expression $a(n),X(n),Y(n)$ and $Z(n)$ are given by

$\displaystyle a(n)$ $\displaystyle=\tfrac{3^{2n+1}-8n-3}{16}.$ $\displaystyle
X(n)$ $\displaystyle=54a(n-1)+25n-18-\tfrac{n(n-1)}{2},$ $\displaystyle Y(n)$
$\displaystyle=18a(n-1)+7n-6-\tfrac{n(n-1)}{2},$ $\displaystyle\text{and
}\quad Z(n)$
$\displaystyle=(2n+1)!\sum_{k=0}^{n}\tfrac{1}{(2k)!}\tfrac{1}{[2(n-k)+1]!}\quantity[X(k-1)+Y(n-k)-k(n-k)],$

All $n$ and $k$’s are integers. The first few numbers in each sequence are
listed in Table I

$n$ | $a(n)$ | $X(n)$ | $Y(n)$ | $Z(n)$
---|---|---|---|---
$-1$ | $1/3$ | 0 | $2/3$ | —
0 | 0 | 0 | 0 | 0
1 | 1 | 7 | 1 | 1
2 | 14 | 85 | 25 | 60
3 | 135 | 810 | 264 | 1552
4 | 1228 | 7366 | 2446 | 29632
5 | 11069 | 66409 | 22123 | 506112
6 | 99642 | 597843 | 199263 | 8289280
Table 1:

One can check that the operators are physical in the sense that the
commutation relations are satisfied to $\mathcal{O}(\chi^{2})$, namely:

$\displaystyle[a_{o},a_{o}^{\dagger}]=[b_{o},b_{o}^{\dagger}]=[c_{o},c_{o}^{\dagger}]$
$\displaystyle=1+\mathcal{O}(\chi^{3}),$ all other commutation relations
$\displaystyle=0+\mathcal{O}(\chi^{3}).$

The vast majority of the terms in the output modes do not contribute to the
expectation values at $\mathcal{O}(\chi^{2})$. Specifically, it turns out the
only second order term in $a_{o}$ that contributes to
$\expectationvalue*{a_{o}^{\dagger}a_{o}}$ is
$\chi^{2}A_{b^{\dagger}}b^{\dagger}$; the only second order terms that
contribute to $\expectationvalue*{a_{o}b_{o}}$ are
$\chi^{2}(A_{b^{\dagger}}b^{\dagger}+A_{a}a)$. Similarly, the only second
order term in $c_{o}$ that contributes to $\expectationvalue*{\delta
c_{o}\delta c_{o}}$ is $\chi^{2}C_{\delta c^{\dagger}}\delta c^{\dagger}$, and
the only second order or above term that contributes to
$\expectationvalue*{c_{o}^{\dagger}c_{o}}$ is $\chi^{3}C_{\alpha}$. All other
terms either annihilate $\bra{0}$ or $\ket{0}$ to give 0 contributions, or
they only contribute to $\mathcal{O}(\chi^{3})$ terms. So for the purpose of
calculating expectation values, we may simply take

$\displaystyle a_{o}$ $\displaystyle=a\quantity[\cosh\chi^{\prime}+(\delta
c+\delta c^{\dagger})\tfrac{\chi}{2}\sinh\chi^{\prime}]$ $\displaystyle-
ib^{\dagger}\big{[}\sinh\chi^{\prime}+\tfrac{\chi}{2}\cosh\chi^{\prime}(\delta
c+\delta c^{\dagger})$
$\displaystyle\qquad+\tfrac{\chi}{2\chi^{\prime}}\sinh\chi^{\prime}(\delta
c-\delta c^{\dagger})\big{]}$
$\displaystyle\qquad+\chi^{2}(A_{a}a+A_{b^{\dagger}}b^{\dagger}),$ (29)
$\displaystyle b_{o}$ $\displaystyle=b\quantity[\cosh\chi^{\prime}+(\delta
c+\delta c^{\dagger})\tfrac{\chi}{2}\sinh\chi^{\prime}]$ $\displaystyle-
ia^{\dagger}\big{[}\sinh\chi^{\prime}+\tfrac{\chi}{2}\cosh\chi^{\prime}(\delta
c+\delta c^{\dagger})$
$\displaystyle\qquad+\tfrac{\chi}{2\chi^{\prime}}\sinh\chi^{\prime}(\delta
c-\delta c^{\dagger})\big{]}$
$\displaystyle\qquad+\chi^{2}(A_{a}b+A_{b^{\dagger}}a^{\dagger}),$ (30)
$\displaystyle c_{o}$ $\displaystyle=\alpha_{o}+\delta c_{o}$
$\displaystyle=\alpha-\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}+\chi^{3}C_{\alpha}$
$\displaystyle\quad+\delta
c-(a^{\dagger}a+b^{\dagger}b)\tfrac{\chi}{2\chi^{\prime}}\sinh^{2}\chi^{\prime}$
$\displaystyle\quad-
ia^{\dagger}b^{\dagger}\tfrac{\chi}{2}(1-\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})$
$\displaystyle\quad-
iab\tfrac{\chi}{2}(1+\tfrac{1}{\chi^{\prime}}\sinh\chi^{\prime}\cosh\chi^{\prime})$
$\displaystyle\quad+\chi^{2}C_{\delta c^{\dagger}}\delta c^{\dagger},$ (31)

These are ‘effective operators’ in the sense that they give the correct
results for
$\expectationvalue*{a_{o}^{\dagger}a_{o}},\expectationvalue*{a_{o}b_{o}},\expectationvalue*{\delta
c_{o}^{\dagger}\delta c_{o}},\expectationvalue*{\delta c_{o}\delta c_{o}}$ and
$\alpha_{o}^{2}$, and therefore all variances calculated from these operators
are correct. As such, one can also check that energy is conserved,

$\displaystyle\alpha^{2}=\expectationvalue*{c_{o}^{\dagger}c_{o}}+\expectationvalue*{a_{o}^{\dagger}a_{o}}=\alpha_{o}^{2}+\expectationvalue*{\delta
c_{o}^{\dagger}\delta c_{o}}+\expectationvalue*{a_{o}^{\dagger}a_{o}}.$

Although not obvious, it turns out the three effective operators above also
give the correct formulae for $\expectationvalue*{ab\delta c}$ and
$\expectationvalue*{ab\delta c^{\dagger}}$ up to $\mathcal{O}(\chi^{2})$.

The only drawback of using the effective operators is that the commutation
relations given by these operator are only correct to $\mathcal{O}(\chi)$, not
the desired $\mathcal{O}(\chi^{2})$. Therefore, we should not use them to
calculate non-normally ordered operator products, for example,
$\expectationvalue*{a_{o}a_{o}^{\dagger}}$ or $\expectationvalue*{\delta
c_{o}\delta c_{o}^{\dagger}}$. Instead we should first normal order them, such
that for example,
$\expectationvalue*{a_{o}a_{o}^{\dagger}}\to\expectationvalue*{1+a_{o}^{\dagger}a_{o}}$
and $\expectationvalue*{\delta c_{o}\delta
c_{o}^{\dagger}}\to\expectationvalue*{1+\delta c_{o}^{\dagger}\delta c_{o}}$,
before evaluating them.

The bottom line is, we may use equation (29) to (31) to calculate any second-
or third-order correlations to $\mathcal{O}(\chi^{2})$, provided the
correlations are normal ordered. If we want our theory to be fully self-
consistent without reordering and capable of predicting any correlations to
$\mathcal{O}(\chi^{2})$, we should use the full equations (27) and (28).