Higher-Order Interactions in Quantum Optomechanics: Analysis of Quadratic Terms
Sina Khorasani
quant-ph





 
Vienna Center for Quantum Science and Technology, University of Vienna, 1090 Vienna, Austria; sina.khorasani@ieee.org




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. 




Higher-Order Interactions in Quantum Optomechanics: Analysis of Quadratic Terms
    Sina Khorasani
    December 30, 2023
===============================================================================

empty



§ INTRODUCTION


The field of quantum optomechanics <cit.> is flourishing as one of the modern applications of quantum physics, where interactions of optical field and mechanical motion inside a confined cavity is being studied. Nonlinear interactions in optomechanics play a critical role in a growing number of studied physical phenomena, which include emergence of second-order mechanical side-bands <cit.>, nonlinear optomechanical induced transparency <cit.>, phonon-laser <cit.>, nonlinear reciprocity <cit.>, and parameteric phonon-phonon coupling for cascaded optical transparency <cit.>. Recently, existence of optomechanical chaos has been confirmed <cit.>, which is contingent on and can be explained only using the nonlinearity of interactions. Meterology with optomechanical symmetry breaking <cit.> as well as Kerr-type nonlinear interactions <cit.> are other examples of optomechanical phenomena which call for nonlinear analysis. Interestingly, it has been also shown using an extensive analysis, that the gravitational constant can be measured by increased precision using nonlinear optomechanics <cit.>, reaching an ideal fundamental sensitivity of 10^-15ms^-2 for state-of-the-art parameters.

Many of the studies in this field utilize linearization of photonic â and phononic  b̂ ladder operators around their mean values as â→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 <cit.> and higher-order interactions where the resulting Langevin equations <cit.> are expected to be strongly nonlinear. Full linearization essentially transforms back every such nonlinear interaction picture into the simple linearized form of (â+â^†)(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 <cit.> as well as others for nonlinear <cit.> and quadratic <cit.> 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 <cit.>. 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 <cit.>, 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 <cit.> 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 <cit.>. 

Growing out of the context of quantum optomechanics, the method of higher-order operators developed by the author <cit.> 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 <cit.>, 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 <cit.> 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̂=â^†â and squared displacement 𝕏^2=(b̂+b̂^†)^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 ℙ^2=(b̂-b̂^†)^2 and squared second quadrature of the electromagnetic field 𝒫^2=(â-â^†)^2 <cit.>, 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 𝒳^2ℙ^2, where 𝒳^2=(â+â^†)^2. 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 <cit.>. 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 <cit.>. The total Hamiltonian is thus given by ℍ=ℍ_0+ℍ_1+ℍ_2+ℍ_d where

    ℍ_0   =   ħω̃â^†â+ħΩ̃b̂^†b̂=ħω̃n̂+ħΩ̃m̂,
    ℍ_1   =   ħ12εâ^†â(b̂+b̂^†)^2=ħ12εn̂𝕏^2, 
    ℍ_2   =   -ħ12β (â-â^†)^2(b̂-b̂^†)^2=-ħ12β𝒫^2ℙ^2, 
    ℍ_d   =   ħγ(α e^iω_d tâ+α^∗ e^-iω_d tâ^†).

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 <cit.> is not to be used here. In (<ref>), 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 ℍ_2=-ħ12β𝒳^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 <cit.>

    β=1/4(π^2/3+1/4)(Ω̃/ω̃)^2ε≡12ρε.

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 <cit.>. The large mechanical frequency regime of standard quadratic interactions has been however recently probed <cit.> 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 Ω̃=1/2√(1/3π^2+1/4)ω̃≈ 0.941ω̃ 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

    ℍ=ħωn̂+ħΩm̂+ħϵ[n̂(m̂+d̂+d̂^†)-ρ(n̂+ĉ+ĉ^†)(d̂+d̂^†-m̂)]+ℍ_d,

in which ĉ=1/2â^2 and d̂=1/2b̂^2 are defined and discussed extensively in the preceding articles <cit.>, 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

    ω=ω̃+12ε+β,   Ω=Ω̃+β.


The importance of the non-standard quadratic term ℍ_2 <cit.> 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. 



 §.§ Basis


Analysis of the Hamiltonian (<ref>) here is going to be based on the six-dimensional space spanned by the basis operators {A}^T={ĉ,ĉ^†,n̂,d̂,d̂^†,m̂} <cit.>. This basis can be easily seen to be the smallest possible set, with closed commutators, capable of describing the system modeled by (<ref>). The commutation properties of this basis <cit.> is here given for the sake of convenience

    [ĉ,ĉ^†]=n̂+12,     
    [ĉ,n̂]=2ĉ,      
    [n̂,ĉ^†]=2ĉ^†, 
    
    [d̂,d̂^†]=m̂+12,     
    [d̂,m̂]=2d̂,      
    [m̂,d̂^†]=2d̂^†.

All commutators among photonic {ĉ,ĉ^†,n̂} and phononic operators {d̂,d̂^†,m̂} are clearly zero. This fact together with the set of commutators (<ref>) 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 <cit.> and quadratic interactions <cit.>.

Furthermore, we will need [ĉ,â]=[ĉ^†,â^†]=0 along with the pair of relationships

    [ĉ,â^†]=[â,n̂]=â,   
    [â,ĉ^†]=[n,â^†]=â^†,

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 (<ref>,<ref>) show that the eight-dimensional basis {â,â^†}⋃{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 <cit.>

    ẑ̇=-i/ħ[ẑ,ℍ]-[ẑ,x̂^†](12γx̂+√(γ)x̂_in)+(12γx̂^†+√(γ)x̂_in^†)[ẑ,x̂],

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̂=â with γ=κ and x̂_in=â_in, or x̂=b̂ with γ=Γ and x̂_in=b̂_in to allow simple and straightforward construction of noise terms. This particular choice avoids appearance of squared noises, which are otherwise required <cit.>.

There are a total of six Langevin equations for the system (<ref>), which can be constructed one by one, corresponding to the six members of {A}. These equations can be obtained using the commutators (<ref>,<ref>) after some straightforward algebra as

    ĉ̇   =   -i2ωĉ-κĉ-iε[2(m̂+d̂+d̂^†)ĉ-ρ(d̂+d̂^†-m̂)(2ĉ+n̂+12)]+iα^∗ e^-iω tâ-√(κ)ââ_in ,
    ĉ̇^†   =   i2ωĉ^†-κĉ^†+iε[2(m̂+d̂+d̂^†)ĉ^†-ρ(d̂+d̂^†-m̂)(2ĉ^†+n̂+12)]-iα e^iω tâ^†-√(κ)â_in^†â^† ,
    ṅ̂̇   =   -i2β(ĉ-ĉ^†)(d̂+d̂^†-m̂)-κn̂+i(α e^iω tâ-α^∗ e^-iω tâ^†)-√(κ)(â^†â_in+â_in^†â), 
    ḋ̂̇   =   -i2Ωd̂-Γd̂-iε[n̂(2d̂+m̂+12)-ρ(ĉ+ĉ^†+n̂)(m̂+12-2d̂)]-√(Γ)b̂b̂_in ,
    ḋ̂̇^†   =   i2Ωd̂^†-Γd̂^†+iε[n̂(2d̂^†+m̂+12)-ρ(ĉ+ĉ^†+n̂)(m̂+12-2d̂^†)]-√(Γ)b̂_in^†b̂^† ,
    ṁ̂̇   =   -i2ε(d̂-d̂^†)[(ρ-1)n̂+ρ(ĉ+ĉ^†)]-√(Γ)(b̂^†b̂_in+b̂_in^†b̂).




 §.§ Steady-State Equilibrium

The special algebraic form of non-standard interaction together with the presence of drive term makes the analysis requiring a bit of care. Ambiguities could be avoided by using explicit decompositions and replacements

    â(t)→a̅e^-iω t+â(t),         â^†(t)→a̅^∗ e^iω t+â^†(t), 
    n̂(t)→n̅+n̂(t),         m̂(t)→m̅+m̂(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

    ĉ(t)→c̅e^-2iω t+ĉ(t),         ĉ^†(t)→c̅^∗ e^2iω t+ĉ^†(t), 
    d̂(t)→d̅e^-2iΩ t+d̂(t),         d̂^†(t)→d̅^∗ e^2iΩ t+d̂^†(t).

It has to be mentioned that the replacements (<ref>,<ref>) 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 (<ref>,<ref>) in Langevin equations (<ref>), taking the derivatives on the left, multiplying both sides of the first and fourth equations of (<ref>) by respectively e^i2ω t and e^i2Ω t, 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 (<ref>) and simplifying, gives the following four nonlinear algebraic equations for steady state terms

    iα^∗√(n̅)   =   12κn̅+iε[12(1+ρ)m̅n̅-ρ(n̅+12)(δ_ω,Ωd̅+δ_ω,-Ωd̅^∗)],
    m̅^2-m̅   =   4|d̅|^2,
    
    -Γd̅   =   iε[2(1+ρ)n̅d̅-ρ(m̅+12)(δ_ω,Ω+δ_ω,-Ω)n̅], 
    κn̅   =   i(α-α^∗)√(n̅)-i2β(d̅^∗-d̅)n̅(δ_ω,Ω-δ_ω,-Ω).

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 (<ref>). 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 (<ref>) reduces to the fairly simple pair 

    |α|^2=ε^2(1+ρ)^2m̅^2n̅,   n̅|α|^2=ε^2(m̅^2-m̅).

These two nonlinear algebraic equations can be numerically solved for non-negative real roots of n̅ and m̅.



 §.§ Integrable System & Stability


The Langevin equations (<ref>) 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 (<ref>,<ref>), 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 <cit.> as well as others <cit.>. 

Following the one-dimensional analysis provided in the preceding study <cit.>, the strength of subsequent nonlinear nth-order terms decrease typically as (x_zp/l)^n 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_zp/l)^n at least. This ratio shall be normally too small to be of any physical significance. Doing this leaves the linearized integrable form

    ĉ̇   =   [i2ω-κ-i2(ε+β)m̅]ĉ+iβ(n̅+12)(d̂+d̂^†-m̂)-iβm̅n̂-√(12n̅κ)â_in, 
    ĉ̇^†   =   [-i2ω-κ+i2(ε+β)m̅]ĉ^†-iβ(n̅+12)(d̂+d̂^†-m̂)+iβm̅n̂-√(12n̅κ)â_in^†, 
    ṅ̂̇   =   -κn̂+i2βm̅(ĉ-ĉ^†)-√(κn̅)(â_in+â_in^†),
    ḋ̂̇   =   [-i2Ω-Γ-i2(ε+β)n̅]d̂+iβ(m̅+12)(ĉ+ĉ^†)-i(ε-β)[n̅m̂+(m̅+12)n̂]-√(2Γ|d̅|)b̂_in,
    ḋ̂̇^†   =   [i2Ω-Γ+i2(ε+β)n̅]d̂^†-iβ(m̅+12)(ĉ+ĉ^†)+i(ε-β)[n̅m̂+(m̅+12)n̂]-√(2Γ|d̅|)b̂_in^†,
    ṁ̂̇   =   -Γm̂+i2(ε-β)n̅(d̂-d̂^†)-√(Γ|d̅|)(b̂_in+b̂_in^†).
 
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 exp(± iω t) and exp(± iΩ t). 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 <cit.> which need a careful and very special treatment to evaluate their corresponding spectral densities. Now, the system of Langevin equations (<ref>) 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

    {A_in}^T={â_in,
    â_in^†, 
    12(â_in+â_in^†), 
    b̂_in,
    b̂_in^†,
    12(b̂_in+b̂_in^†) 
    },

and the diagonal matrix

    [γ]=Diag[n̅κ,n̅κ,4n̅κ,2|d̅|Γ,2|d̅|Γ,4|d̅|Γ].

The Langevin equations now read

    d/dt{A}=[M]{A}-√([γ]){A_in},

where the coefficients matrix [M] is given as

	
    [M]= 
    	[
    	[  i2(ω-ζm̅)-κ            0        -iβm̅   iβ(n̅+1/2)   iβ(n̅+1/2)  -iβ(n̅+1/2);            0 -i2(ω-ζm̅)-κ         iβm̅  -iβ(n̅+1/2)  -iβ(n̅+1/2)   iβ(n̅+1/2);        i2βm̅       -i2βm̅           -κ            0            0            0;   iβ(m̅+1/2)   iβ(m̅+1/2)  -iχ(m̅+1/2) -i2(Ω+ζn̅)-Γ            0        -iχn̅;  -iβ(m̅+1/2)  -iβ(m̅+1/2)   iχ(m̅+1/2)            0  i2(Ω+ζn̅)-Γ         iχn̅;            0            0            0        i2χn̅       -i2χn̅           -Γ ]],


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 (<ref>). 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 (<ref>) shall put the interaction in the well recognized forms of either 𝒳𝕏 or 𝒫ℙ, which reverts back to the basic conventional physics of optomechanics <cit.>. Hence, any linearization process on the first-order operators will eventually discard the physics of higher-order interactions. The method of higher-order operators discussed in this article and previous researches <cit.> is targeted to alleviate this problem.



 §.§ Spectral Density


The input-output relations <cit.> connects the incident and reflection waves of a weak probe beam at frequency w (not to be confused with optical frequency ω). This relationship in Fourier domain together with (<ref>) yield

    {A_out(w)}   =   {A_in(w)}+√([γ]){A(w)},
    {A_out(w)}   =   [S(w)]{A_in(w)},
    
    [S(w)]   =   [I]-√([γ])([M]-iw [I])√([γ]),

where [I] is the identity matrix. This can be now summarized to yield the spectral density of ĉ for positive frequencies w>0 as

    S_CC(w)=|S_11(w)|^2+|S_13(w)|^212+|S_14(w)|^2(m̅+1)+|S_15(w)|^2m̅+|S_14(w)|^2(m̅+12),

where S_ij(w) are frequency dependent elements of the 6× 6 scattering matrix [S(w)] defined in (<ref>). In the above, we have made use of the fact that spectral density of the input vector {Â_in(w)} (<ref>) is 

    { S_in(w>0)}^ T   =   {1,0,12,m̅+1,m̅,m̅+12}, 
    { S_in(w<0)}^ T   =   {0,1,12,m̅,m̅+1,m̅+12}.

There are quite a few assumptions needed to obtain (<ref>), 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 (<ref>) leading to the mean-field approximation for nonlinear multiplicative noise. 

All remains now is to recover the spectral density S_AA(w) of the first-order operator â from the calculated spectral density S_CC(w) of the second-order operator. Assuming that S_CC(w) is composed of well isolated Gaussian peaks as S_CC(w)=∑ s_j exp[-(w-ω_j)^2/Δω_j^2] with s_j and Δω_j respectively being the amplitude and spread of the j-th Lorentzian, S_AA(w) can be related <cit.> to S_CC(w) roughly as

    S_AA(w)≈ 2π√(π)∑_j s_jω_j^2/Δω_jexp[-(w-1/2ω_j)^2/2Δω_j^2].

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 (<ref>) because of the definition ĉ=1/2â^2.



§ DISCUSSION


In this section, we discuss numerical results of solving the system (<ref>) using the scattering matrix formalism (<ref>) and present the computed spectral densities. We base the simulation on the parameters mostly taken from a recent study on superconducting electromechanics <cit.> 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 <cit.>.

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 <ref>. The numbers in Table <ref> are already chosen as similar as possible to the experiment <cit.> 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 (<ref>) for the near-resonant and non-resonant cases ω≠Ω, while for the resonant case with ω=Ω, the set of equations (<ref>) 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. <ref> as logarithmic contour plots for various frequency ratios Ω/ω and input fluxes α. However, making ordinary logarithmic plots at near-resonance ω≈Ω in Fig. <ref>, and off-resonant cases with Ω=1.78ω and Ω=0.56ω respectively in Figs. <ref> and <ref> reveals that both populations n̅ and m̅  start to increase beyond a certain threshold. In all cases, under weak coupling with ε=2π× 5 Hz, this threshold is around α_sat=20Hz. 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 <cit.>. 

The case of resonant excitation with ω=Ω exhibits a very different behavior, as illustrated in Fig. <ref>. One could notice that this time, only the photon population n̅ increases almost linearly with the input power with no 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 <cit.> for zero-point induced spring effect, which is beyond the purpose of current study. However, dynamical stability maps, as shown in Figs. <ref>-<ref> 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. <ref>-<ref> show the calculated spectral density S_CC(w/2), respectively for resonant and non-resonant cases ω=Ω, ω=2Ω, and 2ω=Ω. The frequency axis is already compressed by a factor of 2 because of (<ref>), 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 <cit.>, 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^3Hz shown in Fig. <ref>. 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 <cit.>. 

The behavior at optical frequencies ω>Ω is quite expected. Since from (<ref>) 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. <ref> 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 (ω/Ω)^2≤ 10 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. <ref>. 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 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 ε=2π× 500 Hz. This regime has been studied in the next sets of plots included in Figs. <ref>-<ref>. 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 α_sat=200Hz 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. <ref>. 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. <ref>. 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. <ref>, as opposed to the weak coupling regime in Fig. <ref>. The general trend of behavior from weak to strong coupling regimes at α=10^3Hz 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. <ref> as opposed to the results of weak coupling shown in Fig. <ref>. 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. <ref>. 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 <cit.>.  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. 



§ ACKNOWLEDGMENTS

	
Stimulating discussions of this work with Dr. David Edward Bruschi at the University of Vienna, Prof. Radim Filip at Palacký University, members of the Laboratory of Quantum Foundations and Quantum Information at Nano- and Micro-scale in Vienna Center for Quantum Science and Technology (VCQ), as well as members of the Laboratory of Photonic and Quantum Measurements at École Polytechnique Fédérale de Lausanne (EPFL) is much appreciated. This work is dedicated to the celebrated artist, Anastasia Huppmann.
	


§ ADDITIONAL INFORMATION

The author declares no competing interests.

99
	Kip
	Kippenberg, T. J. & Vahala, K. J. Cavity optomechanics. Science 321, 1172 (2008).
	Aspel
	Aspelmeyer, M., Kippenberg, T. J. & Marquardt, F. Cavity optomechanics. Rev. Mod. Phys. 86, 1391 (2014).
	Bowen
	Bowen, W. P. & Milburn, G. J. Quantum Optomechanics (CRC Press: Boca Raton, 2016).


q1
Xiong, H., Si, L.-G., Zheng, A.-N., Yang, X. & Wu, Y. Higher-order sidebands in optomechanically induced 
transparency. Phys. Rev. A 86, 013815 (2012). 
q12
Jiao, Y. F., Lu, T. X. & Jing, H. Optomechanical second-order sidebands and group delays in a Kerr resonator. 
Phys. Rev. A 97, 013843 (2018). 

q2
Lemonde, M.-A., Didier, N. & Clerk, A. A. Nonlinear interaction effects in a strongly driven optomechanical cavity. Phys. Rev. Lett. 111, 053602 (2013). 

q3
Kronwald, A. & Marquardt, F. Optomechanically induced transparency in the nonlinear quantum regime. Phys. 
Rev. Lett. 111, 133601 (2013). 

q4
Karuza, M., et al. Optomechanically induced transparency in a membrane-in-the-middle setup at room temperature. Phys. Rev. A 88, 013804 (2013). 

q8
Shen, Z., et al. Compensation of the Kerr effect for transient optomechanically induced transparency in a silica microsphere. Opt. Lett. 41, 1249 (2016).

q5
Jing, H., et al. PT-symmetric phonon laser. Phys. Rev. Lett. 113, 053604 (2014).
q13
Zhang, J., et al. L. A phonon laser operating at an exceptional point. Nat. Photon. 12, 479 (2018).

q6
Zhang, J., et al. Giant nonlinearity via breaking parity-time symmetry: A route to low-threshold phonon diodes. Phys. Rev. B 92, 115407 (2015).

q7
Fan, L., Fong, K. Y., Poot, M. & Tang, H. X. Cascaded optical transparency in multimode-cavity optomechanical 
systems. Nat. Commun. 6, 5850 (2015). 


q9
Monifi, F., et al. L. Optomechanically induced stochastic resonance and chaos transfer between optical fields. Nat. Photon. 10, 399 (2016).

q10
Liu, Z.-P., et al. Metrology with PT-symmetric cavities: Enhanced sensitivity near the PT-phase transition. Phys. Rev. Lett. 117, 110802 (2016).

q11
Mikkelsen, M., Fogarty, T., Twamley, J. & Busch, T. Optomechanics with a position-modulated Kerr-type nonlinear coupling. Phys. Rev. A 96, 043832 (2017). 


q14
Qvarfort, S., Serafini, A., Barker, P. F. & Bose, S. Gravimetry through non-linear optomechanics. Nat. Commun. 9, 3690 (2018).


	Quad1
	Zhang, L., Ji, F., Zhang, X. & Zhang, W. Photon-phonon parametric oscillation induced by quadratic coupling in an optomechanical resonator. J. Phys. B 50, 145501 (2017).
	Quad3
	Nunnenkamp, A., Børkje, K., Harris, J. G. E. & Girvin, S. M. Cooling and squeezing via quadratic optomechanical coupling. Phys. Rev. A 82, 021806 (2010).
	Quad4
	Dalafi, A., Naderi, M. H. & Motazedifard, A. Effects of quadratic coupling and squeezed vacuum injection in an optomechanical cavity assisted with a Bose-Einstein condensate.	Phys. Rev. A 97, 043619 (2018).
	Quad5 
	Asjad, M., et al. Robust stationary mechanical squeezing in a kicked quadratic optomechanical system. Phys. Rev. A 89, 023849 (2014).
	Quad6 
	Liao, J. Q. & Nori, F. Photon blockade in quadratically coupled optomechanical systems. Phys. Rev. A 88, 023853 (2013).
	Quad7 
	Zhan, X.-G., Si, L.-G., Zheng, A.-S. & Yang, X. X. Tunable slow light in a quadratically coupled optomechanical system. J. Phys. B 46,  025501 (2013).
	Quad8
	Doolin, C., et al. Nonlinear optomechanics in the stationary regime. Phys. Rev. A 89, 053838 (2014).
	Bruschi
	Bruschi, D. E. & Xuereb, A. Mechano-optics: An optomechanical quantum simulator. New J. Phys. 20, 065004 (2018).
	
	Langevin1 
	Gardiner, C. W. & Zoller, P. Quantum Noise (Springer: Berlin, 2004).
	Langevin2 
	Gardiner, C. W. & Collett, M. J. Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation. Phys. Rev. A 31, 3761 (1985).
	Langevin3
	Gardiner, C. & Zoller, P. The Quantum World of Ultra-Cold Atoms and Light. Book I: Foundations of Quantum Optics (Imperial College Press: London, 2014).
	Langevin4
	Combesa, J., Kerckhoff, J. & Sarovar, M. The SLH framework for modeling quantum input-output networks. Adv. Phys.: X 2, 784 (2017).
	
	Paper1
	Khorasani, S. Higher-order interactions in quantum optomechanics: Revisiting theoretical foundations. Appl. Sci. 7, 656 (2017).
	
	Sala
	Sala, K. & Tufarelli, T. Exploring corrections to the Optomechanical Hamiltonian. Sci. Rep. 8,  9157 (2018).
	
	Thesis
	Pang, B. H. Theoretical Foundations for Quantum Measurement in a General Relativistic Framework. Ph.D. Dissertation, California Institute of Technology; doi:10.7907/dfyy-y188 (2018). 
	
	Quad2
	Bagci, T., et al. Optical detection of radio waves through a nanomechanical transducer. Nature 507, 81 (2014).

	Paper2
	Khorasani, S. Higher-order interactions in quantum optomechanics: Analytical solution of nonlinearity. Photonics 4, 48 (2017).
	Paper3
	Khorasani, S. Method of higher-order operators for quantum optomechanics. Sci. Rep. 8, 11566 (2018).
	Quad9
	Liu, S., Yang, W.-X., Shui, T., Zhu, Z. & Chen, A.-X. Tunable two-phonon higher-order sideband amplification in a quadratically coupled optomechanical system. Sci. Rep. 7, 17637 (2017).
	
	QuadCir
	Kim, E.J., Johansson, J.R. & Nori, F. Circuit analog of quadratic optomechanics. Phys. Rev. A 91, 033835 (2015).
	


