Stability and squeezing of the three-photon degenerate parametric down-conversion

Abstract

Nonlinear multiphoton conversion processes are theoretically and experimentally challenging systems in quantum optics and related areas. These phenomena bring remarkable quantum properties and can be of fundamental relevance in quantum information and computation applications. Based on recent experimental advances and the possible uses in generalized three-photon squeezing, we investigate degenerate three-photon parametric down-conversion generated inside a cavity fed by a classical pump. In addition, we consider that the process is stimulated by a coherent driving mode and subjected to dissipation by spontaneous emission of the cavity. A treatment using phase-space methods and stochastic differential equations is applied and the stability conditions of its steady states are found. We find that the system does not exhibit multistability, instead it has a limited parameter region leading to a single stable steady branch, alongside multiple unstable branches. Focusing on the stable steady states and linearizing the dynamics around them, we calculate the spectral densities of the quadrature variables. We identify and characterize the squeezing of three-photon down-conversion. Moreover, since the quasi-probability differential equations of the higher-order multiphoton processes have higher-order derivatives, hindering their mapping into stochastic differential equations, we make use of a positive Wigner function approach to calculate two-dimensional spectral densities associated with three-time correlation functions, with the goal of studying non-Gaussian properties of the system.

I Introduction

Optical parametric conversion processes are a key resource for implementing quantum states and protocols in Quantum Information and Computation (Gardiner-Zoller 2004; Scully-Zubairy 1997). Besides their origins in Quantum Optics, they are also implemented in circuit QED, for example (cQED21 2021). Recently, the importance of third-order nonlinear processes for completing a full set of quantum logic gates for quantum computation in the context of continuous variables has been highlighted (GKP01 2001; Zheng21 2021; Hahn22 2022). However, mathematical difficulties of the third-order squeezing operators or higher, such as divergent matrix elements in the Fock basis states, were found many years ago (Fisher84 1984). Posterior developments avoided these issues, showing that for evolutions in finite time, the states of these systems can be calculated (Braunstein87 1987; Braunstein90 1990; Ashhab25 2025). Furthermore, it was demonstrated that finite and regular results are obtained for models in which the squeezing terms are fed by quantized pumping operator, that is, exchange Hamiltonians involving interactions between modes with third-order or higher powers (Hillery90 1990; Drobny92 1992; Tanas92 1992; Banaszek97 1997; Olsen02 2002). Moreover, phase-space distribution representations were obtained, which exhibit very suggestive symmetry structures, directly associated with the order of the squeezing operators (Braunstein87 1987; Braunstein90 1990; Tanas92 1992; Banaszek97 1997; Felbinger98 1998). Within the context of the dissipative cavity systems, some works have studied steady state properties, finding interesting bistability and non-Gaussian features (Felbinger98 1998; Bajer91 1991; Denys21 2021).

Despite many experimental difficulties in handling third-order nonlinearities in material media, three-photon nondegenerate down-conversion has been achieved, by stimulating two output modes of a bulk crystal in the traveling-wave regime (Douady04 2004; Bencheikh07 2007; Gravier08 2008). More recently, another version of the experiment has been performed for a single stimulated output mode (Bertrand25 2025). Meanwhile, theoretical studies have modeled the operating conditions and shown the quantum features of these processes (Dot12 2012; Borshchevskaya15 2015; Okoth19 2019; Bencheikh22 2022). Another important achievement was the realization of spontaneous three-photon parametric down-conversion in cavities of superconducting circuits (SandboChang20 2020). In that work, both degenerate and nondegenerate cases were performed. It is worth mentioning active research in engineered optical fibers, waveguides, and metamaterials for producing three-photon parametric down-conversion processes in the traveling-wave regime (Corona11a 2011; Corona11b 2011; Akbari16 2016; Cavanna16 2016; Cavanna20 2020; Banic22 2022; Bacaoco25 2025).

Considering the use for the diverse possibilities in quantum information processing and computing, the three-photon down-conversion processes can produce more than genuine tripartite entanglement, like the Greenberger–Horne–Zeilinger (GHZ) states. In the case of the degenerate down-conversion, the nonlinear unitary transformation depends on the third-order power of the creation and annihilation operators, which has already been proposed as a way to obtain a cubic-phase transformation, necessary for Gottesman-Kitaev-Preskill (GKP) codes in the continuous-variable fault-tolerant quantum computation (GKP01 2001; Zheng21 2021; Hahn22 2022). Given these motivations, we have studied the three-photon degenerate parametric down-conversion performed within a cavity, taking into account certain realistic conditions, such as the cavity dissipation and the stimulation of the generated output photons by a driving coherent mode. Based on the experimental accomplishments (Bencheikh07 2007; Bertrand25 2025; SandboChang20 2020), we have considered the pumping mode to be sufficiently intense, so that depletion effects can be neglected and it can be treated as a classical wave. We have focused on the stationary states of this system, obtaining finite quadrature averages. However, unlike other driven dissipative systems (Felbinger98 1998; Bajer91 1991; Denys21 2021), we do not find multistable steady states, but only a limited region of single stability, with transitions to other unstable regions.

Furthermore, we have studied the spectra of quantum fluctuations of the generated mode around the stable steady states. We have obtained the spectral densities linked to correlation functions of the output quadratures, revealing third-order squeezing and its frequency-dependent structure, unlike the common squeezing. Due to the non-Gaussian nature of the noise, we have also taken into account the third-order noise components, which are commonly truncated, because of the representation issues with the noise in terms of stochastic differential equations (Olsen02 2002; Drummond14 2014). However, we have used an alternative approach, namely the positive Wigner functions, developed by Drummond (Drummond14 2014; Drummond17 2017). To calculate relevant quantities influenced by these higher-order noises, we use three-time correlation functions and calculate the respective two-dimensional frequency spectra, as is already done in the magnetic resonance spectroscopy (Ernst-Bodenhausen-Wokaun 1997). In this way, we propose a theoretical tool to calculate the statistical properties of non-Gaussian systems.

II Model

Many variations of the three-photon down-conversion or triple photon generation of nonlinear optical processes have been studied in the literature. In this article, we focus on the original problem proposed in theoretical works (Fisher84 1984; Braunstein87 1987; Braunstein90 1990; Ashhab25 2025), in which a generalization of squeezing by degenerate k-photon parametric amplifiers was studied, but here we restrict ourselves to the three-photon case. Moreover, we add the dynamics of the driven dissipative systems to get closer to experimental achievements, where third-order down-conversion implementations have been performed in recent years (Douady04 2004; Bencheikh07 2007; Gravier08 2008; Bertrand25 2025; SandboChang20 2020). In the physical model considered here, the nonlinear medium is fed by a pumping laser mode within an optical cavity. The pumping mode is sufficiently intense such that its depletion in down-conversion process is negligible, allowing us to treat it as a classical pump with frequency $\omega_{p}$. The triple photons produced by the parametric down-conversion are degenerate with frequency $\omega_{0}$. A driving mode with complex amplitude $F$ and frequency $\omega_{d}$ is added for stimulating the triple photon mode. So the nonlinear process can be modeled by the Hamiltonian

$$H_{0}=\hbar\omega_{0}a^{\dagger}a+i\hbar(Fe^{-i\omega_{d}t}a^{\dagger}-F^{*}e^{i\omega_{d}t}a)+i\hbar\kappa(e^{-i\omega_{p}t}a^{\dagger 3}-e^{i\omega_{p}t}a^{3}).$$ (1)

where $a$ is the annihilation operator of the subharmonic generated mode a. The real constant parameter $\kappa$ characterizes the nonlinear process and is proportional to the pumping mode amplitude.

To handle a time-independent Hamiltonian, we must perform a transformation to the rotating frame. The frequencies of the pumping and driving modes are externally controllable parameters, so we can set them as $\omega_{p}=3\omega_{d}$. Thus the Hamiltonian (1) can be reduced to

$$H_{0}=\hbar\Delta a^{\dagger}a+i\hbar(Fa^{\dagger}-F^{*}a)+i\hbar\kappa(a^{\dagger 3}-a^{3}),$$ (2)

where we define the detuning $\Delta\equiv\omega_{0}-\omega_{d}=\omega_{0}-\omega_{p}/3$.

In addition to the cavity dynamics described by expression (2), the dissipative part is described by single-photon dissipator at zero temperature,

$$\mathcal{D}(a)=\frac{\gamma}{2}\left(2a\rho a^{\dagger}-a^{\dagger}a\rho-\rho a^{\dagger}a\right).$$ (3)

Thus the density matrix of the cavity modes, $\rho$, is ruled by the Lindbladian master equation (Gardiner-Zoller 2004)

$$\frac{d\rho}{dt}=-\frac{i}{\hbar}[H_{0},\rho]+\mathcal{D}(a).$$ (4)

In what follows, we analyze the stability of the system focusing on the nonlinear interactions, so we must study the intracavity dynamics described by $a$ and $a^{\dagger}$ operators and their combinations. On the other hand, the study of spectral properties and noise levels is more interesting for the modes that are actually detected, therefore a method for calculating the extracavity modes is necessary. For this task, we approached the issue using the input/output fields method of Gardiner and Collett (Gardiner-Zoller 2004; Collett84 1984; Gardiner85 1985). In this way, we consider an input mode in the vacuum state, described by annihilation operator $a^{\mathrm{(in)}}(t)$, the intracavity modes with the respective operator $a(t)$, whose dynamics is ruled by equation (4), and output mode, with respective operator $a^{\mathrm{(out)}}(t)$, which follows the condition

$$a^{\mathrm{(out)}}(t)=a^{\mathrm{(in)}}(t)+\sqrt{\gamma}a(t).$$ (5)

To focus in the nonlinear process, we model the system as a one-sided cavity, that is, the cavity has relevant loss through only one mirror, whereas the pumping and driving modes are injected through other cavity ports without substantial losses.

In the spectral study of the system, we apply the Fourier analysis to the respective stochastic variables of the intracavity mode and we model the input mode with $[a^{\mathrm{(in)}}(t),a^{\mathrm{(in)\dagger}}(t^{\prime})]=\delta(t-t^{\prime})$, so that in next Section we associate it with a Wiener process (or, for the sake of integration, a white noise). Moreover, the input mode is understood as the resultant bath modes incoming to cavity, so we assume it as being in the vacuum state. In this way, we can obtain the spectral densities for the output mode by Eq. (5).

III Equations in phase-space

By standard methods (Gardiner 2004), the master equation (4) can be mapped into a partial differential equation for the Wigner function $W$, known as Kramers-Moyal equation (KME), taking the form

$$\frac{\partial W(\alpha,\alpha^{*})}{\partial t} = \left\{\frac{\partial}{\partial\alpha}\left[\Big(\frac{\gamma}{2}+i\Delta\Big)\alpha-F-3\kappa\alpha^{*2}\right]+\frac{\partial}{\partial\alpha^{*}}\left[\Big(\frac{\gamma}{2}-i\Delta\Big)\alpha^{*}-F^{*}-3\kappa\alpha^{2}\right]\right.$$ (6)

$$+ \left.\frac{\gamma}{2}\frac{\partial^{2}}{\partial\alpha\partial\alpha^{*}}-\frac{\kappa}{4}\left(\frac{\partial^{3}}{\partial\alpha^{3}}+\frac{\partial^{3}}{\partial\alpha^{*3}}\right)\right\}W(\alpha,\alpha^{*})$$ (7)

where the ordinary complex variables $\alpha$ and $\alpha^{*}$ are analogous to the operators $a$ and $a^{\dagger}$ in the domain of the Wigner function. Third or higher-order equations, as equation (6), do not have a direct correspondence with stochastic differential equations, because, in general, their solutions are not always positive to any conditions, so they are not genuine probability densities (Olsen02 2002; Drummond14 2014; Pawula67 1967). A method developed by Drummond to avoid this issue is to extend the domain of the $W$ function, so that the complex variables $\alpha$ and $\alpha^{*}$ become independent complex variables, not complex conjugate to each other, as $\alpha\rightarrow\alpha$ and $\alpha^{*}\rightarrow\beta$. Thus, the dimension of the original domain doubles to become a two-dimensional complex domain, bringing the new $W$ function to be necessarily positive, with analytic moments equivalent to those of the original Wigner function (Drummond14 2014; Drummond17 2017). This procedure follows closely the well-known positive-P function formalism (Drummond80a 1980), for this reason the new $W$ is called positive-W function.

The correspondence between the Fokker-Planck equations and the stochastic differential equations, obtained by the relations of the Itô’s calculus, is well-known (Gardiner 2004). Following Drummond (Drummond17 2017), we associate the KME (6) to a stochastic differential equation, extending the domain to a double complex space and, in particular, transforming the third-order derivatives in Eq. (6) into third-order noises. Thus, the differential equations for respective stochastic variables $\alpha(t)$ and $\beta(t)$ arranged in vector form, with their variations $dx=(d\alpha(t),d\beta(t))^{T}$, are

$$dx=Ddt+R^{(2)}dZ+R^{(3)}d\Xi,$$ (8)

such that the second and third-order noises are $dZ=(dZ_{1},dZ_{2})^{T}$ and $d\Xi=(d\Xi_{1},d\Xi_{2})^{T}$, respectively. The drift vector $D$ is

$$D=\left(\begin{array}[]{c}(-\frac{\gamma}{2}-i\Delta)\alpha+F+3\kappa\beta^{2}\\ (-\frac{\gamma}{2}+i\Delta)\beta+F^{*}+3\kappa\alpha^{2}\end{array}\right)$$ (9)

and the noise coefficients are matrices $R^{(2)}=\sqrt{\frac{\gamma}{4}}\left(\begin{array}[]{cc}1&i\\ 1&-i\end{array}\right)$ and $R^{(3)}=\sqrt[3]{\frac{3\kappa}{8}}\left(\begin{array}[]{cc}1&0\\ 0&1\end{array}\right)$. The noises are characterized by their correlations:

$$\left\langle dZ_{i}\right\rangle=0,$$ (10)

$$\left\langle dZ_{i}dZ_{j}\right\rangle=\delta_{ij}dt,$$ (11)

to the second-order noises and

$$\left\langle d\Xi_{i}\right\rangle=\left\langle d\Xi_{i}d\Xi_{j}\right\rangle=0,$$ (12)

$$\left\langle d\Xi_{i}d\Xi_{j}d\Xi_{k}\right\rangle=\delta_{ij}\delta_{jk}dt,$$ (13)

to the third-order ones. The system proposed in this article brings an advantage in use of the Drummond’s method because, whereas the KME second-order derivatives are non-diagonal but can be treated in the usual way as in the Fokker-Planck equation, the third-order derivatives correspond to the diagonal noise case, whose correlations follow a simple and direct form as in Eqs. (10) and (11). We can also note that the third-order noises are evidently non-Gaussian contributions to the stochastic behavior of the system and they have origin only from the nonlinear three-photon generation.

To study the squeezing and the respective spectrum of the three-photon down-conversion, it is useful to perform a change of variables to quadrature phase amplitudes, namely,

$$q = \frac{1}{2}(\alpha+\beta),$$ (14)

$$p = \frac{1}{2i}(\alpha-\beta),$$

so that the respective stochastic quadrature variables obey the following equations

$$dq=\left(-\frac{\gamma}{2}q+\Delta p+\Re(F)+3\kappa(q^{2}-p^{2})\right)dt+\frac{\sqrt{\gamma}}{2}Z_{1}+\frac{\sqrt[3]{3\kappa}}{2}\Xi_{+}$$ (15)

and

$$dp=\left(-\frac{\gamma}{2}p-\Delta q+\Im(F)-6\kappa qp\right)dt+\frac{\sqrt{\gamma}}{2}Z_{2}-i\frac{\sqrt[3]{3\kappa}}{2}\Xi_{-}.$$ (16)

Note that the second-order noise term is diagonal for the new variables of quadratures, whereas the third-order one is not. Even so, we can redefine the third-order noises as

$$\Xi_{+}=\frac{\Xi_{1}+\Xi_{2}}{\sqrt[3]{2}}$$ (17)

and

$$\Xi_{-}=\frac{\Xi_{1}-\Xi_{2}}{\sqrt[3]{2}},$$ (18)

whose correlation functions can be found as

$$\left\langle d\Xi_{+}d\Xi_{+}d\Xi_{+}\right\rangle=dt$$ (19)

and

$$\left\langle d\Xi_{+}d\Xi_{-}d\Xi_{-}\right\rangle=dt$$ (20)

to the nonvanishing cases. The other null correlations are

$$\left\langle d\Xi_{-}d\Xi_{-}d\Xi_{-}\right\rangle=\left\langle d\Xi_{+}d\Xi_{+}d\Xi_{-}\right\rangle=0$$ (21)

Looking at Eqs. (19)–(21), we note the non-trivial combination of noises for null and non-null correlation functions. This is consequence of the non-diagonal third-order noises for the differential quadrature equations (15) and (16) (Drummond17 2017).

IV Steady state conditions

Deterministic semi-classical equations of the quadrature variables can be obtained by calculating the expectations of the stochastic equations and applying the mean-field approximation, that is, in the present calculation, averages of variable products are factored as products of averages. In the stationary state, the stochastic variations vanish in the mean, so that the equations reduce to

$$-\frac{\gamma}{2}\left\langle q\right\rangle_{ss}+\Delta\left\langle p\right\rangle_{ss}+\Re(F)+3\kappa\left(\left\langle q\right\rangle_{ss}^{2}-\left\langle p\right\rangle_{ss}^{2}\right)=0$$ (22)

and

$$-\frac{\gamma}{2}\left\langle p\right\rangle_{ss}-\Delta\left\langle q\right\rangle_{ss}+\Im(F)-6\kappa\left\langle q\right\rangle_{ss}\left\langle p\right\rangle_{ss}=0,$$ (23)

where the $ss$ subscript denotes the steady state solutions.

We can decouple Eqs. (22) and (23), to obtain two independent quartic algebraic equations for both mean quadratures. Their solutions are very cumbersome and will be omitted. With these results, we can plot bifurcation diagrams and visualize the stable and unstable steady state branches. In Fig. 1, we show the bifurcation diagrams of the mean quadratures $\left\langle q\right\rangle_{ss}$, $\left\langle p\right\rangle_{ss}$, and the semi-classical energy $\left\langle E\right\rangle_{ss}=\left\langle q\right\rangle_{ss}^{2}+\left\langle p\right\rangle_{ss}^{2}$ as functions of $\Re(F)$, for the case with $\Im(F)=0$. In Fig. 2, we show the same diagrams, but the quantities are as functions of $\Im(F)$ and setting $\Re(F)=0$.

Figure 1: Steady state average values of the quadratures $\langle q\rangle_{ss}$ and $\langle p\rangle_{ss}$, and of the respective semi-classical energy $\langle E\rangle_{ss}=\langle q\rangle_{ss}^{2}+\langle p\rangle_{ss}^{2}$, as functions of the real part of the driving mode amplitude, $\Re(F)$, for the four physical solutions. Blue continuous lines represent the unique stable steady branch. The three other dashed lines represent unstable steady branches. The parameters for the three plots are $\gamma=1$, $\Delta=20$, $\kappa=10$, and $\Im(F)=0$.

Figure 2: Steady state average values of the quadratures $\langle q\rangle_{ss}$ and $\langle p\rangle_{ss}$, and of the respective semi-classical energy $\langle E\rangle_{ss}=\langle q\rangle_{ss}^{2}+\langle p\rangle_{ss}^{2}$, as functions of the imaginary part of the driving mode amplitude, $\Im(F)$, for the four physical solutions. Blue continuous lines represent the unique stable steady branch. The three other dashed lines represent unstable steady branches. The parameters for the three plots are $\gamma=1$, $\Delta=20$, $\kappa=10$, and $\Re(F)=0$.

In both Figs. 1 and 2, the stable steady states are represented only by blue lines, all other dotted lines represent unstable steady states. It is practically impossible for these unstable states to be attainable, so that the possibly observable steady states are restricted to a bounded set of parameters. Another aspect notable is the difference and the asymmetry of the diagrams, depending on the real and imaginary values of the driving mode. Such system’s features can be better understood analyzing its stability conditions (Strogatz 2018).

The stability of the branches presented in Figs. 1 and 2 can be established by the Jacobian of the average of the stochastic equations (15) and (16),

$$J=\begin{pmatrix}\partial\langle\dot{q}\rangle/\partial\langle q\rangle&\partial\langle\dot{q}\rangle/\partial\langle p\rangle\\ \partial\langle\dot{p}\rangle/\partial\langle q\rangle&\partial\langle\dot{p}\rangle/\partial\langle p\rangle\end{pmatrix},$$ (24)

hence

$$J=\begin{pmatrix}-\frac{\gamma}{2}+6\kappa\langle q\rangle&\Delta-6\kappa\langle p\rangle\\ -\Delta-6\kappa\langle p\rangle&-\frac{\gamma}{2}-6\kappa\langle q\rangle\end{pmatrix}.$$ (25)

For the identification of the stable or unstable steady states, the Jacobian must obey the following inequalities

$$\mathrm{tr}J=-\gamma\leq 0$$ (26)

and

$$\det J=\left(\frac{\gamma}{2}\right)^{2}+\Delta^{2}-36\kappa^{2}\left(\langle q\rangle^{2}+\langle p\rangle^{2}\right)\geq 0,$$ (27)

calculated to the solutions found from Eqs. (22) and (23). The first condition is clearly always met, but the second condition comprises a limited region. In this way, the stable and unstable steady state branches can be established in Figs. 1 and 2.

Figure 3: Stability region boundary as function of the real and imaginary parts of the driving mode amplitude. The driving mode amplitude values for which the system exhibits stable steady states lie within the region bounded by the contours. In plot (a), we fixed the detuning $\Delta=1$ and the nonlinear coupling constant $\kappa=10$ and varied the dissipative factor $\gamma$. In plot (b), we fixed $\gamma=1$ and $\kappa=10$ and varied $\Delta$. In plot (c), we fixed $\gamma=1$ and $\Delta=1$ and varied $\kappa$.

The boundary of the stability condition is found with $\det(J_{b})=0.$ Rewriting the quadrature variables in terms of polar variables, $q_{b}=r_{b}\cos\theta$ and $p_{b}=r_{b}\sin\theta$, the stability boundary is found as

$$q_{b}^{2}+p_{b}^{2}=r_{b}^{2}=\frac{\left(\frac{\gamma}{2}\right)^{2}+\Delta^{2}}{(6\kappa)^{2}}.$$ (28)

Substituting in the Eqs. (22) and (23), we obtain the stability boundary for the real and imaginary parts of the driving mode amplitude,

$$\Re(F_{b})=r_{b}\left[\frac{\gamma}{2}\cos\theta-\Delta\sin\theta-3\kappa r_{b}\left(\cos^{2}\theta-\sin^{2}\theta\right)\right]$$ (29)

and

$$\Im(F_{b})=r_{b}\left(\frac{\gamma}{2}\sin\theta+\Delta\cos\theta+6\kappa r_{b}\cos\theta\sin\theta\right)$$ (30)

Expressions (29) and (30) can be plotted in a complex plane of the driving mode amplitude, so that the inner region of the boundary line represents the stability region, as we see in Fig. 3, varying the parameters $\gamma$, $\Delta$ and $\kappa$.

The triply folded aspect of the stability region reveals an asymmetric projection onto both components of the driving mode amplitude, $\Re(F)$ and $\Im(F)$, reflecting the asymmetries in Figs. 1 and 2.

To complement our stability analysis, let us investigate the existence of steady states regarding the photon number of the mode generated within the cavity. We follow previous analyses which utilized the Heisenberg equation for the isolated system (Hillery90 1990; Olsen02 2002). However, given the inclusion of the quantum mode dissipation, we must account for a modification to the Heisenberg equation. To this goal, we consider the input-output field method (Gardiner-Zoller 2004; Collett84 1984; Gardiner85 1985). Thus, any operator $\mathcal{O}$ associated with the generated mode must satisfy the equation

$$\frac{d\mathcal{O}}{dt}=-\frac{1}{\hbar}[\mathcal{O},H_{0}]-\left[[\mathcal{O},a^{\dagger}]\left(\frac{\gamma}{2}a+\sqrt{\gamma}a^{\mathrm{(in)}}\right)-\left(\frac{\gamma}{2}a^{\dagger}+\sqrt{\gamma}a^{\mathrm{(in)}\dagger}\right)[\mathcal{O},a]\right].$$ (31)

To analyze the photon number evolution, we consider the operator $\mathcal{O}\rightarrow a^{\dagger}a$ and calculate up to its second time derivative. We obtain terms that depend purely on the intracavity mode operators, or products of them with the input modes, or a pure input term involving $a^{\mathrm{(in)}\dagger}a^{\mathrm{(in)}}$. Assuming the input mode bath is in the vacuum state, and the intracavity operators are reasonably described as non-anticipating quantities, so that they are uncorrelated to the input operators, then all the terms with the input mode operators must cancel out.

Figure 4: Mean photon number to steady states without driving mode, $F=0$. Blue continuous line represents the unique stable steady state branch. Red dashed lines represent the unstable steady states.

Again, making the mean-field approximation, but here assuming that all operator products above the second power factorize, the equation for mean values is

$$\frac{d^{2}\mathcal{N}}{dt} = i\Delta\left[\mathcal{N},\frac{d\mathcal{N}}{dt}\right]-\frac{5\gamma}{2}\frac{d\mathcal{N}}{dt}+2|F|^{2}+\gamma\left(F^{*}\langle a\rangle+F\langle a^{\dagger}\rangle\right)+12\kappa\left(F\langle a^{2}\rangle+F^{*}\langle a^{\dagger 2}\rangle\right)+$$ (32)

$$+ 18\kappa^{2}(3\mathcal{N}^{2}+3\mathcal{N}+2)-\frac{3\gamma^{2}}{2}\mathcal{N}.$$

where we define $\mathcal{N}=\langle a^{\dagger}a\rangle$. Of course this equation is incomplete, because it involves independent terms as $\langle a\rangle$ and $\langle a^{2}\rangle$ and any attempt to add other equations would result an infinite equation system. However, for the sake of comparison let us discard the driving mode, that is $F=0$, since it was not considered in the previous articles. In the steady state, all time derivatives vanish, yielding the algebraic equation for the respective mean photon number

$$18\kappa^{2}(3\mathcal{N}_{ss}^{2}+3\mathcal{N}_{ss}+2)-\frac{3\gamma^{2}}{2}\mathcal{N}_{ss}=0.$$ (33)

We have plotted a bifurcation diagram of Eq. (33), as it is shown in Fig. 4. We find that, unlike in previous studies (Hillery90 1990; Olsen02 2002), the average photon number converges asymptotically to a branch of stable steady states, as shown by the solid blue line in Fig. 4. The difference here lies in the contribution of the dissipative term in Eq. (33), which regularizes the system behavior and keeps the photon number finite. Interestingly, it can be observed that outside the stable steady state region, the average photon number diverges in finite time, revealing a remnant of the highly irregular behavior associated with three-photon down-conversion, a widely studied aspect (Fisher84 1984; Braunstein87 1987; Hillery90 1990; Olsen02 2002).

V Spectra in the stable states

An analysis of the spectral densities of the generated three-photon mode can be done for steady states using standard Fourier transform methods. The dynamical properties of the system in the neighborhood of its fixed points are sufficiently evaluated by linearizing around those points, thus it is pertinent to apply the Fourier analysis for finding stationary spectra near the stable fixed points (Gardiner 2004; Strogatz 2018).

Defining a vector of the quadrature variables centered in the stable steady states (sss),

$$Q=\begin{pmatrix}Q_{q}\\ Q_{p}\end{pmatrix}=\begin{pmatrix}q-\langle q\rangle_{sss}\\ p-\langle p\rangle_{sss}\end{pmatrix},$$ (34)

we can obtain the linearized stochastic differential equation calculating the Jacobian (25) to the stable steady branch, $A=J|_{sss}$. In the following analysis, it is suitable to use the Langevin equation formalism, since it is more natural to the calculations of the Fourier analysis. Hence, we get

$$\frac{dQ}{dt}=AQ+B\zeta+C\xi,$$ (35)

where $B$ and $C$ are matrix coefficients of the noise terms, given by

$$B=\sqrt{\frac{\gamma}{4}}\begin{pmatrix}1&0\\ 0&1\end{pmatrix}$$ (36)

and

$$C=\sqrt[3]{\frac{3\kappa}{8}}\begin{pmatrix}1&0\\ 0&-i\end{pmatrix}.$$ (37)

The noises in the Langevin equation (35) are obtained by matching the mean value calculations of stochastic variables, or using Ito´s integration with Wiener process or using the ordinary calculus with white noise, that is, we follow the prescription done by Gardiner in reference (Gardiner 2004). Roughly speaking, we take $dZ_{i}\rightarrow\zeta_{i}dt$ and $d\Xi_{i}\rightarrow\xi_{i}dt$. In this way, the noise vectors are

$$\zeta=\begin{pmatrix}\zeta_{1}\\ \zeta_{2}\end{pmatrix}\quad,\quad\xi=\begin{pmatrix}\xi_{+}\\ \xi_{-}\end{pmatrix},$$ (38)

characterized by the following mean values, $\langle\zeta_{i}(t)\rangle=\langle\xi_{i}(t)\rangle=0$. The second-order noises obey the correlations

$$\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime}).$$ (39)

To the third-order noise, we must obtain (see Section 7.7 in (Gardiner 2004) and (Gardiner77 1977))

$$\langle\xi_{+}(t)\xi_{+}(t^{\prime})\xi_{+}(t^{\prime\prime})\rangle=\langle\xi_{+}(t)\xi_{-}(t^{\prime})\xi_{-}(t^{\prime\prime})\rangle=\delta(t-t^{\prime})\delta(t-t^{\prime\prime})$$ (40)

and

$$\langle\xi_{+}(t)\xi_{+}(t^{\prime})\xi_{-}(t^{\prime\prime})\rangle=\langle\xi_{-}(t)\xi_{-}(t^{\prime})\xi_{-}(t^{\prime\prime})\rangle=0.$$ (41)

The Fourier transforms of the centered quadrature variables are

$$\tilde{Q}_{i}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}Q_{i}(t)e^{-i\omega t}dt,$$ (42)

and the respective inverse transforms are

$$Q_{i}(t)=\int_{-\infty}^{+\infty}\tilde{Q_{i}}(\omega)e^{i\omega t}d\omega.$$ (43)

Moreover, the relevant spectral densities are not concerning the intracavity variables, but are those related with the output mode. Based on Eq. (5), we have in terms of quadrature variables

$$Q^{\mathrm{(out)}}(t)=Q^{\mathrm{(in)}}(t)+\sqrt{\gamma}Q(t),$$ (44)

where the input variable vector is $Q^{\mathrm{(in)}}(t)=(q^{\mathrm{(in)}}(t),p^{\mathrm{(in)}}(t))^{T}=(\zeta_{1}(t),\zeta_{2}(t))^{T}$, that is, the input mode is just the second-order noise injected into the cavity. Thus the Fourier transform of Eq. (44) is

$$\tilde{Q}^{\mathrm{(out)}}(\omega)=\tilde{Q}^{\mathrm{(in)}}(\omega)+\sqrt{\gamma}\tilde{Q}(\omega),$$ (45)

so that the input mode is considered to be in the vacuum state, as presented in Section II. Hence, $\tilde{Q}_{i}^{\mathrm{(in)}}(\omega)=\tilde{\zeta}_{i}(\omega)$ with $\left\langle\tilde{\zeta}_{i}(\omega)\tilde{\zeta}_{j}^{*}(\omega^{\prime})\right\rangle=\delta_{ij}\delta(\omega-\omega^{\prime})$ (Gardiner-Zoller 2004; Gardiner85 1985). Thus the two-time correlation functions related to the output mode are

$$\phi_{ij}(\tau)=\left\langle Q_{i}^{\mathrm{(out)}}(t)Q_{j}^{\mathrm{(out)}}(t+\tau)\right\rangle,$$ (46)

and their Fourier transforms define the respective spectral densities (Gardiner 2004),

$$S_{ij}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\phi_{ij}(\tau)e^{-i\omega\tau}d\tau.$$ (47)

The spectral densities are also related to the correlation functions of the transformed quantities

$$\left\langle\tilde{Q_{i}}^{\mathrm{(out)}}(\omega)\tilde{Q_{j}}^{\mathrm{(out)}*}(\omega^{\prime})\right\rangle=S_{ij}(\omega)\delta(\omega-\omega^{\prime}),$$ (48)

of which we note that $S_{ij}(\omega)=S_{ji}(\omega)^{*}$. For notation simplicity, we use $S_{ii}(\omega)\equiv S_{i}(\omega)$.

Figure 5: Noise spectra of the quadratures $Q_{q}^{\mathrm{(out)}}$ and $Q_{p}^{\mathrm{(out)}}$ of the output modes from the cavity. Respectively, the blue continuous curve is the spectrum $S_{q}(\omega)$ and the red dashed curve is the spectrum $S_{p}(\omega)$. The parameters for the plots are $\gamma=1$, $\Delta=20$, $\kappa=10$, and $F=i3.9$. SQL is set to 1.

Figure 6: Cross noise spectrum of the quadrature product $Q_{q}^{\mathrm{(out)}}Q_{p}^{\mathrm{(out)}}$ of the output modes from the cavity. The spectrum $S_{qp}(\omega)$ is a complex function, so the blue continuous curve is its real part and the red dashed curve is its imaginary part. The parameters for the plots are $\gamma=1$, $\Delta=20$, $\kappa=10$, and $F=i3.9$.

The Fourier transforms of the intracavity centered quadrature variables are found as

$$\tilde{Q_{q}}(\omega)=\frac{f_{-}N_{1+}+g_{-}N_{2-}}{g_{+}g_{-}-f_{-}f_{+}}$$ (49)

and

$$\tilde{Q_{p}}(\omega)=\frac{f_{+}N_{2-}+g_{+}N_{1+}}{g_{+}g_{-}-f_{-}f_{+}},$$ (50)

where we use terms defined by

$$f_{\pm}(\omega)=-\left(\frac{\gamma}{2}+i\omega\right)\pm 6\kappa\langle q\rangle_{sss}$$ (51)

and

$$g_{\pm}=6\kappa\langle p\rangle_{sss}\pm\Delta,$$ (52)

and the noise functions are written as

$$N_{1+}=\sqrt{\gamma/4}\tilde{\zeta_{1}}(\omega)+\sqrt[3]{3\kappa/8}\tilde{\xi_{+}}(\omega)$$ (53)

and

$$N_{2-}=\sqrt{\gamma/4}\tilde{\zeta_{2}}(\omega)-i\sqrt[3]{3\kappa/8}\tilde{\xi_{-}}(\omega).$$ (54)

With all these pieces, we can finally calculate the spectral densities, including the spectra $S_{q}(\omega)$ and $S_{p}(\omega)$ related to the quadratures $Q_{q}^{\mathrm{(out)}}$ and $Q_{p}^{\mathrm{(out)}}$. Both spectra are real functions and are presented in Fig. (5). We observe the squeezing phenomenon depending on the frequency relative to the cavity resonance frequency. Around zero detuning, the quadrature $Q_{p}^{\mathrm{(out)}}$ exhibits noise squeezing, with its spectrum dropping below the vacuum level (the standard quantum limit - SQL), consequently the quadrature $Q_{q}^{\mathrm{(out)}}$ receives excess noise, remaining above the SQL. However, this scenario is reversed at higher frequency. After both quadratures exhibit noise peaks, it is $Q_{q}^{\mathrm{(out)}}$ that has squeezed noise, and $Q_{p}^{\mathrm{(out)}}$ noise raises. This intricate feature is a manifestation of a new type of squeezing, which is a characteristic trait of the third-order nonlinearity. For consistency, it is possible to show that the product of $S_{q}(\omega)$ and $S_{p}(\omega)$ obeys the Heisenberg inequality, that is, the product of both spectra must remain above the minimum value set by the vacuum state at all frequencies.

The cross spectral density $S_{qp}(\omega)$ is a complex function, so it must be presented in two plots, one to real part and the other to imaginary part, as shown in Fig. (6). The interpretation of $S_{qp}(\omega)$ is less clear. This function is related to covariance between the quadratures. Its value is associated with the angular position of the quasiprobability density in phase-space, reflecting the frequency’s dependence on the phase shift of the mode’s oscillatory components.

We note in both Figs. (5) and (6) resonant noise peaks at symmetric frequencies. These peaks are results of the denominators of expressions (49) and (50), which are complicated functions of the dissipation factor $\gamma$, detuning $\Delta$, and the nonlinearity constant $\kappa$. Another feature of the second-order spectra is the property $S_{ij}(-\omega)=S_{ij}(\omega)^{*}$, obtained for typical real processes (Gardiner-Zoller 2004; Gardiner 2004). This is not the case to the third-order spectra, presented in what follows. At this point, we should highlight that the degenerate three-photon down-conversion with classical pump yields a stochastic process in which the dissipative noise and the third-order nonlinear noise are explicitly separated in Eqs. (8) and (35). Hence the second and third-order statistical properties are very distinct and provide a clear scenario to test experimentally the nonlinear optical process and our theoretical method.

Since the system is fundamentally non-Gaussian, as is explicit from the third-order noise in Eq. (35), we must look for methods beyond second-order mean values. Inspired by the two-dimensional spectroscopy in Nuclear Magnetic Resonance (NMR) (Ernst-Bodenhausen-Wokaun 1997), we propose the introduction of multi-time correlation functions, in particular third-order three-time correlation functions, for studying non-Gaussian processes,

$$\phi_{ijk}(\tau_{1},\tau_{2})=\left\langle Q_{i}^{\mathrm{(out)}}(t+\tau_{1})Q_{j}^{\mathrm{(out)}}(t+\tau_{2})Q_{k}^{\mathrm{(out)}}(t)\right\rangle,$$ (55)

whose Fourier transforms define the respective two-dimensional (2D) spectral densities,

$$S_{ijk}(\omega_{1},\omega_{2})=\frac{1}{(2\pi)^{2}}\int_{-\infty}^{+\infty}d\tau_{1}e^{-i\omega_{1}\tau_{1}}\int_{-\infty}^{+\infty}d\tau_{2}e^{-i\omega_{2}\tau_{2}}\phi_{ijk}(\tau_{1},\tau_{2}).$$ (56)

In NMR, the 2D spectra are obtained as responses to pulse sequences that generally depend on two time parameters. Here, they indeed reflect the existence of non-null higher-order correlation functions in non-Gaussian processes, as well as higher-order cumulants must be non-null in such cases. Similar to the method presented here, a previous study also utilized multi-time correlations to investigate the spectrum generated by three-photon parametric down-conversion, but focused on the entanglement properties of the generated photons (Chekhova05 2020).

The 2D spectral densities, $S_{ijk}(\omega_{1},\omega_{2})$, are related with the third-order correlation function of the Fourier-transformed quantities by

$$\left\langle\tilde{Q_{i}}^{\mathrm{(out)}}(\omega_{1})\tilde{Q_{j}}^{\mathrm{(out)}}(\omega_{2})\tilde{Q_{k}}^{\mathrm{(out)}}(-\omega)\right\rangle=S_{ijk}(\omega_{1},\omega_{2})\delta(\omega_{1}+\omega_{2}-\omega)$$ (57)

which provides a way to calculate the respective spectra. This approach is also applied to the Fourier transforms of the third-order noise correlations (40) and (41). The nonvanishing correlations are transformed in

$$\left\langle\tilde{\xi}_{+}(\omega)\tilde{\xi}_{+}(\omega^{\prime})\tilde{\xi}_{+}(\omega^{\prime\prime})\right\rangle=\left\langle\tilde{\xi}_{+}(\omega)\tilde{\xi}_{-}(\omega^{\prime})\tilde{\xi}_{-}(\omega^{\prime\prime})\right\rangle=\delta(\omega+\omega^{\prime}+\omega^{\prime\prime}),$$ (58)

and the null correlation cases are trivial. Bringing together the transformed quadratures (49) and (50) and Eqs. (57) and (58), we find the third-order spectra $S_{ijk}(\omega_{1},\omega_{2})$ and plot the respective graphs, presented in Figs. 7 and 8.

Figure 7: 2D Noise spectra of the output quadratures $Q_{q}^{\mathrm{(out)}}$ and $Q_{p}^{\mathrm{(out)}}$. The spectra $S_{q^{3}}\equiv S_{qqq}$ and $S_{p^{3}}\equiv S_{ppp}$ are complex functions, so the real and imaginary parts are shown separately. The parameters for the plots are $\gamma=1$, $\Delta=20$, $\kappa=10$, and $F=i3.9$.

Figure 8: 2D Cross noise spectra of the combinations of three output quadrature products for $Q_{q}^{\mathrm{(out)}}$ and $Q_{p}^{\mathrm{(out)}}$. The spectra $S_{q^{2}p}\equiv S_{qqp}$ and $S_{p^{2}q}\equiv S_{ppq}$ are complex functions, so the real and imaginary parts are shown separately. The parameters for the plots are $\gamma=1$, $\Delta=20$, $\kappa=10$, and $F=i3.9$.

As in the magnetic resonance (Ernst-Bodenhausen-Wokaun 1997), the functions $S_{ijk}(\omega_{1},\omega_{2})$ are complex and generally manifest mixed absorptive and dispersive aspects. In the present article, we find four relevant 2D spectral densities: $S_{qqq}(\omega_{1},\omega_{2})$, $S_{ppp}(\omega_{1},\omega_{2})$, $S_{qqp}(\omega_{1},\omega_{2})$, $S_{ppq}(\omega_{1},\omega_{2})$. Other cases are only rotated and rescaled versions, which do not add information.

Third-order spectral densities are all complex functions, as they are Fourier transforms of three-time correlation functions, which are themselves complex functions, given that there is no direct way to express them as symmetric products of conjugate variables, unlike the case for two-time correlation functions. The fact that three-time correlation functions and their respective 2D spectra are complex-valued functions is not a novelty, but an inherent feature of this spectroscopic analysis method (Ernst-Bodenhausen-Wokaun 1997). The repetitive structure observed in Figs. 7 and 8 suggests significant information redundancy in these plots, nevertheless we present all these compositions for the sake of completeness. As with 2D NMR, we can analyze the resonance structures based on their response, whether they are dispersive and/or absorptive, with the latter aspect understood as the dissipative nature of the cavity. It is observed that the real parts of the spectra consist of hybrid responses comprising both dissipative and dispersive components, in turn the imaginary parts are predominantly dispersive. This dispersive predominance is characteristic of the nonlinear process, which inherently involves a phase shift in the system behavior. Furthermore, the frequency values $\omega_{1}$ and $\omega_{2}$ at which noise peaks appear also encode information regarding the nonlinear dynamics, as they are determined by the denominators of Eqs. (49) and (50), following the features of the second-order spectral densities. The triple products of these denominators in Eq. (57) and their off-resonance noise peaks are the reason to the sixfold repetitive structure in the 2D spectra.

VI Discussion

Given the known issues regarding the non-existence of steady states in a pure generalized squeezing Hamiltonian (Fisher84 1984; Braunstein87 1987; Braunstein90 1990; Ashhab25 2025; Hillery90 1990; Olsen02 2002), this article has studied a feasible model in optical cavities or in a circuit QED system for three-photon squeezing, accompanied by dissipative and driving terms. We should note the differences compared to other previously studied systems (Felbinger98 1998; Bajer91 1991; Denys21 2021), which consider quantized pump accompanied by an applied oscillating term. The specific choices made here were guided by experimental studies for which the pumping and driving modes were modeled in a similar way (Douady04 2004; Bencheikh07 2007; Gravier08 2008; Bertrand25 2025; SandboChang20 2020). Furthermore, we have applied our study to stationary fields generated within cavities (whether optical or superconducting circuit-based), which also distinguishes this work from studies focused on traveling-wave regimes (Douady04 2004; Bertrand25 2025; Dot12 2012; Borshchevskaya15 2015; Okoth19 2019; Corona11a 2011; Corona11b 2011; Cavanna16 2016; Cavanna20 2020; Banic22 2022; Bacaoco25 2025). Thus, we believe we have addressed a gap in the study of third-order nonlinear systems.

Using a mean-field approximation, we found that the system does have a single stable steady state branch restricted to a limited region of its parameters. Outside this region, the system does not have stable states, basically exhibiting an evolution towards divergent values of the observable quantities. The stability region expands as the dissipation factor $\gamma$ and the detuning $\Delta$ increase, shrinking to a region of zero measure for an undamped system at resonance ($\gamma=0$ and $\Delta=0$). On the other hand, an inverse relationship holds for nonlinearity factor $\kappa$. This is consistent with previously established results, which showed that there are no steady solutions for a generalized parametric down-conversion Hamiltonian (Fisher84 1984; Braunstein87 1987; Braunstein90 1990; Hillery90 1990; Olsen02 2002). However, we have shown that even with classical pump, but with dissipative and driving terms, we can obtain finite stable steady results.

We have also considered third-order noise originated from the Kramers-Moyal equation of the Wigner function. Although there are difficulties in interpreting the third-order or higher partial differential equations for quasiprobabilities and converting them to stochastic equations, we applied a method that extends the domain of Wigner function, so that it is possible interpret it as a stochastic process (Drummond14 2014; Drummond17 2017). We have considered that this positive-Wigner function method deserves more attention, because it offers a new way to study stochastic quantum systems via Wigner functions, potentially providing a more accessible physical interpretation. Furthermore, it serves as an additional computational tool for studying nonlinear systems, as demonstrated by the decoupling of calculations regarding dissipation and nonlinear noises.

Applying the stochastic differential equations of the system to the input-output field formalism (Gardiner-Zoller 2004; Collett84 1984; Gardiner85 1985), we have obtained a system of equations that has been linearized. From this, we perform a Fourier analysis of the output modes from the nonlinear cavity. Studying the output quadrature spectra, we were able to verify the existence of quantum noise squeezing with unique features distinct from ordinary squeezing. Unlike second-order squeezing, in three-photon squeezing we observe two symmetric noise peaks away from cavity resonance ($\Delta=0$), which separate two frequency regimes with different behaviors: for smaller detunings, one quadrature is squeezed; whereas for larger detunings the squeezed quadrature is inverted. It is also observed that the noise level reduction is limited and less pronounced in three-photon squeezing. This is due to the triply symmetric geometry of squeezing in phase-space and the squeezing limitation imposed by the Heisenberg inequality. It is also noted that the squeezing effect is greater near the stability boundary of the steady states, at the limit of the driving mode amplitude near the bifurcation, with squeezing being zero when the driving mode is null.

In addition to the spectral analysis associated with the two-time correlation functions, we were able to calculate 2D spectral densities, so that it was possible to observe the dispersive and dissipative aspects and the resonant noise peaks off the cavity resonance of the three-photon squeezing. We have followed the mathematical method developed for magnetic resonance spectroscopy in two frequency dimensions (Ernst-Bodenhausen-Wokaun 1997). Since the 2D spectra are Fourier transforms of the three-time correlation functions of the quadratures centered in their stable state, they are explicitly non-Gaussian, because the third-order cumulants are null for Gaussian systems. Hence we have shown a novel analysis for non-Gaussian and nonlinear systems.

Through our study of the driven dissipative three-photon down-conversion, we were able to deepen our understanding of the third-order squeezing operation, which is a potential process for use in continuous-variable quantum information and computing protocols (Minganti23 2023; Labay-Mora23 2023; Labay-Mora24 2024). Here, we provide further insights for implementing such nonlinear squeezing operations in fault-tolerant computing processes or for constructing a complete set of logic gates (GKP01 2001; Zheng21 2021; Hahn22 2022). Although we focused on the degenerate down-conversion case, this study serves as a starting point for enhancing our understanding of the nondegenerate process. We understand the application of this analysis to the generation of tripartite entangled states as a natural next step. Despite the entanglement has already been the subject of numerous studies (Chekhova05 2020; Agusti20 2020; Zhang21 2021; Tian22 2022; Zhang23 2023; Wei24 2024), our work offers additional insight into stability and spectral analysis.