# Self-Consistent Projection Operator Theory in Nonlinear Quantum Optical Systems:

A case study on Degenerate Optical Parametric Oscillators

###### Abstract

Nonlinear quantum optical systems are of paramount relevance for modern quantum technologies, as well as for the study of dissipative phase transitions. Their nonlinear nature makes their theoretical study very challenging and hence they have always served as great motivation to develop new techniques for the analysis of open quantum systems. In this article we apply the recently developed self-consistent projection operator theory to the degenerate optical parametric oscillator to exemplify its general applicability to quantum optical systems. We show that this theory provides an efficient method to calculate the full quantum state of each mode with high degree of accuracy, even at the critical point. It is equally successful in describing both the stationary limit and the dynamics, including regions of the parameter space where the numerical integration of the full problem is significantly less efficient. We further develop a Gaussian approach consistent with our theory, which yields sensibly better results than the previous Gaussian methods developed for this system, most notably standard linearization techniques.

###### pacs:

42.50.-p, 03.65.Yz, 42.65.Sf, 42.65.Yj## I Introduction

Nonlinear optical systems play an important role in the field of optics both in classical Boyd2003 (); Staliunas2002 () and in quantum WMbook (); Scullybook (); Meystre1991 (); Hartmann07 (); Brandao08 () regimes. Quantum mechanical effects, in particular, which are not explainable by classical optics, have triggered substantial research, especially in connection to modern applications such as high-precission measurements Goda08 (); Vahlbruch05 (); Treps03 (); Treps02 () and quantum information communication and processing Braunstein05 (); Weedbrock12 (). Importantly, the nonlinear nature of these systems leads to non-Gaussian states, which typically precludes an analytic treatment and therefore requires elaborate theoretical approaches CarmichaelBook1 (); CarmichaelBook2 ().

In a system where the dynamical degrees of freedom evolve on different time scales, approximate descriptions of reduced complexity may be found. For example, adiabatic elimination techniques can be exploited to derive effective equations of motion Mori65 (); ZollerBook (). In this work, we apply the recently introduced self-consistent projection operator theory Degenfeld14 () to the degenerate optical parametric oscillator, and exemplify how it generalizes adiabatic elimination approaches. This theory takes dynamical back-action between the degrees of freedom into account and therefore does not require any time-scale separation. We expect our method to be directly applicable to other nonlinear quantum optical models such as those for nondegenerate or multi-mode parametric oscillation Reid88 (); Drummond90 (); Navarrete08 (); Navarrete09 (), lasing WMbook (); Scullybook (); MandelWolf (); Breuer07 (), optomechanical parametric oscillation OMPOCarlos (); Nori (), or the dissipative Dicke model DallaTorre (); Parkins (); Ritsch ().

Degenerate optical parametric oscillators (DOPOs) have been extensively studied in the past Meystre1991 (); CarmichaelBook2 () and are one of the paradigm examples of a system subject to a driven and dissipative phase transition. It is formulated as a bosonic problem with two modes, signal and pump, subject to dissipation and interacting nonlinearly. In the adiabatic limit of a fast decaying pump mode, an effective master equation can be derived by means of standard projection operator approaches CarmichaelBook2 () and due to its reduced complexity, the steady state can be found by solving the corresponding Fokker-Planck equations for the positive P distribution Wolinsky88 (); Drummond80 (). Yet away from the adiabatic limit one has to resort to numerical simulations or perturbative treatments Gardiner80 (); Kinsler93 (); Kinsler95 (); Chaturvedi02 (); Chaturvedi99 (); Pope00 (). Non-equilibrium many-body techniques such as the Keldysh formalism have also been employed to study steady-state properties FleischhauerKeldysh (); Mertens93 (); Swain93 (). While the application of all these techniques has allowed to deepen our understanding of DOPOs and phase transitions in driven dissipative quantum systems enormously, it is important to note that they are naturally built to determine the evolution of observables, making the determination of the quantum state of the optical fields very challenging, if not impossible.

Our approach, in contrast, derives a set of coupled equations for the reduced states of the two optical modes of the DOPO. By numerically solving these equations, we find the reduced density matrices of both the pump and the signal modes. We test the accuracy of our method by comparing its results with those of the full DOPO problem in regions of the parameter space where this is numerically tractable. Our findings show that our method is remarkably close to the exact results, both for steady states and dynamics, while being less numerically demanding than the full simulation of the DOPO problem. It thus gives access to the reduced states of the cavity modes in regions of the parameters that are inaccessible to the latter.

The possibly largest reduction of complexity in nonlinear quantum optical systems, however, comes from the application of Gaussian approximations on the state of the system. Within a Gaussian theory one can basically cover the whole parameter space efficiently to determine both steady-state and dynamical quantities such as two-time correlation functions. The simplest and most widely used Gaussian approach is known as the linearization technique Drummond80 (); Lugiato81 (), which consists in assuming that the system configuration is, on average, in its classical state, but is constantly driven out of it by some “small” quantum fluctuations. While this technique provides a good qualitative picture of the physics in many, albeit not all, systems, it leads to unphysical predictions close to the critical points of the classical theory, e.g., to infinite photon numbers in the case of the DOPO Collet84 (). These unphysical predictions can be regularized by applying a more elaborate Gaussian state approximation where the system is not forced to stay in its classical state, but chooses instead an average configuration more consistent with the quantum fluctuations that perturb it Carlos14 (). Motivated by such an idea, we apply a Gaussian approximation within the self-consistent projection operator theory, and show that it gives more accurate quantitative results than any of the usual Gaussian techniques, as it does not assume a Gaussian state for the entire system, but only for the reduced state of one of the modes.

The remainder of the paper is organized as follows. In Sec. II we introduce the DOPO model. We also discuss its symmetries and briefly elaborate on the standard linearization approach in Sec. II.1. Sec. III reviews the main concepts of the self-consistent projection operator theory and introduces the self-consistent Mori projector (c-MoP) equations, which lie at the center of our study. Our theory provides a systematic extension of mean-field approaches as demonstrated in Sec. III.1 and reproduces known results in the adiabatic and the diabatic limits introduced in Sec. III.3. An efficient procedure designed to deal with the non-Markovian structure of the c-MoP equations is provided in Sec. III.2, which we use in Sec. IV to test the accuracy of our method for steady-state quantities and to present quantum states of the signal mode. A Gaussian state approximation on the c-MoP equations is performed in Sec. V, which is shown to lead to highly accurate quantitative results as compared to previous linearization techniques. As a further test, we check in Sec. VI that our method provides the same level of accuracy for the dynamics, as it does for steady states. Finally, we conclude our work and present an outlook in Sec. VII.

## Ii The degenerate optical parametric oscillator

A DOPO consists of a driven optical cavity containing a crystal with second order optical nonlinearity, see Fig. 1. Two relevant resonances at frequencies (signal mode) and (pump mode) exist in the cavity, which are nonlinearly coupled via parametric down-conversion inside the crystal, capable of transforming a pump photon into a pair of signal photons, and vice versa. We assume that the external driving laser is resonant with the pump mode. Including damping through the partially transmitting mirrors at rates and for the pump and signal modes, respectively, the equation governing the evolution of the state of the system in a picture rotating at the laser frequency is given by Meystre1991 (); CarmichaelBook2 ()

(1) |

where is the down-conversion rate and is proportional to square root of the injected laser’s power. We have defined bosonic operators and for the pump and signal modes, respectively, which satisfy canonical commutation relations and . Note that the nonlinear interaction is third order in the field operators, precluding a general analytic solution of Eq. (1) to which we refer as the Liouville-von Neumann equation or simply the full master equation of the DOPO.

### ii.1 Linearization approach and symmetry breaking

The right hand side of Eq. (1) can also be written in a shorthand notation by introducing a superoperator (Liouvillian), such that . For the major part of this work, we will be interested in the steady state , which fullfills the equation . Due to the dissipation acting on both modes and because an arbitrarily large but finite truncation will always provide an arbitrarily good approximation, we expect the steady state to be unique Schirmer (); Rivas12 ().

We further note the invariance of the Liouvillian under a unitary transformation of Ising-type which transforms as . Since the steady state is unique, this implies that it has to be invariant under the transformation too, i.e. . This in turn leads to vanishing steady state expectation values which include odd powers of the signal field operator . In particular , as for example .

However, the most common technique used to analyze Eq. (1), known as the linearization approach, breaks this symmetry Drummond80 (); Lugiato81 (), which has to be restored “by hand” at the end of the calculation, following the procedure that we explain at the end of Sec. IV. Even though this method is more naturally introduced in the Heisenberg picture using the language of quantum Langevin equations, it also admits a Schrödinger picture interpretation in terms of two successive approximations in the master equation. It starts by writing the bosonic operators as , with and hence . In the first approximation, the fluctuation operators are neglected altogether; the evolution equations for (Bloch equations), then provide a set of nonlinear differential equations for the amplitudes , which in the case of the DOPO read

(2) |

These correspond to the classical equations of the system, as they could have been obtained directly from Eq. (1) by assuming a coherent state for , or simply from Maxwell’s equations. Depending on the injection parameter one finds two types of steady-state solutions of Eq. (2). One of them has and , and hence it does not break the symmetry; it is known as the below-threshold solution, and is only stable for . The other solution is bistable and has and , hence breaking the symmetry; it is known as the above-threshold solution, and exists only for . The threshold point marks a critical point where the classical theory predicts a phase transition from a signal-off phase with to a signal-on phase with . In the signal-off phase all injected power goes into the pump mode, while after crossing the critical point all the extra injection is transferred to the signal mode, see the gray thin solid line in Fig. 2.

Once the classical solutions have been identified, the second approximation consists in coming back to the original master equation with the bosonic operators written as , and neglecting any term which goes beyond quadratic order in the fluctuation operators . This leads to a so-called linearized master equation which can be easily solved.

One has to keep in mind that this linearized theory can only be trustworthy when the classical solution is a strong attractor, because only then the quantum fluctuations driving the system out of equilibrium are strongly damped, and quantum noise can be treated as a small perturbation. This means that, in particular, any predictions obtained through this method cannot be trusted in the vicinities of critical points of the classical theory: points of the parameter space where one solution becomes unstable, making way for a new solution to kick in, hence creating non-analytic behaviour in some observable, that is, a classical phase transition. Indeed, this is exactly the case for the DOPO, in which this linearized description breaks down at threshold, offering unphysical predictions such as infinite photon numbers in the signal field (as illustrated by the gray thin line in Fig. 4).

## Iii Self-consistent Mori Projector Approach

To explain the approach employed in our calculations, we will first recapitulate some basic ideas of the self-consistent projection operator theory Degenfeld14 (). The first step is to divide the entire system into subsystems. In the DOPO this naturally amounts to consider the pump mode described by its reduced state and the signal mode described by . In the spirit of open system theory Zwanzig01 (); Breuer07 () we will first treat the pump mode as an “environment” for the signal mode, which then takes the role of the open “system”. Technically this is done by introducing the time-dependent, self-consistent Mori projector whose action on the full state gives the factorized state . The term “self-consistent” is chosen because the state of the pump in is not a time-independent reference state but is rather obtained consistently from the time-evolving state of the full dynamics. Using this projector, we derive a generalized Nakajima-Zwanzig equation which is an exact equation for the reduced state of the signal mode Degenfeld14 (). The effective Liouvillian describing such a Nakajima-Zwanzig equation will depend on the state of the pump . In order to obtain a closed set of equations we need to reverse the scenario and treat the pump mode as the “system” and the signal mode as the environment, see Fig. 1 for an illustration.

Again analogous to open system theory, we split the full Liouvillian from Eq. (1) into three parts. After performing a displacement , where will be chosen later, see Sec. III.1, we write , with

(3) |

where we have defined the standard Lindblad superoperator , with being an arbitrary operator. The displacement moves the large coherent background of the pump field into the free evolution of the signal , keeping only the pump mode’s fluctuations within the nonlinear signal-pump interaction . Such a step is important as our theory expands in powers of the interaction Liouvillian in order to solve the Nakajima-Zwanzig equation. As in reference Degenfeld14 () we will expand to second order in the system-environment interaction. This approximation is known as the Born approximation Breuer07 (). The effective equations of the signal and the pump mode then read,

(4) | ||||

(5) | ||||

where we have defined the Kernel superoperators

(6) | ||||

(7) | ||||

and for any operator acting on the signal () or pump () subspace, we have defined the corresponding fluctuation operator .

The state of the pump mode enters the signal mode’s dynamics, eq. (4), via and the correlation functions

(8) | |||

In turn, the state of the signal mode enters the pump mode’s dynamics, eq. (5), via the expectation value and the correlation functions

(9) | |||

Equations (4) and (5) should be understood as two coupled equations which represent effective equations for the reduced states of the signal and the pump mode. We refer to these two equations as the c-MoP (consistent Mori Projector) equations of the DOPO. They can be thought of as non-Markovian and nonlinear master equations which do not rely on any time-scale separation between the modes. We will elaborate in detail on the limits where time-scale separation is present in Sec. III.3.

The only assumptions made so far are the Born approximation and the assumption of an initially factorized state . The latter seems very reasonable by considering the vacuum as the state of the two modes before the driving laser is switched on. We also emphasize, our approach does not ignore system-environment or rather signal-pump correlations. In fact, it has been shown Degenfeld14 () that the Born term, the term second order in which is here proportional to , clearly takes signal-pump correlations into account. We will show the crucial importance of the Born term in several examples below. Of course, c-MoP theory or any theory based on the concept of projection operators does not give access to explicit expressions for system-environment correlation functions. An example in this context could be the cross-correlation function .

The most striking advantage of projection operator theories and in particular of the c-MoP theory is the reduction of the complexity of the problem. In the example of the DOPO the complexity of the Liouville-von Neumann eq. (1) scales as , where denotes the Hilbert space of the signal/pump modes, while the complexity of the c-MoP equations scale as . The self-consistent Mori-projector theory thus offers a very significant reduction of complexity.

### iii.1 Mean-field Approximation

A merely approximate but very simple way of solving the c-MoP equations is to consider all terms up to first order in the interaction only. Hence we drop all terms proportional to from eqs. (4) and (5). Within this approximation it does not make a difference whether the pump field is displaced or not. For simplicity we put the displacement from eq. (3) to zero and obtain two coupled equations

(10) |

known as mean-field equations FleischhauerMF (). These equations are quadratic in the field operators and therefore it is straightforward to solve them either numerically for the dynamics or analytically for the fixed points FleischhauerMF (); Carlos14 (). The stationary state of the signal mode will be a Gaussian state Braunstein05 (); Weedbrock12 (); CarlosQI () centered around a vanishing field amplitude as the mean-field equations do not break the Ising-type symmetry. The steady state of the pump mode will be a coherent state with an amplitude given by .

Just like the c-MoP equations (4) and (5), the mean-field equations are coupled nonlinear equations which have to be solved self-consistently. Within mean-field theory fluctuations of the pump mode are disregarded. Fluctuations of the signal mode, however, are (at least to some extend) taken into account FleischhauerMF (); Carlos14 (). This leads to the regularization of the divergences appearing in the classical theory or rather the standard linearization approach. For our purposes it is important to note that the pump field amplitude always stays below the classical above-threshold solution, i.e. . In the remainder of the paper we will use it as the displacement in eq. (3), i.e. . This will guarantee a well-behaved Liouvillian for the free system as we will explain in more detail in Sec.III.2.

The mean-field equations can also be found by putting the factorized state Ansatz into the Liouville-von Neumann equation, here given by eq. (1), before tracing out each of the modes separately. This well-known procedure, indeed, neglects all signal-pump correlations. Within the self-consistent projection operator theory, mean-field can be understood as an approximation to linear order in the interaction for the dynamics of reduced density matrices. Our theory therefore provides a systematic generalization of mean-field approaches. It is due to the Born terms, which are second order in , that signal-pump correlations are taken into account. Therefore, we expect a different quality of approximation by going from first order to second order in the interaction.

### iii.2 Born terms

In order to solve the full c-MoP equations including the Born terms we will need to overcome two main difficulties. While the c-MoP equation (5) of the pump mode is quadratic in the field operators, granting us with a closed set of equations including only first and second moments of the pump field, the c-MoP equation (4) of the signal mode is quartic in the field operators. We will therefore either solve the equation of the signal fully numerically, see Sec. IV, or apply a Gaussian state approximation as presented in section Sec. V. In any of these two approaches, we need to overcome the second difficulty which arrises due to the non-Markovian structure of our theory. In the remainder of this section we will thus show how to rewrite an integro-differential equation of first order into a set of coupled ordinary differential equations. For the present problem this step is crucial, as solving the integro-differential equations is significantly more demanding for both numerical and analytical approaches.

We start by evaluating the correlation functions of the pump. By taking derivatives of the pump correlators and with respect to , see eq. (8), considering initial conditions at (note that we understand from the c-MoP equations that ), and exploiting the fact that the operator is traceless, we find

(11) |

Hence, all correlation functions of the pump can be written in a form where the dependence only enters in a simple exponential factor.

A bit more effort is needed in order to simplify the correlation functions of the signal, but the main steps are mainly identical. All the functions in eq. (9) are of the form with a traceless operator depending solely on . Again, we take the derivative of with respect to and find an equation of motion of the form with a column vector

(12) |

where the expectation values with the tilde are defined in the usual way as the trace over the signal mode but with a density matrix given by . The matrix reads

It is straight forward to diagonalize . We write , with a similarity matrix that can be found analytically (but its expression is too lengthy to be reported here), and is the diagonal form of containing its eigenvalues , and . We now solve for the vector , to find

(13) |

where we have defined the initial condition vector

(14) |

and the matrices , where is a projector in the ’th “direction”, that is, a matrix with zeros everywhere but in element which is one.

Note that for the limit to be uniquely defined, and therefore for to be well-behaved, all the eigenvalues of must satisfy , which in turn leads to a requirement for the displacement . This requirement is indeed fulfilled by choosing the mean-field displacement as mentioned above, see Sec. III.1. In contrast, taking the classical solution as the displacement would lead to an ill-behaved above and at the classical threshold point, that is, for .

Coming back to the correlation functions in eq. (9), the general solution (13) allows us to write them all as

(15) |

with (the subscript denoting the third vector component), where denotes any of the correlation functions , for which is taken, respectively, as . Let us emphasize that, just as with the pump mode, we have been able to write all the correlation functions of the signal mode into a form where the dependence only enters in simple exponential factors.

Finally, let us show how this form for the correlation functions allows us to turn the c-MoP equations, which are coupled integro-differential equations, into coupled ordinary differential equations. For this aim, let us rewrite eqs. (4) and (5) as

(16) | ||||

(17) | ||||

where we have defined the operators

(18) |

with the superoperator defined as in eq. (7), but with the correlation functions instead of . Using their definition, and the solutions found for the correlation functions, eqs. (11) and (15), their evolution equations are found to be

(19) | ||||

(20) |

Together with eqs. (16) and (17), these form a closed set of coupled nonlinear ordinary differential equations for the reduced states and , and the traceless operators and . These are the equations that we analyze in the remainder of the paper.

Overall we have shown for the example of the DOPO that it is indeed possible to rewrite the integro-differential c-MoP equations into a set of ordinary differential equations. The steps presented here are quite general and can be pursued for all c-MoP equations describing any physical system. The complexity of the resulting set of coupled equations will depend on the complexity of the subparts of the full quantum system, here given by the complexity of and .

Finally, we remark that the c-MoP equations preserve the trace and the hermiticity but they do not guarantee for the positivity of the density matrix. Such an issue is not unusual for projection operator theories, in fact, the same conditions can be found in the well established Redfield equations Redfield (); Blum (). Obviously whenever the c-MoP equations provide a good approximation, they will yield a positive density matrix. Hence the positivity of the eigenvalues can be used as a consistency test for the approximation.

### iii.3 The adiabatic and the diabiatic limit

The two dissipation rates and set a time scale on which the pump and the signal, respectively, relax to the steady state of their unperturbed Liouvillians and . In standard open system theory one relies on a separation of time scales between the system dynamics and the environment correlations. A similar reasoning is applied in adiabatic elimination approaches, where for the DOPO one relies on a time scale separation between signal and pump. The c-MoP theory can, in fact, be understood as a generalization of adiabatic elimination procedures where one has to consider the back-action of the “system” onto the “environment”. We will now show that the effective equations for the reduced state of the signal known in the adiabatic CarmichaelBook2 () and the diabatic FleischhauerMF () limit can, indeed, be obtained as limiting cases of the c-MoP equations.

The adiabatic limit in which the time scale of the pump mode is much faster than the time scale of the signal mode is defined such that while is kept finite. The diabatic limit describes the opposite scenario where . We proceed by comparing the Born terms with the free evolution operators and , for which we consider the scaling of and , which can be obtained by simple inspection of eqs. (19) and (20) divided by and , respectively.

In the adiabatic limit, we infer from eq. (20)/ that for all and . Introducing this result into eq. (17), we see that the state of the pump will be coherent with a field amplitude obeying the equation of motion

(21) |

On the other hand, eq. (19)/ leads to , where we have used eqs. (6) and (11) and the fact that when the pump is in a coherent state all the expectation values in eq.(11) cancel. Introducing this result into eq. (16), together with the steady-state solution of eq. (21) for , we end up with the effective master equation of the signal mode in the adiabatic limit

(22) |

where is an injection parameter corresponding to a coherent exchange of excitations between the signal and pump modes, while accounts for signal photon pairs that are lost to the strongly damped pump mode. Equation (22) has been extensively studied in the literature Wolinsky88 (); Drummond80 (); Drummond91 (). It can be derived via standard adiabatic elimination which in the language of projection operator theory uses a time-independent projection superoperator projecting out the coherent laser field CarmichaelBook2 (). Its action on the full density matrix is given by , where is a coherent state with . The fast exponential decay of the pump correlation functions allows in this case for a Markovian approximation in the Born terms, that is .

Let us now analyze the c-MoP equations in the diabatic limit. In this case, eq. (19)/ provides us with , which introduced in eq. (16) leads to an effective master equation

(23) |

for the signal state. The pump state only enters this equation trough the amplitude which obeys eq. (21) since is traceless. Noting that this equation is equivalent to eq. (10), we conclude that the diabatic limit reduces the full c-MoP equations to the mean-field equations.

We emphasize that within these limits both eqs. (10) and (22) become exact. We have thus shown that the c-MoP theory provides us with exact equations of motion in the limits (adiabatic) and (diabatic) where it therefore becomes equivalent with well established theories CarmichaelBook2 (); FleischhauerMF (). In the remainder of the paper we will step beyond these cases in which time-scale separation is present and use the c-MoP theory to access the signal state in the scenario.

## Iv Accuracy tests and full quantum states of the signal mode

In the previous section we have shown how to deal with the non-Markovian structure of the c-MoP equations. The only remaining difficulty is given by the quartic structure of the effective equations of motion derived for the signal mode, eqs. (16) and (19). In this section we will treat the problem numerically in the Fock state basis by introducing a truncation for the Hilbertspace of the signal, where is chosen such that the results for the observables we are interested in converge up to some desired accuracy. Thus, the reduced state and the operator will be dimensional matrices. Instead of treating the pump mode in an analogous manner, we exploit the fact that the c-MoP equations of the pump mode (17) and (20) are quadratic in the bosonic operators. As a consequence we are able to describe the pump state by a set of closed equations for five variables only, the mode amplitude plus the fluctuations and (note that the first two are complex variables). At the end, we are thus effectively left with two coupled differential equations for the matrices and , with the pump equations solved either in parallel numerically or analytically as a function of signal observables.

In this section, we will compare the steady states of the classical theory from eq. (2), the steady states of the mean-field equations (10), and the steady states of the c-MoP equations (4) and (5). In order to show the accuracy of the c-MoP equations we will also determine the steady state of the full Liouville-von Neumann equation (1) in parameter regimes where it is numerically tractable. This numerical simulation is done as follows: first, we eliminate the large coherent background of the laser drive from the Liouvillian by writing , where is taken to be the classical steady-state solution of eqs. (2); then, we use the superspace formalism, where the steady-state operator and the Liouville superoperator are represented, respectively, by a vector and a matrix , and can be found as the eigenvector with zero eigenvalue of Drummond91 (); CarlosNumerics (). As the dimension of the matrix is , with denoting the pump mode’s Hilbert space dimension, this exact simulation is limited to small photon numbers.

In all the simulations we consider cases without time-scale separation between the two modes and rescale all units to the dissipation rates, i.e. we put . The only remaining parameters are the nonlinear coupling and the injection parameter .

In Fig. 2 we present results in parameter regimes where the full DOPO equation (1) can be solved numerically. In Figs. 2 and 2 we show different steady-state observables for and , respectively. It can be appreciated how the c-MoP results (blue solid line) coincide almost perfectly with the numerical results from the full master equation (red stars). The observables that we show are the pump mode’s amplitude in Figs. 2 and 2, the signal photon number in Figs. 2 and 2, and the function of the signal in Figs. 2 and 2. We also compare with the mean-field predictions of eqs. (10) (black dashed line), which in this context should be understood as the c-MoP theory up to first order, and with the classical steady-state solutions (gray thin solid line) given after eq. (2). Let us remark that despite the nonlinear nature of the mean-field and the c-MoP equations, we only find one physical solution for each of them.

All four theories agree quite well far below the critical point as the states of the signal and pump modes are close to vacuum and a coherent state induced by the external laser drive, respectively. Far above the threshold point, where the classical theory is expected to be approximately valid, we find that both the c-MoP predictions and the exact numerics agree well with the classical solutions for all observables, but with the fundamental difference that the classical theory breaks the symmetry, while c-MoP and the exact solution preserve it. The mean-field solution, on the other hand, fails to describe the state of the signal above threshold as can be appreciated from the function in Figs. 2 and 2. As expected, mean-field theory and the classical theory break down in the vicinity of the threshold point. Remarkably, this is not true for c-MoP which appears to give quasi exact results for all values of , even in cases where the interaction rate is comparable to all other system parameters.

For the experimentally relevant scenario with , the Hilbert space dimension needs to be so large that we are not able to find the numerical solution of the full master equation (1) for injection parameters close to (or above) threshold. However, we can compare the c-MoP predictions (red stars), see Fig. (4), with the perturbative approach which Drummond et al. (dark yellow dot-dashed line) developed in the vicinities of the critical point, by making a consistent multiple-scale expansion of the system’s stochastic variables within the positive representation Kinsler95 (); Chaturvedi02 (). This procedure has the virtue of being valid for any values of and , and close to threshold, concretely for , it is expected to be quasi-exact. As shown in Figs. 4 and 4, we find perfect agreement between this approach and the c-MoP theory for .

Overall, we have indeed shown the drastic impact of the Born terms, which do not only lead to a quantitative improvement as compared to the classical theory or to mean-field, but to a qualitatively different state of the signal mode. The classical theory predicts a coherent state, while the mean-field theory, i.e. the c-MoP theory up to first order, predicts a Gaussian state of the signal centered around Carlos14 (). The c-MoP theory including the born terms, hence including signal-pump correlations within a projection operator based theory, is capable of finding the full quantum state of the signal which is neither coherent nor Gaussian as shown through the function in Figs. 2 and 2.

In order to illustrate the full quantum state, we plot the Wigner function of the signal density matrix obtained from the c-MoP equations in Fig. 3 for and different values of . Let us remark that in our case in which the Wigner function is positive everywhere in the phase-space formed by the quadratures and , it can be simply interpreted as the joint probability distribution describing the statistics of measurements of these observables Braunstein05 (); Weedbrock12 (); CarlosQI (). From a computational point of view, we evaluate it from the steady-state density matrix following the method detailed in Carlos14PRL (). Far below threshold, the Wigner function shows a perfect vacuum for the signal state, see top panel of Fig. 3 for as a reference. As we cross through the critical point, two significant effects take place. First, approaching the threshold we find the well-known quadrature-noise reduction or squeezing CarmichaelBook2 (); Eberle10 (); Mehmet08 (), which is highest around the critical point Chaturvedi02 (), and reaches its asymptotic value for , with corresponding antisqueezing Carlos14 (). Second, as we cross the threshold we appreciate how the state develops two peaks centered (asymptotically) at the quadrature values predicted by the classical solution. Hence, even though the true quantum state never breaks the symmetry, it does so in two qualitatively different ways depending on whether we are below or above threshold.

In retrospect, we see that the symmetry-breaking states predicted above threshold by the standard linearization approach correspond each to one of the two distinct peaks appearing in the exact state. Far above threshold the two peaks have zero overlap and such states provide reasonable predictions for all observables which are not sensitive to symmetry breaking, that is, all observables containing even numbers of signal field operators. Of course, such a deficit can be corrected by simply using a balanced mixture of the symmetry-breaking states Carlos14 (); this construction will guide us in the next section, where we will perform a Gaussian approximation which necessarily breaks the symmetry. It is then close to the critical point where both linearization and mean-field approaches fail, whereas c-MoP provides an accurate description of the quantum state.

Let us remark that we have compared the Wigner function obtained from the c-MoP theory with the reduced signal states obtained from the full master equation, which was only possible for , and found very good agreement, the differences being completely unnoticeable to the naked eye. We emphasize again that, with the numerical solution of the c-MoP equations we are able to find the full reduced density matrices of the modes away from the adiabatic limit. This is in contrast to other approaches such as stochastic simulations Kinsler95 (); Chaturvedi02 (); Chaturvedi99 (); Pope00 () or the Keldysh formalism FleischhauerKeldysh (); Mertens93 (); Swain93 () which are naturally design to provide expectation values of the system operators.

## V Gaussian state Approximation within the c-MoP theory

Despite the fact that the complexity of solving the c-MoP equations fully numerically scales in a more favorable way than the numerical complexity of the full master equation, it still requires to integrate a number of differential equations that scales quadratically with the dimension of the truncated Hilbert space for the signal field. Therefore, it is very desirable to find an effective description of the underlying theory which is numerically more efficient and can thus cover the whole parameter space. In the remainder of this section, we implement such an idea by applying a Gaussian state approximation (GSA) consistent with the c-MoP equations (4) and (5).

Another great advantage of a Gaussian theory, apart from reaching the whole parameter space, is the efficiency in the evaluation of both steady states and dynamical quantities such as two-time correlation functions. The disadvantage of a Gaussian theory, however, is the lack of quantitative accuracy especially in the vicinity of the critical point. Nonetheless, as we show in the following, a Gaussian theory consistent with the c-MoP equations offers better quantitative accuracy than any of the previously developed Gaussian methods, particularly linearization around the classical solution or the recently-developed self-consistent linearization Carlos14 ().

The general procedure for finding a GSA for the state of a certain bosonic master equation is very simple. In a first step, we write the bosonic operators as , with , such that . In the next step we find the evolution equation for the first and second moments, which will depend on higher-order moments in general. Thus, in the final step we assume the state to be Gaussian at all times, so that all higher order moments factorize into products of first and second order moments CarmichaelBook2 (); Carlos14 (); in particular, we will encounter third order moments such as e.g. which vanish identically within the GSA, and forth order moments which factorize according to, e.g.

(24) |

After this final step, we are then left with a closed set of nonlinear equations for the amplitudes and the second order moments of the fluctuations that have to be solved self-consistently.

The standard linearization theory can be understood as a GSA on the full master equation, but with the exception that the amplitudes are not determined self-consistently, but are obtained from the classical theory. As shown by the gray thin solid line in Fig. 4 the complete suppression of quantum fluctuations when determining these amplitudes leads to unphysical results at the threshold point in the DOPO.

The self-consistent linearization method, as it is coined in Reference Carlos14 (), goes one step beyond standard linearization by consistently finding the amplitudes from the GSA still applied to the full master equation. Due to the nonlinear nature of the resulting equations of motion one can find several solutions in a given point of parameter space. However, it was shown that at the end only two types of solutions were physical, qualitatively similar to the solutions found from standard linearization, but quantitatively regularized in such a way that the unphysical results of the latter disappear. In particular, a below threshold (BT) solution was found, which does not break the symmetry, i.e. , but in contrast to the classical theory exists for all values of the injection parameter, not only for . We also found two above threshold (AT) solutions with opposite phase which break the symmetry, that is, , but appear only above a certain injection parameter which is larger than the classical threshold value. Interestingly, we point out that the BT solution found through this self-consistent linearization is exactly equivalent to the mean-field theory introduced in Section III.1.

Motivated by these findings, we apply a GSA to the c-MoP equations. Concretely, we calculate all first and second order moments of the pump and signal fluctuations from the c-MoP eqs. (16), (17), (19), and (20), and apply the factorization of higher order moments as explained above. In strong contrast to the GSA on the full master equation, we do not need to assume a Gaussian state for the full state but only for the reduced state of the signal . Hence, we expect similar qualitative results but with a higher quantitative accuracy.

Indeed, this is what we find and illustrate in Fig. 4 for and . We plot steady state expectation values for the pump amplitude in Fig. 4, the signal photon number in Fig. 4, and the function of the signal in Fig. 4, all as a function of the injection parameter . The blue solid line in Fig. 4 shows the below threshold solution of the GSA on the c-MoP equation, while the green dashed line illustrates the above threshold solution. The latter fulfills and is therefore more likely to provide physically consistent results than the BT solution whenever they coexist. In Fig. 4 we show how the AT solution indeed gives the correct value for the function, what indicates that each of the AT solutions corresponds to one of the peaks of the Wigner function, see Fig. 3, and considers Gaussian fluctuations around it. In order to illustrate this point even further, we show in Fig. 3 the Wigner function Braunstein05 (); Weedbrock12 (); CarlosQI () corresponding to the GSA on the c-MoP equations (as explained in the previous section, above threshold we take the balanced mixture of the two symmetry breaking solutions, such that the resulting state preserves the symmetry).