Regularized linearization for quantum nonlinear optical cavities: Application to Degenerate Optical Parametric Oscillators

Nonlinear optical cavities are crucial both in classical and quantum optics; in particular, nowadays optical parametric oscillators are one of the most versatile and tunable sources of coherent light, as well as the sources of the highest quality quantum-correlated light in the continuous variable regime. Being nonlinear systems, they can be driven through critical points in which a solution ceases to exist in favour of a new one, and it is close to these points where quantum correlations are the strongest. The simplest description of such systems consists in writing the quantum fields as the classical part plus some quantum fluctuations, linearizing then the dynamical equations with respect to the latter; however, such an approach breaks down close to critical points, where it provides unphysical predictions such as infinite photon numbers. On the other hand, techniques going beyond the simple linear description become too complicated especially regarding the evaluation of two-time correlators, which are of major importance to compute observables outside the cavity. In this article we provide a regularized linear description of nonlinear cavities, that is, a linearization procedure yielding physical results, taking the degenerate optical parametric oscillator as the guiding example. The method, which we call self-consistent linearization, is shown to be equivalent to a general Gaussian ansatz for the state of the system, and we compare its predictions with those obtained with available exact (or quasi-exact) methods.


Introduction
Nonlinear optical cavities, that is, cavities containing some element whose response to an applied optical field is nonlinear, are very important both in classical and quantum optics. In the classical domain they allow for effects such as frequency conversion [1] or spatio-temporal pattern formation [2], while in the quantum domain they allow for the generation of quantum correlations manifesting as squeezing or entanglement [3], basic resources for modern applications such as high-precission measurements [4,5,6,7] and quantum information communication and processing [8,9]. The paradigmatic example of such systems are degenerate optical parametric oscillators (DOPOs); these consist in a resonator containing a crystal with second-order nonlinearity, which, when pumped with an external monochromatic laser, is able to generate photons at the subharmonic frequency through the process known as parametric down-conversion (PDC). The interplay between the nonlinear parametric amplification and the cavity losses sets a threshold power below which no subharmonic field is generated at the classical level, and it is close to this critical point where quantum effects are the largest; in particular, more than 90% of quadrature squeezing has been experimentally generated with such a system [10,11,12,13].
In order to analyze the quantum properties of these systems, the simplest and most widely used technique consists in expanding the field as a classical part plus some quantum fluctuations, and linearize the dynamical equations with respect to the latter [14,15]; from another (equivalent) point of view, this approximation means that the state is taken to be Gaussian, with a mean coinciding with the classical field amplitude, and a covariance matrix accounting for the quantum fluctuations. The problem with such approximation is that it gives unphysical predictions such as perfect quadrature squeezing (which requires infinite energy) close to the critical points of the classical theory [16].
Especially for DOPOs, people have developed more refined techniques which go beyond the linear approximation, correcting this unphysical predictions for quantum correlations. Among these techniques, the ones based on the positive P representation [17] are of especial relevance; this representation allows for an exact mapping of the quantum dynamics onto a set of classical stochastic equations from which a proper perturbation expansion or numerical simulation can be carried even at the critical point [18,19,20,21,22]. Moreover, in the limit in which the pump dynamics can be adiabatically eliminated, exact solutions to the steady-state positive P distribution of the DOPO are known [14,23,24]. The DOPO dynamics close to the critical point has even been analyzed via non-equilibrium many-body techniques such as the Keldysh formalism [25,26,27].
The problem with all these beyond-linear techniques is that they are quite complicated when it comes to the evaluation of two-time correlation functions needed for predictions concerning measurements outside the cavity, functions which, on the other hand, are straightforwardly evaluated within a linearized or Gaussian description.
Motivated especially by this last fact, in this work we offer a linear theory (or, equivalently, a Gaussian ansatz for the state of the system) which regularizes the unphysical predictions offered by the usual linearization procedure.
The article is organized as follows. In the next section we describe the quantum model for DOPOs in the Schrödinger and Heisenberg pictures, using, respectively, the cavity modes' master equation and a set of quantum Langevin equations. Then, in Section 3 we introduce our regularized self-consistent linearization procedure in the Heisenberg picture, to show in Section 4 that it is completely equivalent to using a general Gaussian ansatz for the state in the Schrödinger picture. In Section 5 we compare the quantum correlations (quadrature fluctuations) obtained through this method with previous exact (or quasi-exact) methods, showing that, despite its simplicity, it agrees with the latter not only qualitatively, but also quite well quantitatively. In the last section we conclude and comment on other systems where the method could be applied.

The DOPO model
We consider a cavity containing a χ (2) -crystal, pumped with a laser resonant with a cavity mode at frequency 2ω 0 (pump mode), such that photons can be down-converted inside the crystal to the subharmonic resonance ω 0 (signal mode). Denoting byâ p and a s the annihilation operators for pump and signal photons, respectively, the (interaction picture) Hamiltonian which describes such scenario is given byĤ DOPO =Ĥ inj +Ĥ PDC [14,15,16], witĥ where E p is proportional to the amplitude of the injected laser (whose phase is taken as a reference for any other, what allows taking this parameter as real), and χ is proportional to the nonlinear susceptibility of the crystal as well as the overlapping between the spatial modes involved in the down-conversion process. In addition to these coherent processes, we need to introduce the cavity losses; there are two different ways in which this can be done. In the Schrödinger (interaction-)picture, in which the state of the systemρ evolves while operators are fixed, such irreversible processes are accounted for by extra terms in the master equation as [28,24] where the damping rates γ j are proportional to the transmissivity of the coupling mirror at the corresponding frequency.
On the other hand, in the Heisenberg (interaction-)picture, the bosonic operators evolve according to the quantum Langevin equations [28,24] in which the input operators satisfy correlations and account for the vacuum driving fields entering the cavity through the partially transmitting mirror. In the following we explain our regularized linearization of this nonlinear system in both pictures, since they provide different intuitive ideas of what the procedure means.

Heisenberg picture approach: self-consistent linearization
Before explaining the procedure, let us define the following dimensionless parameters and normalized variables in terms of which the quantum Langevin equations (3) are rewritten as note that the normalized input operators satisfy the correlations (4) but now with respect to the new dimensionless time τ . In order to linearize these equations, we expand the annihilation operators aŝ b j = β j + δb j , with β j some amplitudes which one identifies with the mean-field part of the modes b j , and δb j the operators accounting for quantum fluctuations around such mean-field, with respect to which the theory will be linearized. The mean-field amplitudes can be evaluated by taking the quantum expectation value of the quantum Langevin equations (7), leading to 1 κ In the usual approach, these mean-field equations are solved by assuming that the state is coherent in all modes, hence neglecting the δb 2 s and δb p δb † s terms, what gives rise to the classical equations that would have been obtained from Maxwell's equations; in other words, the mean-field amplitudes β j are taken to be the classical solutions of the system [14,15]. In particular, in the case of equations (8), this coherent mean-field ansatz provides two different types of steady-state solutions depending on the injection parameter σ (see the thin-solid light-grey curve of Figure 1): one known as the belowthreshold solution with (β s = 0,β p = σ) which is the only stable solution for σ ≤ 1, and another (bistable) solution for σ > 1 known as the above-threshold solution with the signal field switched on (β s = ± √ σ − 1,β p = 1). The injection σ = 1 marks a critical point where the below-threshold solution changes from stable to unstable, and, as we argue below, this sudden change in the stability conditions is what generates unphysical predictions in the linear theory [14,15,16]. In order to correct such a problem, we propose to incorporate some information of the quantum dynamics in the quantum-fluctuations-dependent mean-field equations (8), with the purpose of obtaining a better ansatz for the mean-field amplitudes. Concretely, using (8), the quantum Langevin equations are rewritten in terms of the fluctuations δb j as We can proceed then as in the regular linearization by neglecting the nonlinear fluctuations δb 2 s − δb 2 s and δb p δb † s − δb p δb † s , obtaining the linear system , and where the so-called linear stability matrix is defined as the difference now is that, instead of substituting the classical solutions which give rise to a singular stability matrix at threshold (what can be traced as the source of all the unphysical predictions ‡), the mean-field amplitudes β j are left as unknown variables which are found self-consistently by calculating the δb 2 s and δb p δb † s terms from this linear system, and plugging them into the mean-field equations (8). In this case, it is simple but lengthy, for example by finding the (bi-orthonormal) eigensystem of (11), to obtain the following steady-state expressions for these correlators where we have introduced what we call the mean-field intensities I j = |β j | 2 , complemented with the phases ϕ j ∈ R defined from β j = I j exp(iϕ j ). Introducing ‡ In particular, the singular character of the linear stability matrix occurs because one of its eigenvalues becomes zero at σ = 1 (it changes from negative to positive, corresponding to the destabilization of the below-threshold solution when crossing the threshold); on the other hand, it is clear that within this linear description, the correlators of quantum fluctuations are inversely proportional to combinations of these eigenvalues, and hence some might diverge at threshold, what we will show later explicitly.  [19,20] (see Section 5), while in (b) it corresponds to the number of (normalized) signal photons obtained from our self-consistent below-threshold solution, showing how it is not divergent at threshold, in contrast with the predictions found with the usual linearization approach.
these expressions into the mean-field equations (8) in the stationary limit (β j = 0), we get where the bar denotes 'steady-state values'. It is then straightforward to show that the pump phaseφ p is locked to 0, while the signal phaseφ s can take the values 0 or π, just as in the classical solution. As for the intensities, equation (14) yields a third order polynomial in I s , with a trivial root I s = 0 (below threshold solution) and second root (above threshold solution) which can be written in terms of the pump intensity as with the third root is not relevant, since it can be checked that it leads to unphysical results, as commented below. Finally, equation (13) gives an equation for the pump intensity, which can be written as below threshold (Ī s = 0), this gives a third order polynomial in Ī p whose roots can be found analytically (although not much insight is gained from their complicated expression, so we don't give them explicitly), but above threshold it gives a high-order polynomial whose roots we've only been able to find numerically. Of the many solutions forĪ p obtained from this equation, most of them are disregarded because they are either negative or complex, makeĪ s in (15) negative or complex, or lead to negative photon numbers â † jâ j or quadrature fluctuations incompatible with the Heisenberg uncertainty principle δx 2 j δŷ 2 j ≥ 1 (see below for a definition of the quadrature fluctuations). It is quite remarkable that, after discarding all these unphysical solutions, only two solutions forĪ p remain: one of the three roots of the below-threshold polynomial- (17) with I s = 0-, and another one from the above threshold one-(17) withĪ s given by (15)-.
In Fig. 1 we plot the pump (a) and signal (b) intensities associated to these solutions as a function of the injection; it is apparent that they tend to the classical solutions far away from the critical point, but never reach, in particular, the value of the intensities that makes the linear stability matrix become singular (Ī s = 0,Ī p = 1). Hence, the solutions of our self-consistent linearization can be seen as regularized versions of the classical below-and above-threshold solutions, which remove the divergences of the linear theory of quantum fluctuations (see Section 5 for and explicit discussion of the regularized quantum properties). The figures also allow us to see the way in which this regularization occurs: in the case of the below-threshold solution (solid blue curves), the pump intensity does not grow quadratically with the injection as happens with the classical solution (thin-solid light-grey curve), staying below its critical value I p = 1 for any physical value of the injection σ; as for the above-threshold solution (dasheddotted red curve), instead of connecting continuously with the below-threshold solution as the classical solution does, its appearance is delayed a little bit respect to the classical threshold σ = 1, and starts with a nonzero signal intensity I s . It might seem strange that the above-and below-threshold solutions do not connect continuously, but in Section 5 we will argue why the presence or not of this 'jump' is quite irrelevant indeed, as can be intuitively understood from the point of view of the symmetry breaking that occurs above threshold: too close to threshold quantum tunneling between the solutions with opposite signal phase is too fast, and it makes no sense to talk about these solutions independently.
In summary, in this section we have provided a self-consistent linearized theory for the DOPO in which the mean-field amplitudes are found not from the classical nonlinear equations of motion, but from ones including quantum corrections; this method provides regularized versions of the solutions that would be obtained from the bare classical theory, which avoid in particular the critical values of the mean-field amplitudes at which the linear stability matrix becomes singular, thus avoiding the unphysical results related to it. Let us remark that this self-consistent method that we have put forward was already introduced in [27] for the DOPO problem, but only below threshold; moreover, in the next section we give full meaning to the method by showing that, within the Schrödinger picture, it is equivalent to making a general Gaussian-state ansatz consistent with the master equation.

Schrödinger picture approach: general Gaussian ansatz
The self-consistent linearization introduced in the previous section admits a simple interpretation from the point of view of the state of the system: we argue in the following that it is equivalent to making a general Gaussian ansatz for it. In order to show this, our starting point is the master equation (2), which in terms of the normalized parameters (5) and variables (6) can be rewritten as where we have defined an anti-hermitian operatorÂ s )/2. The evolution equation for the expectation value of any operatorB is obtained as Applied to the annihilation operatorsb j , this equation provides exactly the mean-field equations (8) found in the previous section for the amplitudes β j = b j , which depend on the second moments δb 2 s and δb p δb † s . On the other hand, the evolution equations of these second moments depend on third-order moments, and here is where the Gaussianstate approximation enters into play: we assume that the state of the system is Gaussian [9] at all times, meaning that higher order moments can be written as products of first and second moments only. In particular, this has the consequence that third order moments of quantum fluctuations vanish, that is, δb m j δb n k δb p l = 0, where m, n, and p can be either dagger or nothing. Under this assumption, the evolution equations of the second order moments , δb p δb † s * , δb 2 s , δb 2 s * , δb † s δb s ), form the following closed linear system dm dτ = M(β p , β s )m + n(β p ), (21) Figure 2. Variances of the squeezed (a) and anti-squeezed (b) quadratures as a function of the injection parameter σ, for κ = 1 and g = 0.01 (similar figures are found for any other choice). As in Fig. 1, the thin-solid light-grey curve corresponds to the usual linearized description; the solid blue line and dashed-dotted red line correspond to our self-consistent method below and above threshold, respectively; and the dashed yellow line corresponds to the predictions of Drummond and collaborators' perturbative analysis [19,20].
with the matrix M given by and n = g 2 col(0, 0, 0, 0, 0, 0, 0, β p , β * p , 0); (22) it is simple to show that the solutions to this linear system coincide with the moments obtained from the linearized quantum Langevin equations (10), in particular steadystate moments such as (12), and hence, the self-consistent linearization introduced in the previous section is strictly equivalent to the assumption that the state of the system is a general Gaussian state whose moments satisfy the constrains imposed by the master equation. Note that, since we have three possible solutions for the meanfield amplitudes β j , the below-threshold solution withβ s = 0 and two above-threshold solutionsβ s = ± Ī s , we have three Gaussian ansatzes that we can use, which we will denote byρ G,0 andρ G,± , respectively. In the next section we interpret the meaning of these solutions, as well as the jump observed above threshold for the signal intensity, which was partially discussed in the previous section.

Analysis of the results and comparison with previous methods
Let us now analyze the predictions that this self-consistent linear theory makes for the squeezing of the intracavity signal field, and compare it with some known beyond-linear results. Let us then define the quadraturesx s =â † s +â s andŷ s = i(â † s −â s ), and the corresponding fluctuation operators, δx s =x s − x s and δŷ s =ŷ s − ŷ s ; it is simple from (10) or (21) to obtain the following expression for their variance in the steady-state: When the classical solutions forĪ p andĪ s are considered, one immediately sees that δx 2 s → ∞ and δŷ 2 s → 0.5 at threshold σ = 1, which is exactly the unphysical prediction that we were talking about, since infinite quadrature fluctuations imply infinite photon number, that is, â † sâ s → ∞. On the other hand, when we introduce in these expressions the regularized solutions discussed in Section 3, obtained through the self-consistent method, δx 2 s becomes finite in all parameter space, while at the same time, the squeezing level of δŷ 2 s is reduced; this can be appreciated in Fig. 2. More quantitatively, exactly at the critical point σ = 1, it is simple to obtain the following expressions for the quadrature variances to the leading order in g: , and lim t→∞ δŷ 2 s ≈ 0.5 + showing explicitly how the anti-squeezing is regularized, while the variance of the squeezed quadrature is increased with respect to its value obtained with the usual linearization procedure. In order to understand how good the self-consistent linearization is from a quantitative point of view, we now compare these results with the ones obtained from the perturbative approach that Drummond and collaborators developed in the vicinities of the critical point, by making a consistent multiple-scale expansion of the stochastic variables within the positive P representation [19,20]. This procedure has the virtue of being valid for any κ, although it is reliable only close enough to threshold, concretely for |σ − 1| < g/ √ 2; nevertheless, since we are mainly interested in how the self-consistent linearization regularizes the conventional linearization around the critical point σ = 1, where the divergences appear, this will be enough for comparing with our results. For our purposes, their most relevant results concern the steady-state (normalized) pump amplitude and the quadrature fluctuations of the signal field, which read to the leading order in g, where we have defined a real stochastic variable x distributed according to the probability density function with d a suitable normalization constant. The square of (25) corresponds to the curve that we plotted in Fig. 1(a) to compare with the steady-state intensityĪ p that our self-consistent linearization provides. On the other hand, exactly at threshold (σ = 1) the moments derived from this distribution admit very simple expressions in terms of Gamma functions Γ(z), and in particular we have together with (24), allows us to write From these expressions we see that the self-consistent linearization provides a regularization which compares pretty well with Drummond's predictions, at least regarding the order of magnitude. In particular, we bring the reader's attention to the anti-squeezing predicted by the self-consistent method, which is only 50% above the one predicted by Drummond and collaborators. As a last test, we now compare the Gaussian ansatzes proportioned by our selfconsistent linearization against the exact state of the signal field, which is known in the limit κ ≫ 1 (in which the pump field can be adiabatically eliminated) [14,23,24]; in particular, the positive P distribution associated to the reduced steady state of the signal mode is given by [23,24] where α s and α + s are real, and K is a suitable normalization factor. Explicit expressions of the stateρ can be built from such a positive P distribution, but, for our purposes, we only need the fact that it allows for the evaluation of steady-state moments in normal order as Let us define the variables x + = α s + α + s and x − = α s − α + s , noting that x + corresponds directly to the stochastic representation of the anti-squeezed quadraturex s =â s +â † s , while x − corresponds to 'i-times' the squeezed one, that is, to iŷ s =â s −â † s . In Figs. 3(a-e) we show how the marginal p(x + ) = R dx − P (x + , x − ) changes as we cross the threshold. The complementary marginal q(x − ) = R dx + P (x + , x − ) is shown only at one value of σ because it does not change perceptibly around the critical point, and the Gaussian ansatzes adapt to it almost perfectly, as can be appreciated in Fig. 3(f).
The first physically relevant result that this exact steady-state solution predicts is â s = 0 for all σ. This might seem surprising, since it seems to suggest that the signal field is never switched on, that is, that the above-threshold solution with â s = 0 characteristic of DOPOs is incorrect; this, however, is not true, the right answer is a bit more subtle: as σ is increased, the distribution develops two peaks (see Fig. 3) which, individually, correspond to the above-threshold amplitudesβ s = ± Ī s , but, since they are developed together, their contributions to â s average to zero. From a more fundamental point of view, this just reflects the fact that the master equation (2) has the Z 2 symmetryâ s → −â s , and hence, ifρ sol,+ is a symmetry-breaking ansatz with â s = Ī s /g, thenρ sol,− = exp(−iπâ † sâ s )ρ sol,+ exp(iπâ † sâ s ), which has â s = − Ī s /g, is another equally valid ansatz; in the absence of any bias, the state of the system has to be regarded as the balanced mixtureρ sol = (ρ sol,+ +ρ sol,− )/2, which is the one giving the correct experimental statistics: every time the DOPO is switched on, it has to choose the phase 0 or π according to the particular initial fluctuations from which the steady state is built up, that is, according to spontaneous symmetry breaking. It feels natural to think that one can force the system to pick one particular phase at every experimental run (at least in a metastable sense) by adding an explicit symmetry breaking mechanism such as the injection of a weak laser at the signal frequency-a term likeĤ inj,s = i E s (â † s −â s ) in the Hamiltonian-; however, this picture is only correct once enough above threshold, since close to threshold quantum tunneling between the stateŝ ρ sol,± is too fast [14,29,30], and no phase locking can be achieved within the observation time. In other words, from an observational point of view, it only makes sense to analyze the properties ofρ sol,+ andρ sol,− separately once the peaks of the positive P distribution are enough far apart, and this is why we made the statement that the jump seen in the signal intensity above threshold with the self-consistent method is not relevant for real applications, as it just reflects the fact that some distance from threshold is required forρ sol,± to have independent meaning, as otherwise fast tunneling times prevent their existence.
Knowing the exact form of the positive P distribution in the limit κ ≫ 1 allows us to get a pictorial feeling of how good our Gaussian ansatzes are. In particular, superposed to the exact marginals (solid-thin light-grey line), in Fig. 3 we plot the marginals corresponding to the Gaussian ansatzesρ G,0 (solid blue) andρ G,> = (ρ G,+ +ρ G,− )/2 (dashed-dotted red). We can appreciate howρ G,0 adapts very well toρ sol below threshold, except close to the point in which the peaks start developing (σ = 1 − g 2 /4), where the exact distribution flattens in the center, loosing its approximate Gaussianity. Note also that the above-threshold Gaussian ansatzρ G,> appears at σ = 1 + g in this κ → ∞ limit, and it quite rapidly converges to the exactρ sol as we increase σ.

Conclusions and outlook
In summary, we have developed a linearization procedure for nonlinear optical cavities in which the mean-field amplitudes are found self-consistently by introducing some information about the quantum fluctuations in their evolution equations; we have applied it to one of the simplest nonlinear resonators, the DOPO, showing how the method is capable of regularizing the divergences found at the critical point with the traditional linearized analysis. We have also shown that such procedure is completely equivalent to using a general Gaussian ansatz for the state of the system. Finally, we have compared the results derived from the self-consistent linearization with other known exact (positive P distribution under adiabatic elimination of the pump [23,24]) or quasi-exact (Drummond and collaborators' perturbative expansion [19,20]) methods for the DOPO, proving that they are reliable qualitatively, and also roughly quantitatively, with the advantage that two-time correlators (and hence output spectra of observables) are straightforwardly evaluated within the linearized description, what we think will be useful to study quantum correlations in more complicated systems such as nondegenerate or multi-mode OPOs [31,32,33,34]. The application of the method to systems with other types of critical points and bifurcations is an open question that will be interesting to analyze in the future as well; the bistability found in Kerr resonators [35,36] or the Hopf bifurcations found in optomechanical cavities [39] and intracavity second-harmonic generation [37,38], are good candidate scenarios.