Higher-order interactions in quantum optomechanics: Analytical solution of nonlinearity

A method is described to solve the nonlinear Langevin equations arising from quadratic interactions in quantum mechanics. While, the zeroth order linearization approximation to the operators is normally used, here first and second order truncation perturbation schemes are proposed. These schemes employ higher-order system operators, and then approximate number operators with their corresponding mean boson numbers, only where needed. Spectral densities of higher-order operators are derived, and an expression for the second-order correlation function at zero time-delay has been found, which reveals that the cavity photon occupation of an ideal laser at threshold reaches $\sqrt{6}-2$, in good agreement with extensive numerical calculations. As further applications, analysis of the quantum anharmonic oscillator, calculation of $Q-$functions, analysis of quantum limited amplifiers, and nondemoliton measurements.

Usually, the analysis of quantum optomechanics is done within the linearized approximation of photon ladder operators, normally done asâ →ā+δâ with |ā| 2 =n being the mean cavity photon number, while nonlinear terms in δâ are ignored. But this suffers from limited accuracy wherever the basic optomechanical interaction H OM = g 0n (b +b † ) is either vanishingly small or non-existent. In fact, the single-photon interaction rate g 0 can be identically made zero by appropriate design [33][34][35][36], when quadratic or even quartic effects are primarily pursued. This urges need for accurate knowledge of higher-order interaction terms.
Some other optomechanical phenomena such as four-wave mixing, also can be suitably understood by incorporation of higher-order interaction terms [37]. Recent experiments [38,39] have already established the significance and prominent role of such type of nonlinear interactions. In fact, quadratic nonlinear optomechanics [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55] is now a well recognized subject of study even down to the single-photon level [56], for which circuit analogues have been constructed [57,58] and may be regarded as fairly convenient simulators [59][60][61] of much more complicated experimental optomechanical analogues. Dual formalisms of quadratic optomechanics are also found in ultracold atom traps [62,63] as well as optical levitation [64]. Such types of nonlinear interactions also appear elsewhere in anharmonic quantum circuits [65]. Quadratic interactions are in particular important for energy and non-demolition measurements of mechanical states [1,2,4,[66][67][68]. While the simple linearization of operators could be still good enough to explain some of the observations, there remains a need for an exact and relatively simple mathematical treatment. Method of Langevin equations also normally fails, and other known methods such as expansion unto number states and master equation, require lots of computation while giving little insight to the problem.
Perturbative expansions and higher-order operators have been used by other researchers to study noise spectra of lasers [69][70][71][72]. Also, the master equation approach [73,74] can be used in combination with the quasi-probablity Wigner functions [75,76] to yield integrable classical Langevin equations. Nevertheless, a method recently has been proposed [77], which offers a truncation correlation scheme for solution of driven-dissipative multi-mode systems. While being general, it deals with the time evolution of expectation values instead of operators within the truncation accuracy, so the corresponding Langevin equations cannot be analytically integrated.
Alternatively, a first-order perturbation has been proposed to tackle the nonlinear quadratic optomechanics [78]. This method perturbatively expands the unknown parameters of classical Langevin equations for the nonlinear system, and proceeds to the truncation at first order. However, the expansion is accurate only where the ratio of photon loss rate to mechanical frequency κ/Ω is large. This condition is strongly violated for instance in superconductive electromechanical systems.
So far and to the best knowledge of author, no treatment of quadratic interactions using Langevin equations for operators has been reported. This paper presents a perturbative mathematical treatment within the first and second order approximations to the nonlinear system of Langevin equations, which ultimately result in an integrable system of quantum mechanical operators. The trick here is to introduce operators of higher dimensionality into the solution space of the problem. Having their commutators calculated, it would be possible to set up an extended system of Langevin equations which could be conveniently solved by truncation at the desirable order. To understand how it works, one may consider the infamous first order quadratic nonlinear Riccati differential equation [79,80], which is exactly integrable if appropriately transformed as a system of two coupled linear first order differential equations. Alternatively, Riccati equation could be exactly transformed into a linear second order differential equation, too. But this is not what we consider here, since it will result in a much more complicated second-order system of Langevin equations involving derivatives of noise terms.

A. Hamiltonian
A nonlinear quadratic optomechanical interaction in the most general form [108] is here defined as which by defining the photon number operatorn =â †â takes essentially the same algebraic form. Direct expansion of (1) shows that it essentially brings in a different interaction type compared to (2). Doing so, we obtain H = γ(b 2 +b †2 ± 2m ± 1) 2 (â 2 +â †2 ± 2n ± 1) wherem =b †b . Hence, (1) includes interactions of typê a 2b2 ,â 2b †2 , and so on, which are absent in (2). It should be noticed that the widely used standard optomechanical interaction H OM results in nonlinear and linear Langevin equations when expressed respectively in the terms of {â,b} and {n,x}. Hence, this type of interaction is not addressed here. In addition to the above Hamiltonians (1,2), there exist still other types of nonlinear optomechanical interactions [16,109] such as H = g(b ±b † )(â 2 ±â †2 ), which is also not considered explicitly here, but can be well treated using the scheme presented in this article.

B. Linear Perturbation
This approach is being mostly used by authors to solve the systems based on either (1) or (2). To this end, ladder field operators are replaced with their perturbations, while product terms beyond are neglected and truncated. Obviously, this will give rise to interactions of the type (b ±b † ) 2 (qδâ + q * δâ † ), where q = 2γ(ā ±ā * ) for (1) and q = γā for (2) is some complex constant in general, and δâ now represents the perturbation term around the steady state average |ā| = √n . This technique is mostly being referred to as the linearization of operators, and directly leads to an integrable set of Langevin equations if also applied to the mechanical displacement as well.

C. Square Field Operators
Here, we define the square field operators [108] [ĉ,n] = 2ĉ, Defining the phonon number operator asm =b †b , in a similar manner we could write The set of commutator equations (4) and (5) enables us to treat the quadratic nonlinear interaction perturbatively to the desirable accuracy, as is described in the following.

D. Langevin Equations
The input/output formalism [29][30][31][32] can be used to assign decay channels to each of the quantum variables of the system. This will result in the set of Langevin equations Here, [Γ] is supposed to be diagonal for simplicity. From the scattering matrix formalism we also have Hence, taking w as the angular frequency and performing a Fourier transform on (6), the scattering matrix is found by using (7) and (8) wherex is any system operator, which is here taken to be the same asẑ to comply with (8). By setting eitherẑ =ĉ orẑ =d the commutators in (10) by (4) or (5) always lead back to the same linear combination of these forms. Thus, the new set of Langevin equations is actually linear in terms of the square or higher-order operators, if perturbatively truncated at a finite order. So, instead of solving the nonlinear system in linearized 2 × 2 space {A} T = {â,b}, one may employ an expanded dimensional space with increased accuracy. There, truncation and sometimes mean field approximations are necessary to restrict the dimension, since commutators of new operators mostly lead to even higher-orders and are thus not closed under commutation. As examples, a 4 × 4 space {A} T = {â,d,d † ,m} truncated at the first-order, or a 6 × 6 space {A} T = {ĉ,ĉ † ,n,d,d † ,m} truncated at the second-order could be used for (1,2). To illustrate the application of this method, we describe two examples in the next section. It could be extended to the accuracy of the second-order perturbation too, by defining appropriate cross product operator terms between photonic and phononic partitions.

III. EXAMPLES
Here, we describe two examples from the nonlinear interactions of having type (1) or (2).

A. Standard Quadratic Interaction (2)
Analysis of such systems requires analysis in a 4-dimensional space, spanned by {A} T = {â,d,d † ,m}. Taking the plus sign here without loss of generality and after dropping a trivial non-interacting term H 0 = γn, the nonlinear interaction is This can be found by expansion of (2), plugging in (3) and [b,b † ] = 1, and dropping a trivial term γn. Using (5), [â,n] =â and [â † ,n] = −â † in the non-rotating frame of operators, and ignoring the self-energy Hamiltonian H self = (ω + γ)n + Ωm for the moment, Langevin equations becomė So far, the set of equations (12) is exact. However, integration of (12) is still not possible at this stage, and taking Fourier transformation must be done later when arriving at a linear operator system. We present a first-order and second-order perturbative method to deal with this difficulty.
It should be furthermore noticed that using a non-rotating frame with the self-energy Hamiltonian H self not ignored, would have resulted in identical equations, except with the addition of the trivial terms −i∆â, −i2Ωd, and +i2Ωd † respectively to the first three equations, where ∆ = ω + γ − ν is the optical detuning with ν being the cavity optical resonance frequency, and ω and Ω are respectively the optical and mechanical frequencies. Also, the damping coefficient in high mechanical quality factor Q m limit could be estimated as Γ 2 = 2Γ m , where Γ m is the damping rate of theb phononic field. Here, it is preferable not to use the rotating frames since the coefficients matrix [M] becomes time-dependent.
1. First-order Perturbation to (12) Now, if the photon and phonon baths each have a mean boson number respectively as n =n and m =m, we could immediately write down the linear system of equations in the non-rotating frame of operators and neglection of self-energies H self asȧ which is now exactly integrable. Here, we use the linearization 2âd = (ā + δâ)d +â(d + δd) →ād +dâ, whered = 1 2ā 2 and higher-order terms of the form δâδd are dropped, and so on. But this cannot be applied tonm =â †âm sincen andâ † are absent from the basis. Furthermore, any linearization of this expansion would generate termsâm andâ †m which are still nonlinear. Both of these issues can be resolved by a second-order perturbation as follows next. This results in the operator equations where The set of equations (14) is linear and can be easily addressed by standard methods of stochastic Langevin equations used in optomechanics [1,2,4,29,30] and elsewhere. More specifically, one may employ analytical Fourier methods in frequency domain as an matrix algebraic problem to obtain spectra of variables, or integrate the system numerically by stochastic numerical methods in time domain to obtain time dependent behavior of expectation values.
All that remains is to find the average cavity boson numbers for photonsn and phononsm. In order to do this, one may first arbitrate d/dt = 0 in (13) at steady state, and then use the equality of real parts in first equation to find the expression forn. Doing this, results inn = 4|ā in | 2 /Γ 1 where |ā in | represents the amplitude of coherent laser input. Also, the initial cavity phonon occupation number at t = 0 could be estimated simply asm = 1/ [exp( Ω/k B T ) − 1] [29,30], where k B T is the thermal energy with k B and T being respectively the Boltzmann's constant and absolute temperature. Detailed numerical examinations reveal that the system of equations (14) is generally very well stable with {eig[M]} < 0 at sufficiently low optical intensities.

B. Full Quadratic Interaction (1)
Analysis of a fully quadratic system requires analysis in a 6 × 6 dimensional space, spanned by {A} T = {ĉ,ĉ † ,n,d,d † ,m}. Taking both of the plus signs here, the Hamiltonian could be written as where a trivial non-interacting term H 0 = 2 γ(1 +n +m +d +ĉ +d † +ĉ † ) is dropped. The set of Langevin equations can be obtained in a similar manner, and in non-rotating frame of operators with neglection of self-energies H self = (ω + 2γ)n + (Ω + 2γ)m for the moment, results iṅ Similar to (12), the damping rate for sufficiently high optical quality factors Q could be estimated as Γ 1 = 2κ, where κ is the damping rate of theâ photonic field. Quite clearly, should we have not ignored the self-energy Hamiltonian H self , then addition of the diagonal terms −i2∆ĉ, +i2∆ĉ † to the first two where ∆ = ω + 2γ − ν with ν being the optical cavity resonance frequency, and similarly −i2Ωd and +i2Ωd † to the fourth and fifth equations would have been necessary. These are not shown here only for the sake of convenience. Again, it is emphasized that transformation to the rotating frame of operators here would make the coefficients time-dependent in an oscillating manner, and it is far better to be avoided for these classes of nonlinear problems.
1. First-order Perturbation to (16) In a similar manner to (13), we may assume photon and phonon baths each have a mean boson number respectively as n =n and m =m, which giveṡ We here need to assume the redefinition √ ∆ 1 = (n + 1 2 ) √ Γ 1 . Now, without taking H self into account, this will lead to the linear system of matrix Langevin equations which is, of course, integrable now. The initial cavity boson numbersn andm can be set in the same manner which was done for the system of equations (14). Numerical tests reveal that (18) is conditionally stable if the optical intensity is kept below a certain limit on the red detuning, and is otherwise unstable.
C. Second order Perturbation to (12,16) The set of Langevin equations (12,16) can be integrated with much more accuracy, if we first identify and sort out the cross terms as individual operators. For instance, (16) contains the cross operatorsĉd,ĉd † ,ĉm,ĉ †d ,ĉ †d † ,ĉ †m , nd,nd † , as well asnm which is self-adjoint. These constitute an extra set of nine cross operators to be included in the treatment. All these cross operators are formed by multiplication of photonic and phononic single operators, whose notation order, such asĉd =dĉ and so on, is obviously immaterial. Now, one may proceed first to determine the commutators between these terms where relevant, which always result in linear combinations of the other existing terms. This will clearly enable a more accurate formulation of (16) but in a 6 + 9 = 15 dimensional space, which is given by the array of operators The independent non-trivial quadratic commutator equations among cross operators here are found after tedious but straightforward algebra as [ĉd,nm] = (n +m + 4)ĉd, The rest of commutators among cross operators are either adjoints of the above, or have a common term which makes their evaluation possible using either (4) or (5). Commutators among cross operators and single operators can be always factored, such as [ĉd,n] = [ĉ,n]d. Commutators among single operators are already known (4,5). It can be therefore seen that commutators (19) always lead to operators of higher orders yet, so that they do not terminate at any finite order of interest by merely expansion of operators basis. This fact puts the perturbative method put into work. There are, however, nonlinear systems such as semiconductor optical cavities [71,74] in which higher-order operators yield an exact closed algebra and satisfy a closedness property within the original space by appropriate definition.
The set of ten commutators now can be perturbatively linearized as a second-order approximation, by replacing the number operators with their mean values, wherever needed to reduce the set of operators back to the available 15 dimensional space. This will give rise to the similar set of equations after some algebra where the reduction of triple operator products among single and cross operators as 4xŷẑ →xŷẑ +xȳẑ +ȳzx +zxŷ is used where appropriate. For instance, the term 4nm 2 is replaced asmmn + 2mnm +m 2n and so on. Also, similar to (13), products among single operators are reduced as 2xŷ →xŷ +ȳx. This is somewhat comparable to the mean field approach in cross Kerr optomechanics [110].
There are two basic reasons why we have adopted this particular approach to the linearization and cuting off the diverging operators of higher orders. The first reason is that number operators vary slowly in time as opposed to their bosonic counterparts which oscillate rapidly in time, given the fact that the use of rotating frames is disallowed here. Secondly, number operators are both positive-definite and self-adjoint, and thus can be approximated by a positive real number. These properties makes the replacementsn →n andm →m reasonable approximations, and the replacement with mean values needs only to be restricted to the number operators, to yield a closed algebra necessary for construction of Langevin equations. Hence, the correct application of replacements only to the triple operator products appearing in (19) will make sure that no operator having an order beyond than that of cross operators will appear in the formulation.
Anyhow, it can be seen now that all approximate commutators in (20) allow the set of operators {A} T ∪ {1} = {1,ĉ,ĉ † ,n,d,d † ,m,ĉd,ĉd † ,ĉm,ĉ †d ,ĉ †d † ,ĉ †m ,nd,nd † ,nm} to take on linear combinations of its members among every pair of commutations possible, where1 is the identity operator. Obviously, this approximate closedness property now makes the full construction of Langevin equations for the operators belonging to {A} possible. It is noted that1 is not an identity element for the commutation.
We Having therefore these ten commutators (20) known, we may proceed now to composing the second-order approximation to the nonlinear Langevin equations (16), from which a much more accurate solution could be obtained. Here, the corresponding Langevin equations may be constructed at each step by setting bothẑ andx in (10) equal to either of the 15 operators, while the noise input terms for cross operators is a simple product of related individual noise terms. The linear damping rates of higher-order operators is furthermore simply the sum of individual damping rates of corresponding single operators, which completes the needed parameter set of Langevin equations.

A. Optomechanical Interaction & Drive Terms
The method described in the above can be simultaneously used if other terms such as the standard optomechanical interaction H OM is non-zero, or there exists a coherent pumping drive term which can be expressed as H d = k F kb † + F * kb , where F k are time-dependent drive amplitudes. While H d does not appear directly in the Langevin equations, treatment of H OM requires inclusion of additional Langevin equations forâ andb where appropriate, as well as few extra terms in the rest. This can be done in a pretty standard way, and is not repeated here for the sake of brevity [1,2,[29][30][31][32].

B. Multi-mode Fields
The analysis is also essentially unaltered if there are more than one mechanical mode to be considered [19,111,112], and the method is still easily applicable with no fundamental change. Suppose that there are a total of M mechanical modes with the corresponding bosonic operatorsb k andb † k where k ∈ [1, M ]. Then, these modes are mutually independent in the sense that [b j ,b k ] = 0 and [b j ,b † k ] = δ jk . The set of commutators (5) will be usable for all M modes individually and as a result (19) and therefore (20) may be still used. The first and second order perturbations will respectively result in 3 + 3M = 3(M + 1) and 3 + 3M + 9M = 3(4M + 1) equations. The redefined set of operators will be respectively now Similarly, in case of N optical modes satisfying [â j ,â k ] = 0 and [â j ,â † k ] = δ jk , the set of commutators (4) can be used and the operator set should be now expanded as M ]} respectively for first and second order perturbations. Hence, the corresponding dimensions will be now respectively either 3(N +M ) or 3(N + M + 3N M ). Higher-order commutators (19) and (20) can be still used again by only addition of appropriate photonic j and phononic k mode indices to the respective operators contained in the expanded operator basis set {A}.

C. Noise Spectra
The required noise spectra [113] of cross operators is clearly a product of each of the individual terms, since the nature of particles are different. However, the noise spectra of quadratic operators themselves need to be appropriately expressed. For instance,d in actually corresponds to the spectral input noise of the square operatord =bb/2 √ Γ from (3), which clearly satisfiesd in (t) = 1 2b in (t)b in (t)/ √ Γ, ord in (w) = 1 2b in (w) * b in (w)/ √ Γ in the frequency domain, where * merely represents the convolution operation. Therefore, onceâ in (w) andb in (w) are known, all relevant remaining input noise spectra could be obtained accordingly using simple convolutions or products in frequency domain.
As a result, the corresponding spectral density of the noise input terms to the cross operators can be determined from the relevant vacuum noise fluctuations and performing a Fourier transform. For instance, we have S CDCD . Then Isserlis-Wick theorem [38,114] could be exploited to yield the desired expressions. If we assume where the dimensionless correlation integrator runs on phase, instead of time, as then the functions ζ(·), ψ(·), and the operatorv(·) should be all having the dimension off 2 (·) as well. That means if f is dimensionless, which is the case for the choice of ladder operators, then ζ(·), ψ(·), andv(·) become dimensionless, too. The functions ζ(·) and ψ(·) together can cause squeezing or thermal states if appropriately defined [29,32]. By Isserlis-Wick theorem applied to scalars we have gives for the operators Hence, for a given stochastic process where f (t)f (τ ) = 0, f (t)f † (τ ) = Ψ(t−τ ), and having the scalar commutator Now, suppose that we have a coherent field of photons at the angular frequency ω with an initial Gaussian distribution, in which Ψ(t) = exp(−χ 2 ω 2 t 2 /2) exp(−iωt) and Υ(t) = Ψ(t), while having the linewidth ∆f = 1 2π χω. Clearly, χ is a dimensionless and positive real number. In the limit of χ → 0 + , the expected relationship Ψ(t) = √ 2πδ(ωt)/χ is easily recovered.
This particular definition of the correlating function Ψ(t) ensures that the corresponding spectral density is appropriately normalized, that is = 1.
Hence, one may obtain the following spectral density which is centered at the doubled frequency 2ω, has a linewidth of √ 2∆f , and satisfies the property Once the spectral densities of input noise terms are found, spectral densities of all output fields immediately follows (8,9)  Many of the important features of an interacting quantum system is given by its second-order correlation function g (2) (0) at zero time-delay [98-100] defined as It is fairly easy to estimate this function once the spectral densities of all higher order operators of the nonlinear system are calculated. For this purpose, we may first employ the definition (3) to rewrite Estimation of the average within brackets can be done by having S CC [w] = 1 4 S A 2 A 2 [w] corresponding to the higherorder operatorĉ. This can be assumed to has been already found from knowledge of the scattering matrix [S(w)], spectral densities of input fields {A[w]} in , and subsequent derivation of spectral density array of output fields {A[w]} out . Then, S CC [w] will be simply an element of the vector {A[w]} out . Using (24), this results in a fairly brief representation With the assumptions above for an ideal initial Gaussian distribution, we have Ψ(0) = Υ(0) = 1 and thus g (2) (0) = 4( 1 2 −n)/n 2 . One should have in mind that this relationship cannot be readily used for a coherent radiation, since for a practical laser the true statistics is Poissonian and not Gaussian. This analysis thus reveals that the cavity occupation number of such an ideal laser with the threshold defined as g (2) (0) = 1 is exactlyn = √ 6 − 2 ≈ 0.450. This is in contrast to the widely used assumption of quantum threshold conditionn = 1 [115][116][117][118][119][120]. Interestingly, a new study [121] of photon statistics in weakly nonlinear optical cavities based on extensive density matrix calculations [122,123] yields the valuen = 0.4172, which is in reasonable agreement to our estimate. An earlier investigation on quantum-dot photonic crystal cavity lasers [124,125] also gives the valuen = 0.485.

V. ANHARMONIC OSCILLATOR
The quantum anharmonic oscillator appears in many nonlinear systems including quadratic optomechanics [126,127], where our method here is applicable. The anharmonic Kerr Hamiltonian is [128,129] in which ζ is a constant. It is well known that in case of ζ > 2ω this system exhibits an effective bistable potential, and is otherwise monostable. However, we are here much interested in a slightly different but more complicated form given by [130] which is monostable or bistable if both ω and ζ are respectively positive or negative. This type of nonlinearity is of particular importance in fourth-order analysis of qubits [131][132][133][134][135][136][137]. While the Hamiltonian (32) is for a singlemode field, the case of multi-mode electromagnetic field could be easily devised following the existing interaction Hamiltonians [130] and the presented method in this article. Nevertheless, the above expression after some algebraic manipulations can be put into the form where a trivial constant term ζ is dropped. Here, we may proceed with the 8-dimensional basis operator set {A} T = {ĉ,ĉ † ,n,n 2 ,ĉ 2 ,ĉ †2 ,nĉ,ĉ †n }, resulting in a second order perturbation accuracy. Treating this problem using the Langevin equation (10), regardless of the values of ζ and ω, is possible, only if the following non-trivial exact commutators [nĉ,ĉ †n ] = 1 2 4n 2 − 3n + 2 n, are known, which may be found after significant algebra. The rest of required commutators which are not conjugates of those in the above, can either directly or after factorization of a common term be easily found from (4). Again, the set of commutators (34) does not yet satisfy the closedness property within {S} = span({A} ∪ {1}), unless the approximate linearization [ĉ 2 ,ĉ †n ] = 3 (n + 2)nĉ + 6ĉ, is employed. The rest of the process is identical to the one described under (20). Construction of the respective noise terms is also possible by iterated use of the results in §IV C and so on.

A. The Husimi-Kano Q-functions
It is mostly appropriate that moments of operators are known, which are scalar functions and much easier to work with. The particular choice of Q−functions [138] is preferred when dealing with ladder operators, and are obtained by taking the expectation value of density operator with respect to a complex coherent state |α and dividing by π. This definition leads to a non-negative real valued function Q(α) = Q( [α], [α]) of |α . Then, obtaining Q−function moments of any expression containing the ladder operators would be straightforward [138]. However, it must be antinormally ordered, with creators be moved to the right. In {A} T above all operators are actually in the normal form, exceptn 2 . It is possible to put the nontrivial members of {A} in the antinormal order While evaluating Q−function moments,â andâ † are replaced with α and α * respectively as All remains now is to redefine the array of Q−functions bases, using common terms as { A } T = {α 2 , α * 2 , |α| 2 , |α| 4 , α 4 , α * 4 , α 2 |α| 2 , α * 2 |α| 2 } from which the original Q−functions could be readily restored. This translates into a set of scalar differential equations which conveniently could be solved. Fluctuations of noise terms also vanish while taking the expectation values, and only their average values survive. To illustrate this, suppose that the system is driven by a coherent fieldâ in with the normalized electric field amplitude β = α/ √ 2 and at the frequency ω. Then, the Q−function moments of the input fields after defining the loss rates Γ 3 = 2Γ 2 = 4Γ 1 become â in = √ 2Γ 1 β, ĉ in = √ Γ 2 β, n in = √ Γ 2 (2|β| 2 + 1), ĉ 2 in = √ Γ 3 β 2 , and n inĉin = √ Γ 3 β 2 (2|β| 2 + 3).

C. Quantum Nondemolition Measurements
Quantum nondemolition measurements of states require a cross-Kerr nonlinear interaction of the type H = ωâ †â + Ωb †b + χâ †âb †b = ωn + Ωm + χnm, in whichâ andb fields respectively correspond to the probe and signal [105,106]. This system can be conveniently analyzed by the preferred choice [105] of the higher-order operators {A} T = {n,m,Ĉ,Ŝ}, whereĈ  (6) is time-dependent [104,139]. For instance, the ultimate optomechanical cooling limit is a function of system dynamics [140]. Then, integration should be done numerically, since exact analytical solutions without infinite perturbations exist only for very restricted cases. This is, however, beyond the scope of the current study.

VI. CONCLUSIONS
A new method was described to solve quadratic quantum interactions using perturbative truncation schemes, by including higher-order operators in the solution space. Spectral densities of higher-order operators, calculation of the second-order correlation function, as well as the quantum anharmonic oscillator and transformation to scalar forms using Q−functions were discussed. Finally, applications of the presented approach to quantum limited amplifiers, and nondemolition measurements were demonstrated. Amir H. Ghadimi for comments and/or discussions. The author is highly indebted to Dr. Hiwa Mahmoudi at Institute of Electrodynamics, Microwave and Circuit Engineering in Technische Universität Wien, and in particular, the Laboratory for Quantum Foundations and Quantum Information on the Nano-and Microscale in Vienna Center for Quantum Science and Technology (VCQ) at Universität Wien for their warm and receptive hospitality during which the numerical computations and final revisions took place. The huge effort needed in improving the presentation of this article has not been possible without support and encouragement of the celebrated artist, Anastasia Huppmann.

APPENDICES
We show by numerical solution of a stochastic nonlinear operator differential equation, and also taking its expectation values within the Mean Field Approximation, that the proposed analytical scheme in the above referenced manuscript works very well, is convergent, and uniformly converges to the accurate solution.
We also demonstrate the second-order exact solution to the standard optomechanical Hamiltonian, and show that the proposed method can reproduce the side-band asymmetry in quantum optomechanics surprisingly well, which is a rigorous proof of the quantum mechanical capability of the proposed approach.

Appendix A: Nonlinear First-order Circuit
First, consider the infinitely-ordered nonlinear operator equation which models the voltage operator of an RC circuit shunted by a nonlinear ideal diode, driven by a sinusoidal voltage source v (t) = V 0 e −αt sin (ωt), and stochastic noisen (t). We suppose that the noisen (t) is governed by a Weiner process. Here, and without loss of generality, both κ and µ are taken to be positive real parameters. Hence, this model does not include an oscillating part due to an imaginary µ, which could have been otherwise absorbed intô u (t) by a rotating frame transformation. This particular choice also eliminates the imaginary part ofû (t). We also assume here, for the illustrative purpose of this example, that µ = V 0 = τ = 1, ω = 2πα, and ω = 2π × 1kHz. The reason for choosing this particular differential equation is that it is nonlinear to the infinite order, and also the extended basis of higher order operators all commute and therefore trivially form a closed basis.
Using the proposed method in the paper under consideration, this above operator equation can be first put into the infinitely-ordered linear system of ordinary differential equations as . . .
Subsequently, the input terms can be linearized using the proposed method in the paper. Doing this results in . . .
is the time-average of the input. Now, the above system of equations can be exactly integrated, after truncation to a finite-order.
Using an extensive code written in Mathematica, the above system of linear stochastic equations can be treated and integrated as an Itô process, and the results for various orders of truncation between 2 and 6 versus the numerically exact solution are displayed in Fig. 1.
It is still not quite clear that the method is convergent to the exact solution, since the Itô integration of a Weiner process every time is carried over a different sequence of random numbers. This difficulty cannot be avoided in principle, since there is no way to reset the numerically random sequence.
Therefore, as a double check, we take the expectation values, which discards the noise term, and transform a meanfield approximation to reach a similar system of differential equations, however, expressed in terms of the expectation value function û (t) and its higher orders. This is equivalent to solving the nonlinear differential equation given the fact that n (t) = 0. Doing this immediately reveals the convergence property of our proposed method, illustrated in Fig. 2. As it can be clearly verified, the numerical solutions are so rapidly and accurately converging to the exact solution, that they are practically indistinguishable beyond the two lowest truncation orders.

Its standard optomechanical Hamiltonian reads
In order to form a closed basis of operators, we may choose which forms a 6 × 6 system of Langevin equations. It is easy to verify that this system is exactly closed, by calculation of all possible commutation pairs between the elements. Out of the 6! commutators, the non-zero ones are âb,n =âb, âb † ,n =âb † , âb,âb † = [ĉ,n] = 2ĉ, This is obviously a closed basis. Now, one may proceed with composition of the Langevin equations. They are given by d dt where γ = κ + Γ, we have setx =â in (10) in all equations,n in =â †â in +ââ † in andĉ in = 2ââ in . This subsequently can be linearized to get the second-order accurate optomechanical system of equations as d dt in which F = g 0 √n , K = 4nκ, G = g 0n , L = G + g 0 (m + 1), and can be further approximated by L = G under normal experimental conditions of an ultracold cavity. The average population valuem = 1/[exp ( Ω/k B T )], whilen can be obtained from the steady state solution of the first row by replacements of input noise term √ κâ in → α, where α is the input photon flux.
Doing this gives the third order equation in terms of √n as Solution of this third-order algebraic equation gives the three distinct solutions is a complex number. These three roots clearly become two distinct solutions in the limit of lossless cavity with γ = 0 as long as 27g 0 C 2 > 4B 2 r , and very well reproduce the expected bistable behavior of cavity photon numbern with respect to various parameters such as detuning ∆ and input photon flux α. This has been shown for typical normalized variables in Figs. 3 and 4.
For the noise terms, we may insert Here, the spectral density ofâ 2 in has already been calculated in Section IV C. As opposed to the second-order accurate optomechanical Langevin equations, the first-order accurate equations are simply obtained by truncating (B4) which recovers the well-known 3×3 system with the obviously closed basis A T = â,b,b † after simple algebraic manipulation d dt Here, we have used the further replacements g 0n → G and g 0âb → Fb. However, if we had made the replacement g 0n → Fâ we would get d dt Finally, we are able to show that this proposed approach can reproduce the side-band asymmetry of an optomechanical cavity, which is a fundamental quantum mechanical property.
In the context of quantum optomechanics, it is general practice to solve either (B11) or (B12) with asymmetric noise spectral densities such as S BB (+ω) = S BB (−ω). While being a purely nonlinear quantum mechanical effect, this is the only known way at this moment to reproduce the side-band asymmetry. Here, we can show indeed that asymmetric spectral noise terms are not needed to obtain the much expected side-band asymmetry. This result becomes possible by considering that in the system (B5) both of the termsâb andâb † are present in the second-order accurate basis. The first term and its conjugate {âb,â †b † } represent the parametric interactions on the blue side, while the second term and its conjugate {âb † ,âb † } represent the hopping interactions on the red side. By assumption of some test values to the input parameters as well as symmetric input noise terms such as S BB (+ω) = S BB (−ω), the spectral density of the collected output light from the cavity can be calculated.
This has been done for the two cases of pumping on the red-side and blue-side shown in Fig. 5. More interestingly and as the ultimate verification of the approach, the pumping can be tried exactly on the cavity resonance. Carrying out this procedure allows calculation of cavity sidebands at blue and red, the ratio of which is here plotted and shown in Fig. 6. This figure reveals that there actually exists an easily observable asymmetry, which is hallmark of nonlinear quantum mechanical effects in cavity optomechanics. This should have rigorously proven, while leaving no further doubt, that the second-order accurate system of equations (B5) can reproduce the side-band asymmetry in a satisfactory manner.
Detailed numerical simulation of standard and non-standard optomechanical problems using the proposed approach in paper, and verification against experiments remains as a subject of a future study. Side-band asymmetry, defined as the ratio of spectral densities at opposite frequencies, in an optomechanical cavity obtained from numerical solution of (B5) with symmetric input noise terms. Horizontal axis is the normalized detuning with respect to the mechanical frequency.