Higher-Order Interactions in Quantum Optomechanics: Analysis of Quadratic Terms

This article presents a full operator analytical method for studying the quadratic nonlinear interactions in quantum optomechanics. The method is based on the application of higher-order operators, using a six-dimensional basis of second order operators which constitute an exactly closed commutators. We consider both types of standard position-field and the recently predicted non-standard momentum-field quadratic interactions, which is significant when the ratio of mechanical frequency to optical frequency is not negligible. This unexplored regime of large mechanical frequency can be investigated in few platforms including the superconducting electromechanics and simulating quantum cavity electrodynamic circuits. It has been shown that the existence of non-standard quadratic interaction could be observable under appropriate conditions.

Many of the studies in this field utilize linearization of photonic â and phononic b ladder operators around their mean values as → +â a a and → +b b b, where the substituted ladder operators now represent field fluctuations around their respective mean values. This way of linearization is however insufficient for quadratic [18][19][20][21][22][23][24][25] and higher-order interactions where the resulting Langevin equations [26][27][28][29] are expected to be strongly nonlinear. Full linearization essentially transforms back every such nonlinear interaction picture into the simple linearized form of + +ˆˆ † † a a b b ( )( ) . As a result, the behavior of systems under study due to quadratic and higher-order interactions becomes indifferent to that of ordinary linearized optomechanics, implying that in this way of linearization some of the important underlying nonlinear physics may be lost.
Recently, theory of optomechanics has been revisited by the author 30 as well as others for nonlinear 31 and quadratic 32 interactions. It has been shown that a non-standard quadratic term could exist due to momentum-field interaction and relativistic effects, the strength of which is proportional to ω Ω ( / ) 2 with Ω and ω respectively being the mechanical and optical frequencies 30 . Normally, such momentum-field interactions are not expected to survive under the regular operating conditions of large optical frequencies ω Ω  . However, it would be a matter of question that whether they could be observable in spectral response of a cavity, given that the ω and Ω could be put within the same order of magnitude? The answer is Yes.
Superconducting electromechanics provides a convenient means to observe quadratic effects at such conditions. Not only high photon cavity numbers could be attained, but also mechanical and radio or microwave frequencies could be set within the same order of magnitude with relative ease. This conditions have been actually met at least in one reported experiment 33 , where the mechanical and superconducting circuit frequencies are designed to be equal at ω π = Ω = × 2 720kHz. Tuning to the lowest-order odd-profiled mechanical mode, or using the membrane-in-the-middle setup 8,33 in optomechanical experiments could altogether eliminate the standard optomechanical interaction, leaving only the quadratic terms and higher. Apart from experimental considerations for observation of quadratic effects, there remains a major obstacle in theoretical analysis of combined standard and non-standard quadratic interactions. To the best knownledge of the author, this regime has been investigated theoretically for the case of only standard quadratic interaction using the time-evolution operators 25 .
Growing out of the context of quantum optomechanics, the method of higher-order operators developed by the author 34 can address problems with any general combination of nonlinearity, stochastic input, operator quantities, and spectral estimation. Higher-order operators have been already used in analysis of nonlinear standard optomechanics 35 , where its application has uncovered effects known as sideband inequivalence, quantities such as coherent phonon population, as well as corrections to the optomechanical spring effect, zero-point field optomechanical interactions, and a minimal basis with the highest-order which allows exact integration of optomechanical Hamiltonian subject to multiplicative noise input. Also, it has been independently used 36 for investigation of quadratic effects. But this method has not been verified yet for non-standard quadratic optomechanics, which is the central topic of this study.
In this article, we employ a six-dimensional basis of higher-order operators, all being second order, which satisfy an exact closed commutation relations. This basis can be used to analyze the quadratic interactions of both standard and non-standard types, which has been so far not done. It has been shown that the momentum-field interaction, referred to as the non-standard quadratic term, does have observable effects on the spectral response of the optomechanical cavity, if the design criteria could violate ω Ω  . Hence, the effect of this interaction should not be overlooked when the ratio ω Ω/ is non-negligible. The non-standard term appears to survive even under weak quadratic coupling. This study paves the way for probing a previously unexplored domain of quantum optomechanics.

Results
In this section, we discuss the model Hamiltonian for the quadratic interaction in quantum optomechanics. As it will be shown, it is composed of two contributing terms. The first term  1 is the well-known standard quadratic term, resulting from the product of photon number =ˆ † n a a and squared displacement  = +ˆ † b b ( ) 2 2 . The second term  2 , which is also quadratic in order, represents the non-standard term and describes the momentum-field interaction among the squared momentum of the mirror  = −ˆ † b b ( ) 2 2 and squared second quadrature of the electromagnetic field , thus admitting the form   2 2 . It should be mentioned that the phase of â can be arbitrarily shifted by π/2, allowing one to rewrite the latter interaction as  In that sense, the referral to the field quadratures as either  or  is quite arbitrary.
Hamiltonian. For the purpose of this article, we consider that the standard optomechanical interaction  OM vanishes due to appropriate design with = g 0 0 . This not only simplifies the description of the problem and reduces the dimension of basis significantly, but is also favorable from an experimental point of view, since the lowest order surviving interaction now would be quadratic. As mentioned in the above, this criterion can be easily met in superconducting electromechanics by tuning the electroagnetic frequency to the first odd-profiled mechanical mode, or in optomechanics by using arrangements such as the membrane-in-the-middle setup 8,33 . In general, the condition = g 0 0 might not be exactly achieved and it may cause some difficulty in observation of quadratic effects, since standard optomechanical interactions could still mask out the much weaker quadratic interactions. However, for the purpose of current study, we may neglect this effect in the same way is being done by other researchers [18][19][20][21][22][23][24][25] . The total Hamiltonian is thus given by   Here, ω  and Ω ∼ are respectively the bare unperturbed optical and mechanical frequencies. We notice that the relative frequency notation of optical detuning, which is useful in standard optomechanical and quadratic interactions 35 is not to be used here. In (1), furthermore, we assume the existence of only one drive term with the complex amplitude α and frequency ω d . In general, it is possible to assume the existence of multiple drive terms at different frequencies, but this does not alter the mathematical approach under consideration. We furthermore notice that the non-standard quadratic term can be also written as X H P β = − 2 1 2 2 2 with basically no physically significant difference, as mentioned in the above section, too. Hence, the non-standard term can take on either of the sign conventions -as   2 2 , or +− as   2 2 . We shall proceed with the latter. Also, ε is the strength of the standard quadratic interaction and β is the strength of the non-standard quadratic interaction. These are shown to be related as 30 It is easily seen that in the regime of large optical frequency ω Ω ∼ , we get β ≈ 0 and the non-standard term vanishes. This is what has actually been probed in nearly all experiments on quadratic optomechanical interactions so far [18][19][20][21][22][23][24]36 . The large mechanical frequency regime of standard quadratic interactions has been however recently probed 25 and it has been suggested that the roles of optical and mechanical parts are expected to interchange without consideration of the non-standard quadratic effect. Hence, the condition defining the critical mechanical frequency as marks a critical value for the transition border, across which the regimes of large and small mechanical frequency with respect to the given electromagnetic frequency are distinguished. This happens to be remarkably close to the identical frequencies as ω Ω = ∼ , too. It is not difficult to see that the above Hamiltonian with the sign convention taken as +− can be rewritten as are defined and discussed extensively in the preceding articles 30,34,35 , and time-independent non-interacting terms are dropped which are irrelevant to the behavior of system dynamics. Furthermore, the altered effective optical and mechanical frequencies due to the quadratic optomechanical interaction assume different forms, and now read The importance of the non-standard quadratic term  2 30 is that it describes a non-vanishing correction to the field-mirror interaction, beyond simple nonlinear quantum back-action of mirror on the reflected light. In standard picture of quantum optomechanics, the light gets reflected off a displaced mirror which already has shifted the resonance frequency of cavity. Higher-order corrections to combination of these effects give rise to quadratic interactions. But in quadratic quantum optomechanics, either momenta of field and mirror do not apparently come into consideration, or their contributions are somehow lost because of the approximations used in the expansions. It is expected that normally such an interaction should take care of momentum exchange. The exchange and conservation of momenta under standard quadratic interaction is uncertain and actually not quite obvious, since momentum operators do not show up in the interaction. 34 . This basis can be easily seen to be the smallest possible set, with closed commutators, capable of describing the system modeled by (1). The commutation properties of this basis 34 is here given for the sake of convenience All commutators among photonic ˆˆ † c c n { , , } and phononic operators ˆˆ † d d m { , , } are clearly zero. This fact together with the set of commutators (5) establishes the closedness property of our basis. The basis A { } under consideration is called to be second-order, since its operators are all products of two single ladder operators. Other bases of the third-and higher-orders are discussed elsewhere, respectively for standard optomechanics 35 and quadratic interactions 34 .

Basis. Analysis of the Hamiltonian (3) here is going to be based on the six-dimensional space spanned by the basis operators
Furthermore, we will need [ , ] 0 along with the pair of relationships to evaluate the effect of drive and input noise terms later. Assumption of a resonant drive term here requires ω ω = d , which highlights a constant shift in cavity optical frequency because of the presence of quadratic optomechanical interactions. This will greatly simplify the mathematics involved in construction of Langevin equations. These latter relations (5,6) show that the eight-dimensional basis ∪ˆ † a a A { , } { }, which is of the mixed first and second-order, also constitutes a closed commutation relationships. However, for the purpose of this article, the original six-dimensional basis A { } is sufficient.

Langevin Equations. The task of construction of Langevin equations proceeds with the original equation 1-3,26-29
[ , ], where ẑ is an arbitrary operator belonging to the basis, x is any system annihilator associated with the decay rate γ, and x in is the corresponding input field including contributions from both of the deterministic drive and stochastic noise terms. Also, we may set either =x a with γ κ = and =x a in in , or =x b with γ = Γ and =x b in in to allow simple and straightforward construction of noise terms. This particular choice avoids appearance of squared noises, which are otherwise required 34  There are a total of six Langevin equations for the system (1), which can be constructed one by one, corresponding to the six members of A { }. These equations can be obtained using the commutators (5,6) after some straightforward algebra as i t i t where = | | n a 2 is the average cavity photon number and m is the average cavity phonon number. From now on, n and m will represent only deviations from the average steady-state populations. In a similar manner, we can employ the decompositions and replacements It has to be mentioned that the replacements (9,10) are not linearization, but rather an algebraic convention which further allows distinction of steady-state cavity photon and phonon populations under resonant drive. The first separated scalar terms represent the major non-oscillating parts n and m, and oscillating parts for the rest.
One may now proceed to construct the equations for steady-state populations n and m. This requires employing the replacements (9,10) in Langevin equations (8), taking the derivatives on the left, multiplying both sides of the first and fourth equations of (8) by respectively ω e i t 2 and Ω e i t 2 , and discarding all remaining time-dependent terms. This is equivalent to a Rotating-Wave Approximation (RWA) analysis, but carried out at double frequencies. Obviously, all subsequent derivations and calculations in the steady-state will retain their validity within the constraint of RWA condition.
Under resonance conditions where optical and mechanical frequencies are the same ω = ±Ω, there remain a few extra terms. After discarding stochastic noise input, application of these steps to the first, third, fourth, and sixth equations of (8) and simplifying, gives the following four nonlinear algebraic equations for steady state terms Here, δ θ ϕ , represents the Kronecker's delta, and equals 1 if θ ϕ = and 0 otherwise. Also, α is the complex drive amplitude as already defined in the above under (1). These are four equations in terms of n, m, d , and α ∠ , and we notice that  α η ω | | = P / op is the photon flux incident on the cavity due to the optical power P op , where η is the input coupling efficiency. When there is no resonance between optics and mechanics with δ = ω ± ±Ω 0 , , the system (11) reduces to the fairly simple pair Integrable System & Stability. The Langevin equations (8) are still nonlinear and thus non-integrable, but they are instead expressed in terms of second-order operators. Hence, even after taking out the constant oscillating parts using the replacements (9,10), and ignoring the remaining fourth-and higher-order nonlinear terms, the basic nonlinear quadratic interaction still survives. This advantage in using higher-order operators has been noticed by the author 34 as well as others 36 .
Following the one-dimensional analysis provided in the preceding study 30 , the strength of subsequent nonlinear nth-order terms decrease typically as x l ( / ) n zp where x zp is the mechanical zero-point fluctuations of vacuum and l is the typical length scale of the cavity. Having that said, even and odd powers grow almost proportionally. So, once the membrane-in-the-middle or any equivalent setup which allows vanishing g 0 is employed, there would be no non-zero odd-order optomechanical interaction at all. Similarly, the lowest non-vanishing nonlinear term after the quadratic term would be the sextic (6th) order interaction which is expected to be weaker by a factor of x l ( / ) n zp at least. This ratio shall be normally too small to be of any physical significance. Doing this leaves the linearized integrable form  Here, we have furthermore employed the white noise approximation to the Weiner processes â in and b in . This Markovian approximation makes these noise processes insensitive to any frequency shift, or multiplication by purely oscillating terms such as ω ±i t exp( ) and ± Ω i t exp( ). This fact facilitates the construction of noise processes, avoiding the burden of higher-order noise terms. Otherwise, terms such as â in 2 and b in 2 enter the Langevin equations 34 which need a careful and very special treatment to evaluate their corresponding spectral densities. Now, the system of Langevin equations (13) is fully linearized and can be integrated to obtain the spectral densities of each variable. In order to do this, we first define the input vector The Langevin equations now read in where the coefficients matrix M [ ] is given as Here, ζ ε β = + and χ ε β = − . This system is fully integrable and linearized, describing a nonlinear Hamiltonian, and can be Fourier transformed to find the spectral densities. Dynamical stability can be easily determined by inspection of the loci of eigenvalues of (17). The optomechanical system is unconditionally stable if all six eigenvalues have negative real values.
It is remarkable that any linearization on the operators â and b done on the Hamiltonian (1) shall put the interaction in the well recognized forms of either  or , which reverts back to the basic conventional physics of  There are quite a few assumptions needed to obtain (19), including validity of Gaussian white noise processes for both photons and phonons, complete independence of stochastic processes for phonons and photons, and ignoring higher-order noise terms arising from multiplicative noise terms in (8) leading to the mean-field approximation for nonlinear multiplicative noise.
All remains now is to recover the spectral density S w ( ) AA of the first-order operator â from the calculated spectral density S w ( ) CC of the second-order operator. Assuming that S w ( ) CC is composed of well isolated Gaussian peaks as CC j j j 2 2 with s j and ω ∆ j respectively being the amplitude and spread of the j-th Lorentzian, S w ( ) AA can be related 34 720 kHz π × 2 720 kHz π × . 2 5 5kHz π × . 2 2 4Hz 40mK π × 2 5Hz Table 1. Numerical values of resonant studied superconducting electromechanical system parameters with quadratic interaction, in weak coupling regime. Numbers are taken where available from an experiment 33 with membrane-in-the-middle design. In non-resonant cases, ω is varied while Ω is kept constant. For strong and ultrastrong couplings, ε is increased respectively by a factor of 10 2 and 10 3 .    Therefore, the spectral density of the requested system operator can be approximately recovered from the calculated spectral density of the higher-order system operator. A factor 4 has been already been absorbed in (21) because of the definition =ĉ a

Discussion
In this section, we discuss numerical results of solving the system (16) using the scattering matrix formalism (18) and present the computed spectral densities. We base the simulation on the parameters mostly taken from a recent study on superconducting electromechanics 33 where radio-and mechanical frequencies are accurately tuned and set to equal values. However, we study resonant ω = Ω, near-resonant ω ≈ Ω, and off-resonant cases ω > Ω and ω < Ω, in weak, strong, and ultrastrong coupling regimes. The cases of strong and ultrastrong coupling are here somehow idealized, and there might be difficulty in attaining those conditions under practical experimental situations. Nevertheless, these can be obtained using appropriately designed circuit cavity-electrodynamics quantum simulation of quadratic optomechanics 37 . Since no quadratic interaction is investigated in the article under consideration, we proceed to assume that the interaction with the second mechanical mode is being studied so that g 0 may be effectively set to zero. Furthermore, ε is taken to be one to two orders of magnitude below the typical measured g 0 in the system, where g 0 is the single-photon optomechanical interaction rate. Simulation parameters are given in Table 1. The numbers in Table 1 are already chosen as similar as possible to the experiment 33 where available. This is the only known experimental configuration to the author, where mechanical frequency Ω and electromagnetic circuit frequency ω are designed to be within the same order of magnitude, here identical at ω π = Ω = × 2 720kHz, using a carefully designed membrane-in-the-middle setup. Having that said, not only it is actually possible to probe the regime of large mechanical frequency and even ω = Ω, as already has been shown, but also a membrane-in-the-middle configuration could provide direct access to quadratic interactions.
The steady-state populations of photon n and phonons m are the first step to solve for. This can be done by numerical solution of (12) for the near-resonant and non-resonant cases ω ≠ Ω, while for the resonant case with ω = Ω, the set of equations (11) has to be solved for positive real roots. The general behavior is that both n and m increase with input photon flux α as illustrated in Fig. 1 as logarithmic contour plots for various frequency ratios ω Ω/ and input fluxes α. However, making ordinary logarithmic plots at near-resonance ω ≈ Ω in Fig. 2, and off-resonant cases with ω Ω = .
0 56 respectively in Figs 3 and 4 reveals that both populations n and m start to increase beyond a certain threshold. In all cases, under weak coupling with ε π = × z 2 5H , this threshold is around α = 20 Hz sat . While n effectively saturates at high illumination intensities for all non-resonant cases ω ≠ Ω, m tends to increase indefinitely. This behavior due to quadratic interaction is remarkably different from what one would expect in standard optomechanics, where ∝ m n 2 35 . The case of resonant excitation with ω = Ω exhibits a very different behavior, as illustrated in Fig. 5. One could notice that this time, only the photon population n increases almost linearly with the input power with no with (solid black) and without (dashed blue) momentum interaction  2 , for the non-resonant case with ω = Ω 2 and strong coupling ε π = × z 2 500H . Figure 10. Calculated spectral density S w ( /2) CC with (solid black) and without (dashed blue) momentum interaction  2 , for the non-resonant case with ω = Ω 2 and strong coupling ε π = × z 2 500H .
theoretically ideal limit, while phonon population m gets saturated around ≈ m 1 above roughly ten times the same threshold of 10α sat . Quite obviously, too much photon population will cause heating and shift in the working temperature, which respectively will cause a temperature-induced drift in all working parameters. This effect is partially studied in a previous article 35 for zero-point induced spring effect, which is beyond the purpose of current study. However, dynamical stability maps, as shown in Figs 2, 3, 4 and 5 reveal no serious concern in the regimes of interest and all systems are unconditionally stable under the studied range of parameters.
The next sets of plots in Figs 6, 7 and 8 show the calculated spectral density S w ( /2) CC , respectively for resonant and non-resonant cases ω = Ω, ω = Ω 2 , and ω = Ω 2 . The frequency axis is already compressed by a factor of 2 because of (21), which estimates such a halving. In order to study the effect of momentum interaction term  2 , which its existence is predicted in an earlier study 30 , we estimate the spectral densities under two conditions. First, we exclude the momentum interaction and only keep the standard quadratic term  1 . Curves plotted in dahsed blue demonstrated this regime of having only standard quadratic interaction. Curves plotted in solid black correspond to simulation results where the non-standard quadratic interaction with momentum conservation  2 has been included. We observe here the markedly prominent effect of  2 under resonant case with high drive α = 10 Hz 3 shown in Fig. 6. At lower intensites, spectral responses differ in magnitude and  2 has caused an intensifying in quadratic interactions. This significant change in optical spectral response across the mechanical frequency at ω = Ω can be taken as a primary signature of existence of momentum interaction term  2 , predicted earlier 30 .
The behavior at optical frequencies ω > Ω is quite expected. Since from (2) we have β ε ω ∝ Ω / ( / ) 2 , one would estimate that at increased optical frequencies, the effect of non-standard interaction  2 would start to fade out. This is quite correct, indeed, as at high optical frequencies ω Ω  , we notice the gradually disappearing effect of non-standard quadratic term  2 , as illustrated in Fig. 7 for ω = Ω 2 . When larger values for the ratio ω Ω / is selected, the effect of  2 is completely wiped out. This fact also confirms why the non-standard momentum interaction  2 has not been observed in experiments so far, as all of the quadratically studied optomechanical interactions fall in this regime of ω Ω  . Also, we may identify the approximate limit of ω Ω ≤ ( / ) 10 2 beyond which there would be no observable non-standard quadratic interaction.
However, when the mechanical frequency gets large with respect to the optical frequency ω Ω > , the effect of non-standard interaction  2 becomes quite evident. This has been shown in Fig. 8. First of all, momentum interactions tend to reduce the overall quadratic interaction under studied conditions. Secondly, at higher illumination and photon flux rates, we notice a clear displacement of the peak toward higher frequencies because of the with (solid black) and without (dashed blue) momentum interaction  2 , for the resonant case with ω = Ω and strong coupling ε π = × z 2 500H . The non-standard quadratic interaction completely wipes out the sharp resonance and peak in the spectrum. non-standard quadratic interaction  2 . This fact can serve as a secondary and unmistakable signature of existence of momentum interaction terms in an experiment. Next, we turn to the case of strongly coupled quadratic optomechanics with ε π = × z 2 500H . This regime has been studied in the next sets of plots included in Figs 9, 10 and 11. This time, the effect of strong coupling on the quadratic response is much pronounced, as opposed to the weak coupling. We notice that the behavior of steady state photon n and phonon m populations is essentially similar to those of weakly coupled regime, with the difference that the flux threshold for saturation or increase is shifted to the higher value of roughly α = 200Hz sat by an order of magnitude. The corresponding graphs are not shown for the sake of brevity. Again, one could observe that at higher optical frequencies ω > Ω, the effect of momentum exchange term  2 is much smaller than the standard quadratic interaction  1 . This fact is shown in plots of Fig. 9. Meanwhile, at small optical frequencies where mechanical frequency prevails ω Ω > , non-standard interactions  2 can have a dominant effect on the response. This is more or less similar to the weak coupling regime show in Fig. 7. However, the overall spectral response is actually less intense, since the associated peaks are smaller. Hence, increasing the strength of quadratic interaction does not necessarily increase the intensity of exhibited response. Furthermore, one could notice a marked shift in spectrum toward higher frequencies at sufficiently high photon input flux.
The case of large mechanical frequency ω Ω > under the strong coupling regime is also studied in Fig. 10, as opposed to the weak coupling regime in Fig. 8. The general trend of behavior from weak to strong coupling regimes at α = 10 Hz 3 has not essentially changed, and there is an expected shift in the peak in spectrum toward higher frequencies because of non-standard interactions. However, continuing to increase the input optical intensity causes complete and distinct separation among the corresponding peaks, while both move toward lower frequencies. This reveals the fact that a tertiary signature of momentum exchange term  2 can be taken as the central location of peak in response, which is dependent on the possible existence of  2 .
Finally, we studied the resonant case ω = Ω under strong coupling, shown in Fig. 11 as opposed to the results of weak coupling shown in Fig. 6. We notice a very distinct difference between these two cases. It can be observed that under strong coupling, the standard quadratic response is expected to exhibit a very sharp resonance exactly at ω = w , which is actually placed a bit higher than the central peak frequency of the output spectrum. The existence of these resonances in clearly shown by zoomed out plots of Fig. 12. One could tell, that existence of a tiny squeezing even could have been expected under standard quadratic interaction  1 , which is again entirely wiped out by the overriding presence of non-standard momentum interaction term  2 . This can be taken as the existence of quaternary signature of momentum-field interactions, which eliminates such squeezing and sharp resonances from the spectral response. The behavior under ultrastrong coupling has been calculated studied, but since it is essentially similar to the strongly coupled case, the results are not shown and dropped from the article.
The results of this research could be relevant to sensing applications where nonlinear optomechanics could in principle produce pronounced sensitivities, such as those reported recently for quantum gravimetry 17,32 . Also possible generation of continuous squeezed electromagnetic radiation seems to be feasible by carefully optimized design of a membrane-in-the-middle superconducting electromechanics setup. Further applications could be still plausible but will need more in depth study. Experiments can unambiguously determine the existence of non-standard quadratic optomechanics.

Conclusions
In summary, we presented a detailed theoretical and numerical analysis of the quadratic optomechanical interactions, using the method of higher-order operators. We have studied both types of standard quadratic and the predicted non-standard quadratic interactions, and through extensive simulations have established clear signatures for existence of non-standard interactions, while stability of the system under study has also been established. Such types of quadratic interactions can be probed in a carefully designed experiment where the optical and mechanical frequencies fall within the same order of magnitude.