Method of Higher-Order Operators for Quantum Optomechanics
Sina Khorasani
quant-ph

	
	Variational Bayes Estimation of Discrete-Margined Copula Models with Application to  Time Series
    Rubén Loaiza-Maya and Michael Stanley Smith
    First Version December 2017
This Version December 30, 2023
================================================================================================

empty
	
	

§ INTRODUCTION

	
	Nonlinear quantum interactions with stochastic noise input stand among the most difficult analytical challenges to solve in the context of stochastic differential equations. While linearized interactions remain accurate for description of many experiments, a certain class of quadratic and higher-order physical phenomena cannot be normally understood under linearized approximations. While in classical problems the resulting Langevin equations are scalar functions, in quantum problems one has to deal with nonlinear operator differential equations. If expanded unto base kets, bosonic operators can assume infinite-dimensional matrix forms, rendering the solution entirely intractable.
	
	Such classes of nonlinear operator problems can be addressed by construction of Fokker-Planck or nonlinear Schrödinger equations, among which there exists a one-to-one correspondence. The Fokker-Planck equation <cit.> is actually equivalent to the nonlinear Schrödinger equation with bosonic operator algebra, and its moments <cit.> translate into nonlinear Langevin equations. The method of master equations <cit.> also can be used in combination with the quasi-probablity Wigner functions <cit.> to deal with nonlinear quantum interactions. The master equation approach is reasonably accurate as long as Born and Markov approximations are not employed <cit.>. But none of these methods is probably as convenient as the method of Langevin equations <cit.>, which has found popularity in the context of quantum optoemchanics <cit.>. 
	
	Being an inherently nonlinear interaction among photonic and phononic baths <cit.>, the standard quantum optomechanics is normally described by linearized Langevin equations <cit.>. This will suffice to address a majority of complex experimental situations such as optomechanical-induced transparency <cit.> and polaron anti-crossing <cit.>, but effects such as non-classical states of light <cit.>, optomechanical emission of real photons from vacuum <cit.>, photon blockade <cit.>, nonlinear self-oscillations <cit.>, and chaos <cit.> are all among manifestations of nonlinear regimes in standard optomechanics, which need description using nonlinear algebra. Also, biquadratic interactions (mostly referred to as quadratic interactions) among bosonic baths remain a hurdle. In quadratic optomechanics <cit.>, which is a topic of growing interest in the recent year, having an analytical tool capable of addressing such kinds of nonlinearity is advantageous. A perturbation technique based on the expansion of time-evolution operators <cit.> is employed to investigate quadratic interactions and it has been shown that for mechanical frequencies exceeding optical frequencies a new unexplored regime appears in which the roles of optical and mechanical partitions are interchanged.
	
	Recently, the author has reconsidered the theoretical description of optomechanics <cit.> and shown that quadratic interactions are subject to two corrections resulting from momentum conservation and relativistic effects. Such types of quadratic corrections become significant when the mechanical frequency is within the same order of or exceeds electromagnetic frequency. Furthermore, an analytical approach is proposed to tackle nonlinear quantum interactions <cit.> and a method of expansion unto higher-order operators is proposed and investigated in details.
	
	In this article, the higher-order operator approach recently proposed by the author <cit.> is employed to address the standard optomechanics, and it is shown that there exists a minimal choice of higher-order operator basis which leads to exactly linear and fully separable Langevin equations with multiplicative input noise terms <cit.>. We also present a full mathematical treatment of multiplicative noise terms, which turn out to play a crucial rule in higher-order quantum optomechanics. This allows one to provide an exact and explicit solution using an operator-based method to solve the optomechanical interactions in the nonlinear regime. There exists higher-order effects appearing at high optical pump rates, and can be predicted using the method discussed here. These include inequivalent red and blue detunings, higher-order resonance shift and spring effects, and also zero-point-field induced optomechanical shift of mechanical frequency. The inequivalency of red and blue detuned side-bands, which appears as a counter-intuitive difference in their respective frequency shifts, is different from the well-known anomalous Stokes-Anti-Stokes symmetry breaking <cit.> which is connected to different scattering amplitudes. The same method of higher-order operator algebra has been recently used independently as well <cit.>. 
	
	We also show for the first time that the introduced method of higher-order operators can be used to estimate the coherent population of phonons in the optomechanical cavity, here referred to as the coherent phonon number. This quantity can not only be calculated explicitly in terms of optomechanical parameters, but also, can be found by fitting the expressions of corrected spring effect to the experimental observations. Also, dynamic linear and nonlinear stability phases in red and blue-detuned drives can be well computed and estimated using the method of higher-order operators.
	
	

§ RESULTS

	
	The standard optomechanical Hamiltonian reads <cit.>
	
    ℍ_OM=ħΩm̂-ħΔn̂-ħ g_0n̂(b̂+b̂^†),

	where n̂=â^†â and m̂=b̂^†b̂ are photon and phonon number operators with â and b̂ respectively being the photon and phonon annihilators, Ω is the mechanical frequency, Δ is optical detuning from cavity resonance, and g_0 is the single-photon optomechanical interaction rate. The interaction ℍ_OM is not quadratic, but is still cubic nonlinear. It is normally solved by a straightforward linearization <cit.>, but can be also solved at the second-order accuracy using the higher-order operators described in the preceding article <cit.>. 
	
	In order to form a closed basis of operators, we may choose either the higher-order operators
	
    {A}^T={â,âb̂,âb̂^†},

	of the second-degree, which forms a 3× 3 system of Langevin equations, or
	
    {A}^T={â,b̂,âb̂,âb̂^†,n̂,ĉ},

	which forms a 6× 6 system of Langevin equations. Here, we adopt the definition ĉ=1/2â^2 <cit.>. 
	
	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 [â,n̂]=-[âb̂^†,b̂]=â, [âb̂,n̂]=âb̂, [âb̂^†,n̂]=âb̂^†, and [âb̂,âb̂^†]=[ĉ,n̂]=2ĉ, which is obviously a closed basis. Now, one may proceed with composition of the Langevin equations. 
	
	The applicability of the basis (<ref>) becomes readily clear by calculating the braket [â,ℍ_OM] as appears in the corresponding Langevin equation. The terms involving the second-degree operators âb̂ and âb̂^† immediately show up. The key in the method of higher-order operators is to keep these operator pairs, triplets and so on together, as each combination has a clear corresponding physical process. While â and b̂ refer to individual ladder operators, âb̂ and âb̂^† respectively construct the blue and red 1-photon/1-phonon processes. For this reason, it is probably more appropriate to call these higher-degree operator combinations as processes. 
	
	The Langevin equations for the blue âb̂ and red âb̂^† processes do not close on themselves, because of the appearance of third-order blue- and red-like processes âb̂^2 and âb̂^† 2, describing 1-photon/2-phonon processes. Similarly, every j-th order blue- or red-like process such as âb̂^j and âb̂^† j will lead to the j+1-order process. Hence, the infinite-dimensional basis {â}∪{∀âb̂^j,âb̂^† j;j∈𝒩} can provide an exact solution to the optomechanics. Furthermore, the convergence of solutions basis on such expansions would be questionable when g_0<<Ω is violated. In general, the j-th order processes correspond to the 1-photon/j-phonon interactions and contribute to the j+1-order sidebands. In this article, it has been shown that under practical conditions, it is unnecessary to take account of the processes j≥ 2 and the 3× 3 basis (<ref>) is rather sufficient for most practical purposes. Nonetheless, j=2 processes contribute significantly to nonlinear stability and second-order mechanical sidebands. While the use of an infinite-dimensional basis is surprisingly unnecessary in still a higher-order formulation, using the compact minimal basis to be discussed in the following can lead to the mathematically exact solution.
    The choice of basis is not unique, and every non-degenerate linear combination of bases leads to another equivalent form. One may for instance arbitrate the three-dimensional linear basis {A}^ T={â,b̂,b̂^†} or the four-dimensional linear basis {A}^ T={â,b̂,â^†,b̂^†} as is taken in the context of linearized standard optomechanics <cit.>, the five-dimensional all-Hermitian basis {A}^ T={n̂,m̂,n̂^2,n̂(b̂+b̂^†),in̂(b̂-b̂^†)} <cit.>, and ultimately the minimal three-dimensional basis 
	
    {A}^ T={n̂^2,n̂b̂,n̂b̂^†}={N̂,B̂,B̂^†},

	assumed here, which is of the fourth-degree. We shall later observe that while (<ref>) is necessary to construct the closed Langevin equations, a second-order linearization will be needed to decouple three operators, leaving only the basis (<ref>) in effect. Quite remarkably, however, and in a similar manner, the use of minimal basis (<ref>) turns out to be fairly convenient to construct the optomechanical Langevin equations. This is not only since the Langevin equations take on exactly linear forms, but also eventually the equation for N̂ and B̂^† will decouple. This leaves the whole standard optomechanical interaction exactly solvable through integration of only one linear differential equation in terms of B̂. The main difference between using various choices of higher-order operator bases <cit.> is the noise terms. It turns out that the definition and higher-order operators lead to multiplicative noise inputs, which once known, the problem will be conveniently solvable. Full mathematical treatment of multiplicative noise terms is necessary for description of some various phenomena and this will be discussed in S10 of supplementary information. 
	
	

 §.§ Side-band Inequivalence

	
	Defining Δ_b and Δ_r respectively as the blue and red frequency shifts of sidebands, it is possible to show that these two quantities do not necessarily agree in magnitude, such that Δ_b+Δ_r≠ 0. As shown in S4 of supplementary information, an explicit relation for the side-band inequivalence δΔ=1/2(Δ_r+Δ_b) can be found through series expansion of the eigenvalues of the coefficient matrix from (S19). With some algebra, it is possible to show that for g_0<<Ω  correct to the fourth-order, we get
	
    δΔ/Ω≈(g_0/Ω)^2(n̅+1/2)-2(g_0/Ω)^4(n̅+1/2)(m̅+1/2).

	Here, n̅(Δ) is the intracavity photon population and m̅(Δ) is the coherent phonon population given by
	
    m̅(Δ)≈32 g_0^2 Ω ^2 (γ ^2+γΓ +4 Δ ^2)/(γ ^2+4 Δ ^2) (Γ ^2+4 Ω ^2)^2n̅^2(Δ )=g_0^2ζ(Δ)n̅^2(Δ),

	where Γ is the mechanical decay rate, and γ=κ+Γ is the total optomechanical decay rate with κ being the optical decay rate, as proved in details in S6 of supplementary information using the method of higher-order operators. Also, n̅(Δ) can be found from numerical solution of a third-order algebraic equation (S9). The relationship m̅∝n̅^2 signifies the fact that mechanical oscillations are nonlinearly driven by optical radiation pressure. A typical behavior of this phenomenon is illustrated in Fig. <ref>.
	
	There is a related polaritonic splitting effect <cit.>, as a result of anti-crossing between the optomechanically interacting optical and mechanical resonances generated across either of the mechanical side-bands, amount of which happens to be exactly 2δΔ. This has nothing to do with the side-band asymmetry, which happens to occur on the two opposite sides of the main cavity resonance. It should be mentioned that observation of this phenomenon in superconducting electromechanics <cit.> as well as parametrically actuated nano-string resonators <cit.> can potentially yield the most clear results due to various experimental conditions. In fact, intracavity photon numbers as large as 10^6 and 10^8 and more are attainable respectively in superconducting electromechanics and optically-trapped nano-particle optomechanics. 
	
	A close inspection of a very high-resolution measurement on a side-band resolved microtoroidal disk<cit.> yields a side-band inequivalence of δΔ≈ 2π×(142± 36) Hz, which perfectly complies to (<ref>) if n̅=(5.1± 1.3)× 10^3. Unfortunately, further such a high resolution measurements on deeply side-band resolved optomechanical cavities are not reported elsewhere to the best knowledge of authors. Nevertheless, clear signatures of side-band inequivalence can be easily verified in few other experiments <cit.>. Remarkably, recent measurements on Stokes-Anti-Stokes scattering from multi-layered MoTe_2 exhibits a difference in frequency shift as large as 7% for the five-layered sample <cit.>, which corresponds to 0.88 ± 0.11cm^-1.

	Also, a recent landmark experiment on room-temperature quantum optomechanical correlations <cit.> has reported measurements which coincidentally exhibit a sideband inequivalence up to 4kHz and roughly agree to the approximation δΔ≈ g_0^2n̅/Ω. This issue remains, nevertheless, as an open problem in the context of experimental quantum optomechanics. 
	
	In any experimental attempt to measure this phenomenon, a side-band resolved cavity could be driven on resonance and noise spectra of the two mechanical side bands be measured with extreme precision in a heterodyne setup. Even in case of well-known thermo-optical effects and two-photon dispersion or absorption which cause drifts in the optical resonance and other optomechanical parameters <cit.>, this effect should be still observable in principle. The reason is that the amount of inequivalence is actually independent of the exact pump frequency as long as intracavity photon population does not change significantly. So, it should be sufficient only if the cavity is driven on or close to the optical resonance for the side-bands to be sufficiently different in their frequency shifts.
	
	

 §.§ Higher-order Resonance Shift

	
	The contribution of the off-diagonal terms to the mechanical frequency Ω in the coefficients matrix of optomechanical Langevin equations (S20) of supplementary information, can be ultimately held responsible for the so-called optomechanical spring effect <cit.>. As the result of optomechanical interaction, both of the optical and mechanical resonance frequencies and damping rates undergo shifts. Even at the limit of zero input optical power α=0 and therefore zero cavity photon number n̅=0, it is possible to show that there is a temperature-dependent shift in the mechanical resonance frequency, markedly different from the lattice-expansion dependent effect. This effect is solely due to the optomechanical interaction with virtual cavity photons, which completely vanishes when g_0=0. In close relationship to the shift of resonances, we can also study the optomechanical spring effect with the corrections from higher-order interactions included.
	
	The analysis of spring effect is normally done by consideration of the effective optomechanical force acting upon the damped mechanical oscillator, thus obtaining a shift in squared mechanical frequency δ(Ω^2), whose real and imaginary parts give expressions for δΩ and δΓ. Corrections to these two terms due to higher-order interactions are discussed in <ref>. Here, we demonstrate that the analysis using higher-order operator algebra can recover some important lost information regarding the optical and mechanical resonances when the analysis is done on the linearized basis {A}^T={â,â^†,b̂,b̂^†}. 
	
	To proceed, we consider finding eigenvalues of the matrix M as defined in (S19) of supplementary information. Ignoring all higher-order nonlinear effects beyond the basis {A}^T={â,âb̂,âb̂^†}, we set s=0. This enables us to search for the eigenvalues of the coefficients matrix M as
	
    eig[M]=eig[ 
    	[      iΔ-κ/2        ig_0        ig_0;    i(G+f^+) -i(Ω-Δ)-γ/2           0;   -i(G-f^-)           0  i(Ω+Δ)-γ/2 ]]=i{[ Δ+λ_1+iγ_1; Δ+λ_2+iγ_2; Δ+λ_3+iγ_3 ]}=i{[ Δ+η_1(Δ,T); Δ+η_2(Δ,T); Δ+η_3(Δ,T) ]},

	in which G=g_0n̅, f^±=g_0(m̅+1/2)±1/2g_0, λ_j=[η_j] and γ_j=[η_j] with j=1,2,3 are real valued functions of Δ and bath temperature T. The temperature T determines m̅ while n̅ is a function of Δ as well as input photon rate α. 	In general, the three eigenvalues η_j=λ_j(Δ,T)+iγ_j(Δ,T), j=1,2,3 are expected to be deviate from the three free-running values ψ_1=i1/2κ, ψ_2=-Ω+i1/2γ, and ψ_3=Ω+i1/2γ, as η_j≈ψ_j-Δ because of non-zero g_0. Solving the three equations therefore gives the values of shifted optical and mechanical frequencies and their damping rates compared to the bare values in absence of optomechanical interactions with g_0=0, given by δΩ=-1/2[η_2-η_3]-Ω, δω=-1/2[η_2+η_3], δΓ=[-2η_1+η_2+η_3]-Γ, and δκ=2[η_1]-κ. This  method to calculate the alteration of resonances, does not regard the strength of the optomechanical interaction or any of the damping rates. In contrast, the known methods to analyze this phenomenon normally require g<<κ and Γ+δΓ<<κ <cit.>.
	
 	

 §.§ Corrections to Spring Effect

	
	As shown in S7 of supplementary information, the full expression for corrected spring effect is given as
	Put together combined, we get
	
    δΩ(w,Δ)   =   g_0^2n̅Ω/w[Δ+w/(Δ+w)^2+1/4κ^2+Δ-w/(Δ-w)^2+1/4κ^2]
       +   g_0^2[μ(w)]Ω/w[Δ+w/(Δ+w)^2+1/4κ^2+Δ-w/(Δ-w)^2+1/4κ^2]+g_0^2[μ(w)]Ω/w[κ/(Δ+w)^2+1/4κ^2-κ/(Δ-w)^2+1/4κ^2],
    δΓ(w,Δ)   =   g_0^2n̅Ω/w[κ/(Δ+w)^2+1/4κ^2-κ/(Δ-w)^2+1/4κ^2]
       +   g_0^2[μ(w)]Ω/w[κ/(Δ+w)^2+1/4κ^2-κ/(Δ-w)^2+1/4κ^2]-g_0^2[μ(w)]Ω/w[Δ+w/(Δ+w)^2+1/4κ^2+Δ-w/(Δ-w)^2+1/4κ^2].

	
	Here, the second and third terms on the rights hand sides of both equations are corrections to the spring effect due to the higher-order interactions, resulting from the temperature-dependent expressions 
	
    [μ(w)]   =   w/Ω(m̅+1/2)+1/2, 
    [μ(w)]   =   Γ/2Ω(m̅+1/2).

	The temperature-dependence of (<ref>) causes dependence of the spring effect on temperature as well. The influence of additional terms in (<ref>) due to higher-order interactions can strongly influence any measurement of g_0 through spring effect, as most easily can be observable in the weak coupling limit for Doppler cavities. 

	


	

 §.§ Weak Coupling Limit

	
	In the weakly coupled operation mode and far Doppler regime where g_0<<Ω and κ>>Ω>>Γ hold <cit.>, using (S9) of supplementary information with n̅≈(Δ^2+1/4κ^2)^-1|α|^2, the spring equations are obtained from (<ref>) by setting w=Ω as
	
    δΩ(Ω,Δ)≈ 2Δ g_0^2n̅(Δ)+m̅(Δ)+1/Δ^2+1/4κ^2≈ g_0^2[2Δ|α|^2/(Δ^2+1/4κ^2)^2]+ g_0^4[2Δζ(Δ)|α|^4/(Δ^2+1/4κ^2)^3].

	Here, |α| is photon input rate to the cavity with α being complex drive amplitude, and [μ(Ω)]=m̅+1 and [μ(Ω)]≈ 0 from (<ref>). The importance of this equation is that the optical spring effect is actually proportional to δΩ∝ g_0^2(n̅+m̅)∝ g_0^2n̅(1+g_0^2ζn̅) where ζ(Δ) is already defined in (S30). This shows that if g_0 is to be determined from experimental measurement of the optical spring effect, then the experiment should be done at the lowest optical power possible, otherwise the term m̅∝ g_0^2n̅^2 becomes large and would result in an apparent change in g_0. This fact also can explain why the measured g_0 through optical spring effect using uncorrected standard expressions (S42) is always different from the design value, which could be attributed to the absence of the second term proportional to g_0^4 in the corrected optical spring effect using the higher-order algebra.
	
	The above equation together with the fact that on the far red detuning Δ→+∞ we have n̅(Δ)→ 0, m̅(Δ)→ 0, and δΩ(Ω,Δ)→ 0, provides an alternate approximation for the resonant coherent phonon number m̅(0) at zero-detuning as
	
    m̅(0)   ≈   κ^2/8g_0^2[∂δΩ(Ω,Δ)/∂Δ]_Δ=0-4|α|^2/κ^2-1 ≈32 g_0^2 Q^2_m/Γ^2n̅^2(0) ≈512 g_0^2 Q^2_m/Γ^2κ^4|α|^4,

	where Q_m=Ω/Γ is the mechanical quality factor, and the expression within the brackets can be measured experimentally, and represents the slope of frequency displacement due to the spring effect  versus detuning. The second expression proportional to n̅^2(0) follows (S30) from S6 of supplementary information where an explicit and accurate formula for m̅(Δ) is found. 
	
	Noting |α|^2∝ P_op reveals that while the intracavity photon number is propotional to the optical power as n̅(0)∝ P_op, the coherent phonon population is proportional to the square of the optical power as m̅(0)∝ P_op^2. This implies that the effects of coherent mechanical field gets important only at sufficiently high optical powers, and also marks the fact that in the low optical power limit where linear optomechanics is expected to work well, effects of coherent phonons do not appear. This also explains why this quantity has not been so far noticed in the context of quantum optomechanics. Because it does not show up anywhere in the corresponding fully linearized Langevin equations. 
	
	

§ DISCUSSION

	
	
	
	As shown in S10 of supplementary information, a fairly convenient but approximate solution to the symmetrized spectral density of output optical field due to multiplicative noise is given as
	
    S(ω)=|Y_11(ω)|^2 S_AA(ω)+1/γ^2|[Y_12(ω)+Y_13(ω)]∗a̅(ω)|^2 S_BB(ω) +1/θ^2|[Y_14(ω)∗ab(ω)+Y_15(ω)∗ab^∗(ω)]|^2 S_BB(ω),

	where spectral power densities S_AA and S_BB are already introduced in (S11) of supplementary information and convolutions ∗ take place over the entire frequency axis. In practice it is far easier to use numerical integration, however, this can cause numerical instabilities when |ω-Δ|>1/2Ω. Owing to the fractional polynomial expressions for the elements of scattering matrix elements as well as the multiplicative terms, it is possible to evaluate the integrals exactly using complex residue techniques. Here, we proceed using numerical integration of the convolution integrals. The third term involving the functions ab(ω) and ab^∗(ω) are unnecessary for the 3× 3 second-order formalism, and arise only in the 5× 5 third-order formalism.
	
	In the above equation, every term adds up the contribution from linear, second-order, and third-order optomechanics. These respectively are due to the processes of photon creation-annihilation {â,â^†}, the 1-photon/1-phonon blue {âb̂,â^†b̂^†} and red {âb̂^†,â^†b̂} processes, and the 1-photon/2-phonon second-order blue-like {âb̂^2,â^†b̂^† 2} and red-like {âb̂^† 2,â^†b̂^2} sideband processes. Apparently the Hermitian conjugate operators do not exist in the original 5× 5 higher-order formalism (S70) of supplementary information, since they are completely uncoupled from their Hermitian counterparts. However, calculation of the noise spectral densities necessitates their presence, so that a real-valued and positive definite spectral density has actually already taken care of these conjugate processes. Obviously, the first term contributes to the w=Δ resonance, while the second term contributes to the first-order mechanical side-bands at w=Δ=±Ω. Similarly, the third term constitutes the second-order mechanical sidebands at w=Δ=±2Ω.
	
	It has to be mentioned that the spectral density (<ref>) is not mathematically exact, since the multiplicative operators appearing behind Weiner noise terms, are approximated by their time-averaged frequency-dependent terms (S68) of supplementary information.
	
	As an application example, we simulate the noise spectrum across the red mechanical sideband and optical resonance generated in an optomechanical experiment on the whispering galley mode of an optical micro-toroid, reported in a very remarkable experiment <cit.>. The pump is set around the red mechanical side band for various detuning values ranging from Δ=2π× 60MHz to Δ=2π× 90MHz, and noise spectra are observed. In Fig. <ref>, the simulation results using linearized and higher-order optomechanics are illustrated. Here, the left panel shows the simulations using 3× 3 linear optomechanics (color fills) with the basis {â,b̂,b̂^†} and 4× 4 linearized optomechanics (black lines) using the basis {â,â^†,b̂,b̂^†}. While the linear 4× 4 formalism is expected to be more accurate than the linear 3× 3 formalism, there exists a notable difference between the two approaches. Here, the optomechanical parameters were taken from the same article with some adjustment to resemble the actual experiment <cit.> as  T=65mK, P_op=1.4mW, Ω=2π× 78MHz, g_0=2π× 3.4kHz, Γ=2π× 407kHz, κ=2π× 3.54MHz η=0.5, and λ=775nm. 
	On the right panel of Fig. <ref> the simulations using 3× 3 higher-order optomechanics (color fills) with the basis {â,âb̂,âb̂^†} and 5× 5 higher-order optomechanics (black lines) with the basis {â,âb̂,âb̂^†,âb̂^2,âb̂^† 2} are shown. Agreement between the second-order 3× 3 and third-order 5× 5 formalisms is remarkably good. 
	
	
	
	
	
	It is possible to employ the method of higher-order operators to investigate the dynamic stability of optomechanical systems in the side-band resolved operation limit. A stable optomechanical system can be still perturbed by thermal effects and they appear to be dominant in driving the cavity into instability for Doppler samples. However, for side-band resolved samples, thermal effects are much less pronounced and the major contribution to the instability comes from inherent nonlinear dynamics of the optomechanical interactions. That implies that optomechanical interactions are linearly stable, but they can become nonlinearly unstable at a certain interaction order to be discussed below.
	
	The dynamic stability can be done by inspecting eigenvalues of the coefficients matrix [M]. If the real part of at least one of the eigenvalues is positive, then the system is unstable and its response to any perturbation grows indefinitely in time. The linear formalisms of optomechanics fails to describe this phenomenon, since they always yield constant eigenvalues. Even the second-order higher-operator method with 3× 3 formalism, which describes the nonlinear 1-photon/1-phonon processes, fails to reproduce the correct expected stability phases. This only can be understood by employing at least the third-order 5× 5 operator method, which includes the nonlinear 1-photon/2-phonon processes. Hence, surprisingly enough, it is the 1-photon/2-phonon process and beyond, which contributes to the unstability of an optomechanical system.
	
	To illustrate this, we calculate the stability phases of the side-band resolved system investigated in S10 of supplementary information. Illustrated in Fig. <ref>, the stable and unstable regions of this systems across blue Δ<0 and red Δ>0 detunings versus input optical power P_op are illustrated. The v-shaped region in violet color, maps the unstable phase, while the red and blue colors correspond to the stable operation phases. 
	
	We also have calculated the linear stability from 4× 4 full linear formalism, which appears on the right of Fig. <ref>. Not surprisingly, the linear and nonlinear stability diagrams remarkably are different. The linear instability starts at moderate resonant pump powers, while it starts rapidly growing exactly over the blue detuning at a slightly higher power. 
	
	Firstly, it can be seen that across almost the entire domains of linear stability, the system is also nonlinearly stable. Secondly, by observation of the nonlinear stability in the left and linear stability on the right, it can be seen that in most of the domain of linear stability, the system is already nonlinearly stable. Hence, any attempt to drive the system within the region of linearly unstable but nonlinearly stable, ultimately results in significant growth of mechanical amplitude and therefore side-bands. Any further increase in the amplitude of side-bands become limited due to nonlinear stability. Hence, four possible stability scenarios could be expected:
	
		
  * Linearly and Nonlinearly Stable: The intersect of the linear and nonlinear stability domains, marks a shared domain of unconditionally stable optomechanical interaction. Unless the system is influenced by thermal or other nonideal effects, the stability is always guaranteed.
		
  * Linearly Unstable, but Nonlinearly Stable: By inspection, an optomechanical system can be linearly unstable while nonlinearly stable. This corresponds to the domains where any attempt to drive the system in these regions causes immediate but limited growth in the amplitude of mechanical oscillations.
		
  * Linearly and Nonlinearly Unstable: Under this scenario, the optomechanical cavity is always unstable regardless of the other nonideal effects. This happens only at remarkably high drive powers.
		
  * Linear Stable, but Nonlinearly Unstable: The unlikely and surprising case of linear stability and nonlinear unstability is also possible according to the stability maps at some portions of non-resonant high drive powers. This strange behavior corresponds to the case when the system remains stable only at infinitesimal optical powers. Any fluctuation beyond tiny amplitudes shall drive the system into unstable growing and large amplitudes.
	
	
	There is a threshold power P_th at which instabilites start to appear. For the side-band resolved case, this happens at P_th=5mW on the blue side. As it expected and in agreement to experimental observations, the unstable domain mostly covers the blue domain with Δ<0. However, at higher optical powers than P_th, instability phase can diffuse well into the red detunings as well Δ>0. For the cavity under consideration, this happens at much higher optical power of P=14P_th=70mW. Therefore, the general impression that instability always occurs on the blue side at every detuning above a certain threshold power is not correct. Interestingly, the boundary separating the dynamically stable and dynamically unstable phases for this side-band resolved sample with Ω>>κ>>Γ can be well estimated using
	
    n̅(Δ)>n_cr,

	in which n_cr is a critical cavity photon number. Extensive numerical tests for cavities within deep side-band resolved Ω>>κ, deep Doppler Ω>>κ, and intermediate Ω∼κ regimes reveal the existence of such a critical intracavity photon number limit n_cr, beyond which dynamical instability takes over. However, for a cavity in deep side-band resolved regime, it can be estimated in a phenomenological way, and is roughly given by
	
    n_cr∼4/𝒞_0(Ω/κ)^2.

	
	
	Here, 𝒞_0 is the single-photon cooperativity given by
	
    𝒞_0=g_0^2/κΓ.

	If the cavity is not side-band resolved, (<ref>) cannot be used, but numerical computations can still yield the limiting number n_cr. 
	
	For Doppler cavities, no such dynamic instability can be observed, and therefore thermal effects should dominate over dynamical effects in driving a Doppler cavity toward instability. Meanwhile, it is the nonlinear optomechanical dynamics which seems to be dominant in driving a side-band resolved cavity into instability. As a result, the existence of such a critical maximum intracavity photon number is not related to thermal effects, but rather to the nonlinear stability. 
	
	

§ METHODS

	
	For an extensive description of theoretical methods, refer to the Supplementary Information.
	
	

§ CONCLUSIONS

	
	A new analytical method was shown to solve the standard optomechanical interaction with cubic nonlinearity interaction, based on the higher-order operators. It was demonstrated that not only the higher-order operator method can reproduce the linear optomechanics, but also it can predict and provide estimates to unnoticed effects such as a new type of symmetry breaking in frequency, here referred to as side-band inequivalence, and yield new explicit expressions for quantities such as the coherent phonon population and higher-order spring effect. Corrections to the standard spring effect due to higher-order interactions have been found, and it has been shown that such corrections arise mainly because of the coherent phonons and can significantly influence measurement of single-photon optomechanical interaction rate through spring effect. A minimal basis has been defined which allows exact and explicit solution to standard nonlinear optomechanics, using the method of higher-order operators. This method can be finally used to investigate the dynamic nonlinear stability of optomechanical systems, and it has been demonstrated that at least the third-order nonlinear processes are prerequisite for occurrence of dynamic instability. We have shown that there is a reasonable correspondence between the onset of nonlinear dynamic instability and a critical intracavity photon number limit, which remains independent of thermal effects. 
	
	

§ ADDITIONAL INFORMATION

	
	The author declares no competing financial and/or non-financial interests in this work.
	
	

§ ACKNOWLEDGEMENT

	
	Discussions with Georg Arnold at IST Austria and group members of the laboratory of Quantum Foundations and Quantum Information on the Nano- and Microscale are appreciated. This paper is dedicated to the celebrated artist, Anastasia Huppmann.
	
	99
		Fokker1
		Risken, H. The Fokker-Planck Equation: Methods of Solution and Applications (Springer: Berlin, 1996).
		Fokker2
		Shapovalov, A. V., Rezaev, R. O. & Trifonov, A. Y. Symmetry operators for the Fokker–Plank–Kolmogorov equation with nonlocal quadratic nonlinearity. Sigma 3, 005 (2007).
		Fokker3
		Pavliotis, G. A. Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations (Springer: New York, 2014).
		Fokker4
		Carmichael, H. J. Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations (Springer: Berlin, 2002).
		Fokker5
		Colmenares, P. J. Fokker-Planck equation of the reduced Wigner function associated to an Ohmic quantum Langevin dynamics. Phys. Rev. E 97, 052126 (2018).
		
		Moments
		Kim, K. I. Higher order bias correcting moment equation for M-estimation and its higher order efficiency. Econometrics 4, 48 (2016).
		
		Master1 
		Ludwig, M., Kubala, B. & Marquardt, F. The optomechanical instability in the quantum regime. New J. Phys. 10,  095013 (2008).
		Master2 
		Hamerly, R. & Mabuchi, H. Quantum noise of free-carrier dispersion in semiconductor optical cavities. Phys. Rev. A 92, 023819 (2015).
		Master3 
		Nunnenkamp, A., Børkje, K. & Girvin, S. M. Single-photon optomechanics. Phys. Rev. Lett. 107, 063602 (2011).
		Master4 
		Rips, S., Kiffner, M., Wilson-Rae, I. & Hartmannnew, M. J. Steady-state negative Wigner functions of nonlinear nanomechanical oscillators. New J. Phys. 14, 023042 (2012).
		Master5
		Boyanovsky, D. & Jasnow, D. Heisenberg-Langevin vs. quantum master equation. Phys. Rev. A 96, 062108 (2017).
		
		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).
		
		Law
		Law, C. K. Interaction between a moving mirror and radiation pressure: A Hamiltonian formulation. Phys. Rev. A 51, 2537 (1995).
		Macri
		Macrì, V., et al. Non-perturbative dynamical Casimir effect in optomechanical systems: Vacuum Casimir-Rabi splittings. Phys. Rev. X 8, 011031 (2018).
		Kip1
		Kippenberg, T. J. & Vahala, K. J. Cavity optomechanics. Science 321,  1172 (2008).
		Aspel1
		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).
		Opto1
		Chang, Y., Shi, T., Liu, Y.-X., Sun, C. P. & Nori, F. Multistability of electromagnetically induced transparency in atom-assisted optomechanical cavities. Phys. Rev. A 83, 063826 (2011). 
		Opto2
		Jing, H., et al. Optomechanically-induced transparency in parity-time-symmetric microresonators. Sci. Rep. 5, 9663 (2015). 
		Opto3
		Cirio, M., Debnath, K., Lambert, N. & Nori, F. Amplified optomechanical transduction of virtual radiation pressure. Phys. Rev. Lett. 119, 053601 (2017). 
		Polariton
		Ian, H., Gong, Z. R., Liu, X.-Y., Sun, C. P. & Nori, F. Cavity optomechanical coupling assisted by an atomic gas. Phys. Rev. A 78, 013824 (2008). 
		Multi1
		Jing, H., Özdemir, Ş. K., Lü, H. & Nori, F. High-order exceptional points in optomechanics. Sci. Rep. 7, 3386 (2017). 
		Multi2
		Liao, J.-Q. & Nori, F. Spectrometric reconstruction of mechanical-motional states in optomechanics. Phys. Rev. A 90, 023851 (2014). 
		Multi3
		Liao, J.-Q., Law, C. K., Kuang, L.-M. & Nori, F. Enhancement of mechanical effects of single photons in modulated two-mode optomechanics. Phys. Rev. A 92, 013822 (2015). 
		Multi4
		Liao, J.-Q., Wu, Q.-Q. & Nori, F. Entangling two macroscopic mechanical mirrors in a two-cavity optomechanical system. Phys. Rev. A 89, 014302 (2014). 
		Multi5
		Liu, Y.-L., et al. Controllable optical response by modifying the gain and loss of a mechanical resonator and cavity mode in an optomechanical system. Phys. Rev. A 95, 013843 (2017).
		Kerr1
		Gong, Z. R., Ian, H., Liu, Y. X., Sun, C. P. & Nori, F. Effective Hamiltonian approach to the Kerr nonlinearity in an optomechanical system. Phys. Rev. A 80, 065801 (2009).
		Kerr2 
		Lu, X. Y., Zhang, W.-M., Ashhab, S., Wu, Y. & Nori, F. Quantum-criticality-induced strong Kerr nonlinearities in optomechanical systems. Sci. Rep. 3, 2943 (2013). 
		Kerr3
		Wang, H., Gu, X., Liu, Y.-X., Miranowicz, A. & Nori, F. Tunable photon blockade in a hybrid system consisting of an optomechanical device coupled to a two-level system. Phys. Rev. A 92, 033806 (2015).  
		Kerr4
		Johansson, J. R., Johansson, G. & Nori, F. Optomechanical-like coupling between superconducting resonators. Phys. Rev. A 90, 053833 (2014). 
		Kerr5
		Wang, H., Gu, X., Liu, Y.-X., Miranowicz, A. & Nori, F. Optomechanical analog of two-color electromagnetically-induced transparency: Photon transmission through an optomechanical device with a two-level system. Phys. Rev. A 90, 023817 (2014). 
		Kerr6
		Monifi, F.,  et al. Optomechanically induced stochastic resonance and chaos transfer between optical fields. Nat. Photon. 10, 399 (2016).
		
		Duffing
		Lü, X.-Y., Liao, J.-Q., Tian, L. & Nori, F. Steady-state mechanical squeezing in an optomechanical system via Duffing nonlinearity. Phys. Rev. A 91, 013834 (2015). 
		X2
		Lü, X.-Y., et al. Squeezed optomechanics with phase-matched amplification and dissipation. Phys. Rev. Lett. 114,  093602 (2015). 
		
		Girvin
		Børkje, K., Nunnenkamp, A., Teufel, J. D. & Girvin, S. M. Signatures of nonlinear cavity optomechanics in the weak coupling regime. Phys. Rev. Lett. 111, 053603 (2012).
		
		Lemonde2
		Lemonde, M.-A., Didier, N. & Clerk, A. A. Nonlinear interaction effects in a strongly driven optomechanical cavity. Phys. Rev. Lett. 111, 053602 (2013).
		
		OMIT1
		Weis, S.,  et al. Optomechanically induced transparency. Science 330, 1520 (2010).
		
		OMIT2
		Teufel, J. D.,  et al.  Circuit cavity electromechanics in the strong-coupling regime. Nature 471, 204 (2011).
		Polaron
		Verhagen, E., Deléglise, S., Weis, S., Schliesser, A. & Kippenberg, T. J. Quantum-coherent coupling of a mechanical oscillator to an optical cavity mode. Nature 482, 63 (2012).
		Rabl
		Rabl, P. Photon blockade effect in optomechanical systems. Phys. Rev. Lett. 107, 063601 (2015).
		
		Lemonde1
		Lemonde, M.-A., Didier, N. & Clerk, A. A. Enhanced nonlinear interactions in quantum optomechanics via mechanical amplification. Nat. Commun. 7, 11338 (2016).
		
		Lemonde3
		Lemonde, M.-A. & Clerk, A. A. Real photons from vacuum fluctuations in optomechanics: The role of polariton interactions. Phys. Rev. A. 91, 033836 (2015).
		
		Self1
		Kippenberg, T. J., Rokhsari, H., Carmon, T., Scherer, A. & Vahala, K. J. Analysis of radiation-pressure induced mechanical oscillation of an optical microcavity. Phys. Rev. Lett. 95, 033901 (2005).
		Self2
		Marquardt, F., Harris, J. G. E. & Girvin, S. M. Dynamical multistability induced by radiation pressure in high-finesse micromechanical optical cavities. Phys. Rev. Lett. 96, 103901 (2006).
		Self3
		Vahala, K.,  et al.  A phonon laser. Nat. Phys. 5, 682 (2009).
		Self4
		Bagheri, M., Poot, M., Li, M., Pernice, W. P. H. & Tang, H. X. Dynamic manipulation of nanomechanical resonators in the high-amplitude regime and non-volatile mechanical memory operation. Nat. Nanotechnol. 6,  726 (2011).
		Self5
		Krause, A. G., et al. Nonlinear radiation pressure dynamics in an optomechanical crystal. Phys. Rev. Lett. 115, 233601 (2015).
		
		Chaos1
		Carmon, T., Cross, M. & Vahala, K. Chaotic quivering of micron-scaled on-chip resonators excited by centrifugal optical pressure. Phys. Rev. Lett. 98, 167203 (2007).
		Chaos2
		Navarro-Urrios, D., et al. Nonlinear dynamics and chaos in an optomechanical beam. Nat. Commun. 8, 14965 (2017).
		
		
		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).
		Quad2
		Sankey, J. C., Yang, C., Zwickl, B. M., Jayich, A. M. & Harris, J. G. E. Strong and tunable nonlinear optomechanical coupling in a low-loss system. Nat. Phys. 6, 707 (2010).
		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
		Lei, C. U., et al. Quantum nondemolition measurement of a quantum squeezed state beyond the 3 dB limit. Phys. Rev. Lett. 117, 100801 (2016).
		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
		Vanner, M. R. Selective linear or quadratic optomechanical coupling via measurement. Phys. Rev. X 1, 021011 (2011).
		Quad9
		Doolin, C.,  et al.  Nonlinear optomechanics in the stationary regime. Phys. Rev. A 89, 053838 (2014).
		Quad10
		Kim, E.-J., Johansson, J. R. & Nori, F. Circuit analog of quadratic optomechanics.  Phys. Rev. A 91, 033835 (2015). 
		Quad11
		Liao, J.-Q. & Nori, F. Single-photon quadratic optomechanics. Sci. Rep. 4, 6302 (2014). 
		
		
		Bruschi
		Bruschi, D. E. & Xuereb, A. Mechano-optics: An optomechanical quantum simulator. New J. Phys. 20, 065004 (2018).
		
		Paper1
		Khorasani, S. Higher-order interactions in quantum optomechanics: Revisiting the theoretical foundations. Appl. Sci. 7, 656 (2017).
		Paper2
		Khorasani, S. Higher-order interactions in quantum optomechanics: Analytical solution of nonlinearity. Photonics 4, 48 (2017).
		Paper4
		Khorasani, S. Higher-order interactions in quantum optomechanics: Analysis of quadratic terms. arxiv, 1712.9615 (2017).
		Chapter
		Khorasani, S. Momentum-field interactions beyond standard quadratic optomechanics. arxiv, 1801.1800 (2018).
		Paper5
		Khorasani, S. Solution of cross-Kerr interaction combined with parametric amplification. arxiv, 1804.7053 (2018).
						
		Noise
		Rubin, K. J., Pruessner, G. & Pavliotis, G. A. Mapping multiplicative to additive noise. J. Phys. A: Math. Theor. 47, 195001 (2014).
		
		Stokes1
		Neuhaus, R., Sellars, M. J., Bingham, S. J. & Suter, D. Breaking the Stokes–anti-Stokes symmetry in Raman heterodyne detection of magnetic-resonance transitions. Phys. Rev. A 58, 4961 (1998).
		Stokes2
		Goldstein, T., et al. J. Raman scattering and anomalous Stokes–anti-Stokes ratio in MoTe_2 atomic layers. Sci. Rep. 6, 28024 (2016).
		Stokes3
		Kittlaus, E. A., Otterstrom, N. T. & Rakich, P. T. On-chip inter-modal Brillouin scattering. Nat. Commun. 8,  15819 (2017).
		
		Liu
		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).
		
		Mika
		Ockeloen-Korppi, C. F., et al.  Entangled massive mechanical oscillators. Nature 556, 478 (2018).
		
		Weig1
		Gajo, K., Schüz, S. &  Weig, E. M. Strong 4-mode coupling of nanomechanical string resonators. Appl. Phys. Lett. 111, 133109 (2017).
		Weig2
		Seitner, M. J., Abdi, M., Ridolfo, A., Hartmann, M. J. &  Weig, E. M. Parametric oscillation, frequency mixing, and injection locking of strongly coupled nanomechanical resonator modes. Phys. Rev. Lett. 118, 254301 (2017).
		
		Kip3
		Schliesser, A., Rivière, R., Anetsberger, G., Arcizet, O. & Kippenberg, T. J. Resolved-sideband cooling of a micromechanical oscillator. Nature Phys. 4, 415 (2008).
		
		Teufel
		Clark, J. B., Lecocq, F., Simmonds, R. W., Aumentado, J. & Teufel, J. D. Sideband cooling beyond the quantum backaction limit with squeezed light. Nature 541, 191 (2017).
		
		Vivishek
		Sudhir, V.,  et al.  Quantum correlations of light from a room-temperature mechanical oscillator. Phys. Rev. X 7, 031055 (2017).
		
		Painter
		Barclay, P. E., Srinivasan, K. & Painter, O. Nonlinear response of silicon photonic crystal microresonators excited via an integrated waveguide and fiber taper. Opt. Express 13, 801 (2005).
		
		Spring1
		Teufel, J. D., Harlow, J. W., Regal, C. A. &  Lehnert, K. W. Dynamical backaction of microwave fields on a nanomechanical oscillator. Phys. Rev. Lett. 101, 197203 (2008).
		Spring2
		Eichenfield, M., Camacho, R., Chan, J., Vahala, K. J. &  Painter, O. A picogram- and nanometre-scale photonic-crystal optomechanical cavity. Nature 459, 550 (2009).
		Spring3
		Safavi-Naeini, A. H.,  et al. Squeezed light from a silicon micromechanical resonator. Nature 500,  185 (2013).
		Spring4
		Deotare, P. B.,  et al.  All optical reconfiuration of optomechanical fiters. Nat. Commun. 3, 846 (2011).
		Spring5
		Sarabalis, C. J., Dahmani, Y. D., Patel, R. N., Hill, J. T. &  Safavi-Naeini, A. H.. Release-free silicon-on-insulator cavity optomechanics. Optica 4, 1147 (2017).
		
		
		Transduction
		Leijssen, R. & Verhagen, E. Strong optomechanical interactions in a sliced photonic crystal nanobeam. Sci. Rep. 5, 15974 (2015).
		
		Verhagen
		Leijssen, R., La Gala, G. R., Freisem, L., Muhonen, J. T. & Verhagen, E. Nonlinear cavity optomechanics with nanomechanical thermal fluctuations. Nat. Commun. 8, 16024 (2017).
		
		
	
	/S21


§ SUPPLEMENTARY INFORMATION: THEORETICAL METHODS



	Method of Higher-order Operators for Quantum Optomechanics
	
	Sina Khorasani









 








§ OPTOMECHANICAL HAMILTONIAN


The Langevin equations for the Hamiltonian ℍ_OM=ħΩm̂-ħΔn̂-ħ g_0n̂(b̂+b̂^†) with the basis {A}^T={â,b̂,âb̂,âb̂^†,n̂,ĉ} are given exactly by 

    [ [ [               iΔ-κ/2                    0                 ig_0                 ig_0                    0                    0;                    0            -(iΩ+Γ/2)                    0                    0                 ig_0                    0;        ig_0(m̂+n̂+1)                    0   -i(Ω-Δ-g_0 b̂)-γ/2                    0                    0                    0;          ig_0(m̂-n̂)                    0                    0   i(Ω+Δ+g_0b̂^†)-γ/2                    0                    0;                    0                    0                    0                    0                   -κ                    0;                    0                    0                ig_0â                ig_0â                    0 2i[Δ+g_0(b̂+b̂^†)]-κ ] ]]
    ×{[     â;    b̂;   âb̂; âb̂^†;    n̂;     ĉ ]}
    -{[       √(κ)â_in;      √(Γ)b̂_in;   √(γ)(âb̂)_in; √(γ)(âb̂^†)_in;     √(2κ)n̂_in;      √(2κ)ĉ_in ]}=d/dt{[     â;    b̂;   âb̂; âb̂^†;    n̂;     ĉ ]},

in which γ =κ +Γ. We have set x̂=â in all equations except the second and third where both of the bath operators x̂=â and x̂=b̂ are taken separately to construct the noise terms, √(2)n̂_in=â^†â_in+â^†_inâ, (âb̂^†)_in=â_inb̂^†+âb̂_in^†, (âb̂)_in=â_inb̂+âb̂_in and √(2)ĉ_in=ââ_in. 

The system (<ref>) is still nonlinear and non-integrable because of the dependence of the coefficients matrix on the operators. But it can be simplified by first noting that from the fifth equation we could expect any disturbance in n̂ would decay as δn̂(t)∼exp(-κ t) on time scales smaller than κ^-1. This can be further approximated as n̂∼n̅ at steady input. Similar argument goes with δm̂∼exp(-Γ t) in response to a disturbance on time scales smaller than Γ^-1, which enables us to make the approximate replacement  m̂∼m̅ at equilibrium. 

For the phononic mechanical operators b̂ and b̂^† appearing within the brackets, approximate decays δb̂(t)∼exp[-(iΩ+1/2Γ) t] and δb̂^†(t)∼exp[(iΩ-1/2Γ) t] in response to disturbances hold, making the coefficients matrix time-dependent. But these can be nevertheless dropped in whole if we notice that g_0b̅<<Ω which is the normal experimental condition of weakly-coupling in optomechanics. Otherwise, they can replaced by constant amplitudes b̅ and b̅^∗ given below in (<ref>) on sufficiently longer time scales than Γ^-1 for strongly-coupled systems. 

Such types of approximations are in fact quite highly in use within the context of continuous wave standard optomechanics. Therefore, once the steady state solution to (<ref>) around the equilibrium values due to optical drive ⟨â_in⟩=α is sought, the coefficients matrix can be kept time-independent, keeping only the fluctuations of input terms as the only source. The case of time dependent drive α=α(t) for pulsed experiments shall be discussed later in the article.

Having said that, all the operators n̂, m̂, b̂, and b̂^† in the coefficients matrix can be replaced by their respective average values to proceed with the second-order accurate optomechanical system of equations as


    [ [ [         iΔ-κ/2              0           ig_0           ig_0              0              0;              0      -(iΩ+Γ/2)              0              0           ig_0              0;           iL^+              0  -i(Ω-Δ-s)-γ/2              0              0              0;           iL^-              0              0 i(Ω+Δ+s^∗)-γ/2              0              0;              0              0              0              0             -κ              0;              0              0             ig             ig              0   2i(Δ+2[s])-κ ] ]]
    ×{[     â;    b̂;   âb̂; âb̂^†;    n̂;     ĉ ]}
    -{[       √(κ)â_in;      √(Γ)b̂_in;   √(γ)(âb̂)_in; √(γ)(âb̂^†)_in;     √(2κ)n̂_in;      √(2κ)ĉ_in ]}=d/dt{[     â;    b̂;   âb̂; âb̂^†;    n̂;     ĉ ]},


in which g=g_0√(n̅), s=g_0b̅ with [b̅]=x̅/2x_ zp and x_ zp being the zero-point displacement,  L^+=g_0(m̅+n̅+1), and L^-=g_0(m̅-n̅). These can be further approximated by L^±≈± g_0n̅=± F under normal experimental conditions of an ultracold cavity with sufficiently high pumping. The average mirror displacement x̅ is due to the average radiation pressure. The fact that L^+≠-L^- provides the quantum mechanical asymmetry between blue and red sidebands. 

It is easy to verify that this way of linearization decouples the state operators and reduces the space into a 3-dimensional one spanned by {A}^ T={â,âb̂,âb̂^†}. This will be discussed in further details later.

In the absence of red-side-band optical cooling tone as well as any other interaction, the average population value is m̅=1/[exp(ħΩ/k_BT)-1], while n̅ can be obtained from the steady state solution of the first row by replacements of input noise term √(κ)â_in→α+√(κ)â_in. Here, α is the input photon flux originally due to an undisplayed resonant drive term ℍ_d=ħ(αâ+α^∗â^†) added to the Hamiltonian ℍ_OM. Furthermore, α has some non-zero phase taken from the cavity population n̅ away. Now that the drive term ℍ_d has been dropped from ℍ_OM, and â_in now only contains the fluctuations with zero-average ⟨â_in⟩=0. 

As it will be shown later, the quantity m̅ here being referred to as the coherent phonon population, can enter the optomechanical interaction processes due to higher-order effects, where its value normally needs to be fitted for cavities in the Doppler regime. Hence, the coherent phonon population m̅ is independent of the simple thermal equilibrium value m, and actually represents those number of phonons who take part in the optomechanical interaction. We will observe that the negative detunings with the blue process can actually lead to a relatively constant phonon population, whereas on the red detunings it starts to decrease with the detuning.

Defining  K=n̅κ we may use the substitutions for the noise and input terms as

    √(γ)(âb̂)_in   →   √(Γn̅)b̂_in+√(κ)b̅â_in+b̅α,
    √(γ)(âb̂^†)_in   →   √(Γn̅)b̂_in^†+√(κ)b̅^∗â_in+b̅^∗α,
    √(κ)n̂_in   →   √(K)â_in+√(K)â_in^†+2√(n̅)[α],
    √(κ)ĉ_ in   →   √(K)â_in+√(n̅)α.

These substitutions follow the fact that terms such as ââ_ in which contain the interaction of a time-dependent operator â(t) and a purely white Weiner noise process with zero average ⟨â_ in⟩=0, can be fairly well approximated by noting first that â(t)∼a̅exp(iΔ t) around the equilibrium, and then noting that shifting the noise process â_ in in frequency to the amount of Δ has essentially no effect by definition. Hence, the sinusoidal time dependence exp(iΔ t) is irrelevant and can be dropped. Similar arguments hold for the phononic operator b̂(t)∼b̅exp(-iΩ t) and their Hermitian adjoints interacting with a white noise term with uniform spectrum.

This allows us to ultimately rewrite the Langevin equations (<ref>) as 

    [ [ [         iΔ-κ/2              0           ig_0           ig_0              0              0;              0      -(iΩ+Γ/2)              0              0           ig_0              0;           iL^+              0  -i(Ω-Δ-s)-γ/2              0              0              0;           iL^-              0              0 i(Ω+Δ+s^∗)-γ/2              0              0;              0              0              0              0             -κ              0;              0              0             ig             ig              0   2i(Δ+2[s])-κ ] ]]
    ×{[     â;    b̂;   âb̂; âb̂^†;    n̂;     ĉ ]}
    -[ 
    [     √(κ)        0        0        0;        0        0     √(Γ)        0;   √(κ)b̅        0   √(Γn̅)        0; √(κ)b̅^∗        0        0   √(Γn̅);     √(K)     √(K)        0        0;     √(K)        0        0        0;          ]]
    {[    â_in;  â_in^†;   b̂_in; b̂_in^† ]}-
    [ 
    [     1     0;     0     0;    b̅     0;  b̅^∗     0; √(n̅) √(n̅); √(n̅)     0;       ]]
    {[   α; α^∗ ]}
    =d/dt{[     â;    b̂;   âb̂; âb̂^†;    n̂;     ĉ ]}.

The second term on the right is the noise fluctuations due to the optical and mechanical fields with zero average ⟨â_in⟩=⟨b̂_in⟩=0, and the last term in the above is proportional to the input photon flux |α|=√(ϵκ P/ħω) where P is the incident radiation power and ϵ is the coupling efficiency. As it will be mentioned briefly later, the average values n̅ and x̅ have to be solved by setting d/dt=0 on the left and taking average values, which eliminates the noise fluctuations, causing the replacements â→√(n̅), b̂→b̅, âb̂→b̅√(n̅), âb̂^†→b̅^∗√(n̅), n̂→n̅, and ĉ→n̅/2.

Hence, the average values a̅=√(n̅) and b̅ get nonlinearly coupled to the input flux α through the system of algebraic relations as


    [ [ [                iΔ-κ/2                     0                  ig_0                  ig_0                     0                     0;                     0             -(iΩ+Γ/2)                     0                     0                  ig_0                     0;                  iL^+                     0     -i(Ω-Δ-g_0b̅)-γ/2                     0                     0                     0;                  iL^-                     0                     0    i(Ω+Δ+g_0b̅^∗)-γ/2                     0                     0;                     0                     0                     0                     0                    -κ                     0;                     0                     0                    ig                    ig                     0 i[Δ+g_0(b̅+b̅^∗)]-κ/2 ] ]]
    ×{[     a̅;     b̅;   a̅b̅; a̅b̅^∗;   a̅^2;   a̅^2 ]}
    ={[    1    0;    0    0;   b̅    0; b̅^∗    0;   a̅   a̅;   a̅    0;      ]}{[   α; α^∗ ]}.


With a given input photon flux |α|, this system can be now solved to obtain the phase ∠α in such a way that ∠a̅=0. Then a̅ and b̅ can be obtained in an algebraic manner. This sets up a system of equations in terms of the total of four unknowns ∠α, a̅=√(n̅), b̅, and b̅^∗. 

In the above system, the second equation is independent of α, while together the fifth they yield 

    b̅   =   ig_0/iΩ+1/2Γa̅^2,
    a̅   =   -1/κ(α+α^∗).

This also already solves b̅^∗ in terms of a̅. Plugging in the results into the first equation leads to the third-order algebraic equation which can be now solved. Doing this and some algebraic manipulation gives the equation

    ig_0^22Ω/Ω^2+1/4Γ^2a̅^3+(iΔ-κ/2)a̅=α.

This equation in general is expected to yield only real-valued a̅. Separating the real and imaginary parts gives

    [α]   =   -κa̅/2
    [α]   =   g_0^22Ω/Ω^2+1/4Γ^2a̅^3+Δa̅.

The first of these is the same as the second of (<ref>). The above two equations can be now iteratively solved to yield ∠α and a̅ for a given |α|. One may also discard ∠α by combining the above two, resulting in 

    |α|^2=[κ^2/4+(2g_0^2Ω/Ω^2+1/4Γ^2n̅+Δ)^2]n̅.

Only real and positive-valued roots of (<ref>) for n̅ are acceptable. Sufficiently large blue-detuning with Δ<Δ_b<0 causes the well-known bistability. It is easy to find the negative blue detuning Δ_b<0 at which bistability starts to appear, by looking for the only negative real root of the cubic equation

    -Δ_b(Δ_b^2+9/4κ^2)=27g_0^2Ω/Ω^2+1/4Γ^2|α|^2.


These two noise terms we assume have the flat shot-noise uncorrelated spectral power densities 

    S_AA(ω)=1/2, 
    
    S_BB(ω)=m+1/2,

which are identical on both positive and negative frequencies. The ultimate difference of noise power spectral densities will be later maintained by the asymmetry caused by L^+ -L^-=2F+g_0 > 0. Here, m=1/[exp(ħΩ/k_BT)-1] is the population of incoherent phonons under thermal equilibrium, which contribute to the random fluctuations of thermal noise. This quantity is not to be mistaken with m̅ which here denotes the population of coherent phonons, contributing coherently to the optomechanical interaction, and are driven by the optical radiation pressure. This shall be discussed later in <ref>.



 §.§ Perturbative Solution


At this moment, the system of equations (<ref>) can be perturbed around equilibrium values found above. This procedure and taking a Fourier transform gives out the solution. Let us define first

    M=[ [ [         iΔ-κ/2              0           ig_0           ig_0              0              0;              0      -(iΩ+Γ/2)              0              0           ig_0              0;           iL^+              0  -i(Ω-Δ-s)-γ/2              0              0              0;           iL^-              0              0 i(Ω+Δ+s^∗)-γ/2              0              0;              0              0              0              0             -κ              0;              0              0             ig             ig              0   2i(Δ+2[s])-κ ] ]],

as well as

    {δ A(ω)}^T   =   {δâ(ω), 
    δb̂(ω),  
    δ(âb̂)(ω),
    δ(âb̂^†)(ω), 
    δn̂(ω), 
    δĉ(ω) 
    },
    {A_in(ω)}^T   =   {â_in(ω) ,
    â_in^†(ω) , 
    b̂_in(ω) ,
    b̂_in^†(ω) 
    },
    [√(Γ)]   =   [ 
    [     √(κ)        0        0        0;        0        0     √(Γ)        0;   √(κ)b̅        0   √(Γn̅)        0; √(κ)b̅^∗        0        0   √(Γn̅);     √(K)     √(K)        0        0;     √(K)        0        0        0;          ]].


Then, taking I_j as the j× j identity matrix, we get

    {A_out(ω)}   =   {A_in(ω)}-[√(Γ)]^T{δ A(ω)}=Y(ω){A_in(ω)}, 
    {δ A(ω)}   =   Z(ω){A_in(ω)}, 
    Z(ω)   =   [M-iωI_6]^-1[√(Γ)], 
    Y(ω)   =   I_4-[√(Γ)]^TZ(ω).

Here, Y(ω) is the scattering matrix connecting the input and output ports. Now, the spectral density of reflected light from the cavity can be found using (<ref>) by the expression

    S(ω)=[|Y_11(ω)|^2+|Y_12(ω)|^2]S_AA(ω)+[|Y_13(ω)|^2+|Y_14(ω)|^2]S_BB(ω),

as long as the noise processes of â_in and b̂_in have zero cross-correlation <cit.>.



§ LINEARIZED OPTOMECHANICS


It is fairly easy to see that the system of equations (<ref>) when simplified and rewritten for the basis {â,b̂,b̂^†} reproduces the widely used linearized optomechanical equations <cit.>. To demonstrate this, we ignore the perturbation matrix δN, as well as b̅/√(n̅) in the noise terms, and then employ the substitutions 

    âb̂   →   exp[(iΔ-1/2κ)t]a̅b̂=exp[(iΔ-1/2κ)t]√(n̅)b̂,
    âb̂^†   →   exp[(iΔ-1/2κ)t]a̅b̂^†=exp[(iΔ-1/2κ)t]√(n̅)b̂^†.

This will immediately result in rewriting (<ref>) as

    d/dt{[    δâ;   δb̂; δb̂^† ]}=[ 
    [  iΔ-κ/2    ig_0    ig_0;       0 -iΩ-Γ/2       0;       0       0  iΩ-Γ/2 ]]
    {[    δâ;   δb̂; δb̂^† ]}+[ 
    [ √(κ)    0    0;    0 √(Γ)    0;    0    0 √(Γ) ]]{[    â_in;   b̂_in; b̂_in^† ]},

which is nothing but exactly the linearized state equations of optomechanics. Hence, the method of higher-order operators <cit.> is mathematically able to reproduce the less approximate linearized optomechanics.



§ PULSED DRIVE


Under the situation of pulsed drive, one may assume the input photon rate α(t) to be a function of time. If the input drive varies on a time-scale or longer than the mechanical period with |dα(t)/dt|<Ωα, then one may assume n̅(t) is solved through (<ref>) at each moment with updated momentary mechanical frequency Ω(t) and linewidth Γ(t) to yield an effective time dependent coefficients matrix M(t). This offers the solution 

    {A(t)}   =   exp[∫_0^tM(τ)dτ]{A(0)} +∫_0^texp[∫_0^t-τM(ν)dν][β(τ)]{α(τ)}dτ
       +   ∫_0^texp[∫_0^t-τM(ν)dν][Γ(τ)]{A_in(τ)}dτ
    
    [β(t)]^T   =   [ 
    [        1        0    b̅(t)  b̅^∗(t) √(n̅(t)) √(n̅(t));        0        0        0 √(n̅(t))        0 ]], 
    {α(t)}^T   =   {α(t), 
    α^∗(t)
    }.




§ SIDE-BAND INEQUIVALENCE


Let us go back to the set of equations (<ref>) and only retain the first, third, and fourth equations.  This reduction gives a 3× 3 system of equations, identical to (<ref>) with the redefinitions 

    M   =   N+δN
    N   =   [ 
    [      iΔ-κ/2        ig_0        ig_0;          iF -i(Ω-Δ)-γ/2           0;         -iF           0  i(Ω+Δ)-γ/2 ]],
    {δ A(ω)}^T   =   {δâ(ω), δ(âb̂)(ω), δ(âb̂^†)(ω) },
    {A_in(ω)}^T   =   {â_in(ω),b̂_in(ω), 
    b̂_in^†(ω) },
    [√(Γ)]   =   [ 
    [     √(κ)        0        0;   √(κ)b̅   √(Γn̅)        0; √(κ)b̅^∗        0   √(Γn̅) ]].

Here, the perturbation matrix δN is defined through the relation

    δN=[ 
    [    0    0    0; if^+   is    0; if^-    0 is^∗ ]],

in which f^+=g_0(m̅+1) and f^-=g_0m̅. It is quite apparent that the second and third rows of N in (<ref>) are complex conjugates. 

By setting Δ=0 in (<ref>) one would expect identically displaced sidebands at ±Ω. However, this is contingent on the fact that the eigenvalues of N be either complex conjugates as [η]=±Ω corresponding to the frequencies of the sidebands, or [η]=0 corresponding to the resonant pump. However, the presence of perturbation matrix δN breaks this symmetry between the sidebands. This causes a very tiny displacement of sidebands so that Δ_r+Δ_b≠0. First figure of the main article illustrates the side-band asymmetry for various intracavity photon numbers n̅=(g/g_0)^2 and coherent phonon numbers m̅, when g_0/Ω=10^-3. This effect is actually due to the higher-order optomechanical spring effect analyzed in the following. 

It has to be noticed that the horizontal axes are nonlinear functions of the incident light intensity and therefore α. Typically, an inequivalence would be observable in a heterodyne side-band resolved experiment if the effect is large enough to allow clear and measurable motion of side-bands. This condition requires |Δ_r+Δ_b|>Γ=Ω/Q_m, in which Q_m is the mechanical quality factor. If Q_m>10^5, then an intracavity occupation number of n̅>10^4 should be sufficient to detect any such inequivalence.

The side-band inequivalence should not be mistaken with the fundamental energy conservation and time reversal symmetry. Firstly according to these, the spectral density on the negative frequencies of the spectrum should be mirror symmetric with respect to the positive frequencies. Normally, the actual optical frequency ω is much larger than the mechanical frequency Ω, so that the observed red- and blue-detuned sidebands within the range Δ∈(-Ω,+Ω) actually entirely correspond to the positive absolute frequencies. So the speculation that δΔ could be non-zero has nothing to do with the time-reversal symmetry. Secondly, side-band inequivalence is a purely nonlinear effect and is therefore strictly forbidden in any linearized approximation of the optomechanical Hamiltonian. 



§ RESONANCE SHIFT

The contribution of the terms ± iF+if^± to the mechanical frequency Ω in the second and third equations of (<ref>), can be held responsible for the so-called optomechanical spring effect <cit.>. As the result of optomechanical interaction, both of the optical and mechanical resonance frequencies and damping rates undergo shifts. Even at the limit of zero input optical power α=0 and therefore zero cavity photon number n̅=0, it is possible to show that there is a temperature-dependent shift in the mechanical resonance frequency, markedly different from the lattice-expansion dependent effect. This effect is solely due to the optomechanical interaction with virtual cavity photons, which completely vanishes when g_0=0. In close relationship to the shift of resonances, we can also study the optomechanical spring effect with the corrections from higher-order interactions included.

The analysis of spring effect is normally done by consideration of the effective optomechanical force acting upon the damped mechanical oscillator, thus obtaining a shift in squared mechanical frequency δ(Ω^2), whose real and imaginary parts give expressions for δΩ and δΓ. Corrections to these two terms due to higher-order interactions are discussed in the main article. Here, we demonstrate that the analysis using higher-order operator algebra can recover some important lost information regarding the optical and mechanical resonances when the analysis is done on the linearized basis {A}^T={â,â^†,b̂,b̂^†}. 

To proceed, we consider finding eigenvalues of the matrix M as defined in (<ref>). Ignoring all higher-order nonlinear effects beyond the basis {A}^T={â,âb̂,âb̂^†}, we set s=0. This enables us to search for the eigenvalues of the coefficients matrix M as

    eig{M}   =   eig[ 
    [      iΔ-κ/2        ig_0        ig_0;    i(G+f^+) -i(Ω-Δ)-γ/2           0;   -i(G-f^-)           0  i(Ω+Δ)-γ/2 ]]
       =   i{[ Δ+λ_1(Δ,T)+iγ_1(Δ,T); Δ+λ_2(Δ,T)+iγ_2(Δ,T); Δ+λ_3(Δ,T)+iγ_3(Δ,T) ]}
       =   i{[ Δ+η_1(Δ,T); Δ+η_2(Δ,T); Δ+η_3(Δ,T) ]},

in which λ_j=[η_j] and γ_j=[η_j] with j=1,2,3 are real valued functions of Δ and bath temperature T. The temperature T determines m̅ while n̅ is a function of Δ as well as input photon rate α. 

In general, the three eigenvalues η_j=λ_j(Δ,T)+iγ_j(Δ,T), j=1,2,3 are expected to be deviate from the three free-running values ψ_1=i1/2κ, ψ_2=-Ω+i1/2γ, and ψ_3=Ω+i1/2γ, as η_j≈ψ_j-Δ because of non-zero g_0. Solving the three equations therefore gives the values of shifted optical and mechanical frequencies and their damping rates compared to the bare values in absence of optomechanical interactions with g_0=0 

    δΩ   =   -1/2[η_2-η_3]-Ω, 
    δω   =   -1/2[η_2+η_3], 
    δΓ   =   [-2η_1+η_2+η_3]-Γ, 
    δκ   =   2[η_1]-κ.

This  method to calculate the alteration of resonances, does not regard the strength of the optomechanical interaction or any of the damping rates. In contrast, the known methods to analyze this phenomenon normally require g<<κ and Γ+δΓ<<κ <cit.>.

The above values can be calculated numerically for a typical optomechanical cavity, whose parameters are displayed in Table <ref>. The selected values of the four example optomechanical set ups result in very different configurations. System A is very strongly coupled with g/(κ/2)=1.2 and in far Doppler limit (κ/2)/Ω=15. Meanwhile, Systems B and C with (κ/2)/Ω=0.15, and System D with (κ/2)/Ω=2.2× 10^-6 are all in the resolved-side band regime. System B is ultrastrongly coupled with g/(κ/2)=1.07×10^4, while for Systems C and D we have respectively g/(κ/2)=1.2 and g/(κ/2)=3.48.



A very simple way to estimate the shift in eigenvalues is by separating the real and imaginary parts of the optomechanical interaction, as Ω→Ω-[s] and γ→γ+2[s]. These shifts in mechanical frequency and damping rates can be approximated using (<ref>) as

    δΩ+δω   ≈   -g_0[b̅]=-g^2Ω/Ω^2+1/4Γ^2,
    δΓ+δκ   ≈   g_0[b̅]=g^2Γ/2(Ω^2+1/4Γ^2),

where g has been defined under (<ref>). This approximation requires the optomechanical processes {âb̂,âb̂^†} being independent of the other state variables. Since this decoupling is not exact, relations (<ref>) also will remain approximate. However, the accuracy of these are still quite remarkable as is demonstrated here. 

First of all, it is noticed through extensive numerical computations that the shifts in optical and mechanical frequencies take place primarily in the optical part.  That implies the resonance shift is typically much stronger in the optical partition of the system instead of the mechanical partition, leading to marked change in reflection spectra of optomechanical cavities. We calculate and plot each of the four individual components of (<ref>) along with the analytical expressions (<ref>) for the four systems A, B, C and D described in Table <ref>, respectively illustrated in Figs. <ref>, <ref>, <ref> and <ref> for the shifts in mechanical and optical frequencies and dissipation-decay rates. 

Extensive numerical calculations for various configurations establish the fact that it is actually the optical resonance frequency which receives the optomechanical interaction effect. The asymmetry of this shift in cavity optical frequency across the zero-detuning Δ=0, is well exhibited in Fig. <ref> for System B, and in Fig. <ref> for System D, both of which are taken to have relatively large intracavity photon numbers around 10^7 to 10^8. This clear numerical signature underlines the fact that the well-known asymmetry of cavity optical response at high intensities should be actually a result of this higher-order spring effect, rather than thermally induced instabilities. 







Summarizing, any such higher-order resonance shift will cause a change in mechanical frequency δΩ and decay rate δΓ, as well as optical detuning δω=-δΔ and decay rate δκ. While all these four components are non-zero, it is δΔ which is ultimately dominant over the three others in the bistability relation (<ref>). This will make the cavity response to follow the bistability and therefore appear to be asymmetric at high illumination drive intensities.

As a result of higher-order spring effect and δΔ, a shift in intracavity photon number follows δn̅, which immediately shifts the higher-order spring effect and therefore δΔ. The infinite cycle of shifts in intracavity photon number and optical resonance frequency establishes a deterministic chaotic behavior, which is also a well known experimental observation in the community.

Furthermore, for System E, the zero-point optical field can change the mechanical frequency as large as 3.3kHz/K as illustrated in Fig. <ref>. Being in close relationship with the Dynamical Casimir effect <cit.>, this value could be in principle measured if temperature-induced expansion and the resulting change of mechanical frequency is much smaller. The thermal expansion coefficient of Silicon is roughly 2.6×10^-6K^-1, roughly equivalent to 2.6kHz/K. The contribution of zero-point field can be therefore larger or at least within the same order of magnitude. Here, we have assumed that the thermal expansion coefficient of Silicon is independent of temperature and also Ω shifts linearly with temperature. 

The same calculation for System D gives out a value of 0.57kHz/K which is much less than the temperature expansion drift of 13.7kHz/K for the same structure. This phenomenon has been also noticed and referred to as the Nonlinear Transduction <cit.> where the photon-phonon coupling can induce a temperature-dependent change in the resonance frequency of the cavity, even on the order of cavity linewidth. 




§ COHERENT PHONON POPULATION


It is here shown that the method of higher-order operators allows one to find an explicit expression for m̅. In order to do this, we need to write down the 3× 3 reduced set of higher-order optomechanical equations with the fluctuations terms dropped, which reads

    d/dt{[     â;   âb̂; âb̂^† ]}=[[       iΔ-1/2κ          ig_0          ig_0; ig_0(m̅+n̅+1)  -i(Ω-Δ)-1/2γ             0;   ig_0(m̅-n̅)             0   i(Ω+Δ)-1/2γ ]]{[     â;   âb̂; âb̂^† ]}-{[       α;     b̅α; b̅^∗α^∗ ]}.

Here, α should be taken as a complex number from (<ref>), b̅ is substituted from (<ref>) in terms of a̅, where a̅=√(n̅) is taken as a real number and n̅ can already by found from the solution of the third-order algebraic equation (<ref>). 

We are here interested in the steady state solutions, so that the time derivative d/dt=0 can be set to zero. Then we arrive at the system of equations

    [[       iΔ-1/2κ          ig_0          ig_0; ig_0(m̅+n̅+1)  -i(Ω-Δ)-1/2γ             0;   ig_0(m̅-n̅)             0   i(Ω+Δ)-1/2γ ]]{[   a̅;   ab; ab^∗ ]}={[       α;     b̅α; b̅^∗α^∗ ]}.

In the above system of equations, ab corresponds to the time-average of the operators ⟨âb̂|$⟩, whileab^∗corresponds to the time-average of the operators⟨âb̂^†|$⟩. Quite obviously, ab≈a̅b̅ and  ab^∗≈a̅b̅^∗ can approximately hold based on the mean-field approximation. We shall here furthermore observe that this approximation does not any longer hold for the coherent phonons as m̅=⟨b̂^†b̂|≠⟩b̅^∗b̅ for the reasons discussed below.

We now can rearrange (<ref>) in terms of the unknown quantities m̅, ab, and ab^∗ as 

    [[            0         ig_0         ig_0;       ig_0a̅ -i(Ω-Δ)-1/2γ            0;       ig_0a̅            0  i(Ω+Δ)-1/2γ ]]{[   m̅;   ab; ab^∗ ]}={[     α-(iΔ-1/2κ)a̅;  b̅α-ig_0(n̅+1)a̅; b̅^∗α^∗+ig_0 n̅a̅ ]}.

This linear system of equations after appropriate substitutions from (<ref>) and (<ref>) now can be solved to find

    m̅(Δ)=64 g_0^2 Ω ^2 n̅^2(Δ ) (γ ^2+γΓ +4 Δ ^2)-[(Γ ^2+4 Ω ^2)^2 (γ ^2+4 Δ  (Δ +Ω )]/2 (γ ^2+4 Δ ^2) (Γ ^2+4 Ω ^2)^2.

Here, a small imaginary part remains which has to be dropped and results from the inexactness of (<ref>) and (<ref>) coming from linearized optomechanics, and not being in complete consistency with the higher-order formalism. 

In a similar manner, one may find

    ab   =   i √(n̅)[g_0^2 (8 n̅+4)+(2 Δ +i κ ) (2 (Δ +Ω )+i γ )]-2 i α[8 Γ  g_0^2 n̅/Γ ^2+4 Ω ^2+γ -2 i (Δ +Ω )]/4 g_0 (γ -2 i Δ ), 
    ab^∗   =   2 α[8 i Γ  g_0^2 n̅/Γ ^2+4 Ω ^2-i γ -2 Δ +2 Ω]-√(n̅)[4 i g_0^2 (2 n̅+1)+(κ -2 i Δ ) (i γ +2 Δ -2 Ω )]/4 g_0 (γ -2 i Δ ).



The expression (<ref>) for m̅ is accurate within half a quanta ±1/2, so that in order to satisfy the zero limits at infinite detuning 

    lim_Δ→∞m̅(Δ)=0,

a half-quanta must be added to (<ref>). Then it will read

    m̅(Δ)   =   32 g_0^2 Ω ^2 (γ ^2+γΓ +4 Δ ^2)/(γ ^2+4 Δ ^2) (Γ ^2+4 Ω ^2)^2n̅^2(Δ ) -2 ΔΩ/γ ^2+4 Δ ^2 ±1/2
       ≈   32 g_0^2 Ω ^2 (γ ^2+γΓ +4 Δ ^2)/(γ ^2+4 Δ ^2) (Γ ^2+4 Ω ^2)^2n̅^2(Δ ) 
       =   g_0^2ζ(Δ)n̅^2(Δ).

The approximation holds well if n̅ is well above unity. Hence, we can infer from (<ref>) that m̅∝n̅^2. In the lossless limit, where γ≈ 0 and Γ≈ 0, one may even further simplify (<ref>) to obtain the simple expression m̅≈ 2g_0^2n̅^2/Ω^2 which is typically accurate within 10% of the actual value or better. It is not difficult to check the resonant coherent phonon number m̅(0). In the practical limit of κ>>Γ, it is easy to verify that (<ref>) actually simplifies to m̅(0)≈ 32 [g_0 Q_mn̅(0)/Γ]^2 with Q_m=Ω/Γ being the mechanical quality factor. In the next section, we point out a straightforward method to measure this quantity through experiment on the well-known optical spring effect.

In practice, the expression (<ref>) is sensitive to the choice of optomechanical parameters and in particular g_0. A slight variation in the basic optomechanical parameters {g_0,ω,κ,Ω,Γ} with γ=κ+Γ as small as few percent can make a pronounced effect in expected behavior in m̅.

For the side-band resolved systems in the lossless limit, it is within 10% of the relationship m̅∝ 2|b̅|^2, meanwhile for Doppler cavities, the agreement is roughly within 3% or better. This result perfectly agrees to the large-amplitude oscillation limit of b̅(t)≈δb̂(t)+[b̅+b̅exp(-iΩ t)] where δb̂ represents the random fluctuations in the mechanical field with ⟨δb̂|=⟩0. This also tells that the coherent oscillations of the mechanical field are not differential in amplitude, and can vary in the range (0,2|b̅|). So, the amplitude of coherent mechanical oscillations is just as big as their average. This large-amplitude coherent mechanical wave is driven and waked by the optical coherent field inside the cavity, through optomechanical interactions. The mean field approximations ab≈a̅b̅ and ab^∗≈a̅b̅^∗ seem however to always hold better than 0.1% for Doppler cavities. This accuracy breaks down for side-band resolved cavities.

If there are more than one mechanical fields available j=1,2,⋯, the coherent phonon population of each mode m̅_j shall be determined with the corresponding sets of optomechanical parameters {g_0,j,ω_j,κ_j,Ω_j,Γ_j}, with the expected approximate result m̅_j≈ 2g_0,j^2n̅_j^2/Ω_j^2 as long as the mechanical modes are almost uncorrelated. The case of coherent phonon numbers of two or more correlated mechanical modes needs a separate study.

It is here again stressed out that the oscillations of the mechanical field can be decomposed into the incoherent and coherent parts. The incoherent part results from random thermal fluctuations with the thermal occupancy m, as well as half a quanta contributing from the quantum noise of the coherent part, while the coherent oscillations correspond to the coherent phonon number m̅. The same also should be true for the optical field, however, the random fluctuations of a coherent light is only half a quanta, and the thermal optical occupancy n of optomechanical cavity is normally negligible under practical considerations and working temperatures. 



§ HIGHER-ORDER SPRING EFFECT


It is possible to calculate the optomechanical spring effect due to the standard linearized and higher-order interactions. In order to do this, we start from the matrix [M] given in (<ref>), and after dropping the noise and drive input terms we notice the expansion of first Langevin equation for the operator â. That reads

    d/dtâ=(iΔ-1/2κ)â+ig_0â(b̂+b̂^†).

From the second and third equations we get

    âd/dtb̂+[(iΔ-1/2κ)â+ig_0â(b̂+b̂^†)]b̂   =   i(f^+ +F)â-[i(Ω-Δ-g_0b̂)+1/2γ]âb̂, 
    âd/dtb̂^†+[(iΔ-1/2κ)â+ig_0â(b̂+b̂^†)]b̂^†   =   i(f^- -F)â+[i(Ω+Δ+g_0b̂)-1/2γ]âb̂^†.

This is equivalent to 

    d/dtb̂   =   i(f^+ +F+b̂^†b̂)-(iΩ+1/2Γ)b̂, 
    d/dtb̂^†   =   i(f^- -F+b̂b̂^†)+(iΩ-1/2Γ)b̂^†.

These two equations can be now combined after dropping the nonlinear terms by addition and subtraction, and then taking the Fourier transform to yield the system

    [
    [ -iw+1/2Γ       iΩ;       iΩ -iw+1/2Γ ]]{[ b̂+b̂^†; b̂-b̂^† ]}=2ig_0{[ m̅+1/2; n̅+1/2 ]}.

This can be solved now to yield the expression for δx̂=x_zp(b̂+b̂^†) as 

    δx̂(w)=2ix_zpg_0(-iw+1/2Γ)(m̅+1/2)-iΩ(n̅+1/2)/(-iw+1/2Γ)^2-(-iΩ)^2.

A rearrangement of this expression yields

    δx̂(w)=2x_zpg_0Ω/-(w+i1/2Γ)^2+Ω^2{n̅+[(w/Ω+iΓ/2Ω)(m̅+1/2)+1/2]}.

It is straightforward now to see that the term within brackets contributes to the necessary corrections to the spring effect. This will change the mechanical response function Σ(w) <cit.> as

    Σ(w,Δ)=2Ω g_0^2 [1/(Δ+w)+i/2κ+1/(Δ-w)-i/2κ][n̅+μ(w)],

where a term with the dimension of mass in the numerator, which in the following calculation ultimately cancels out, and is equal to the effective motion mass m_eff, is not shown for simplicity. We have also

    μ(w)=1/Ω(w+i/2Γ)(m̅+1/2)+1/2,

represents corrections to yield the effective cavity photon number n̅_eff=n̅+μ(w) because of higher-order interactions. This corrections is easy to see that are important if the pump level is not too high. Typically, for n̅<10^2 higher-order spring effects are quite significant, and when n̅>10^3 the higher-order effects are suppressed by the standard spring effect.

The spring effect modifies the measured mechanical frequency Ω and linewidth Γ as

    δΩ(w,Δ)   =   1/2w[Σ(w,Δ)], 
    δΓ(w,Δ)   =   -1/w[Σ(w,Δ)].

Put together combined, we get

    δΩ(w,Δ)   =   g_0^2n̅Ω/w[Δ+w/(Δ+w)^2+1/4κ^2+Δ-w/(Δ-w)^2+1/4κ^2]
       +   g_0^2[μ(w)]Ω/w[Δ+w/(Δ+w)^2+1/4κ^2+Δ-w/(Δ-w)^2+1/4κ^2] 
       +   g_0^2[μ(w)]Ω/w[κ/(Δ+w)^2+1/4κ^2-κ/(Δ-w)^2+1/4κ^2],
    δΓ(w,Δ)   =   g_0^2n̅Ω/w[κ/(Δ+w)^2+1/4κ^2-κ/(Δ-w)^2+1/4κ^2]
       +   g_0^2[μ(w)]Ω/w[κ/(Δ+w)^2+1/4κ^2-κ/(Δ-w)^2+1/4κ^2]
       -   g_0^2[μ(w)]Ω/w[Δ+w/(Δ+w)^2+1/4κ^2+Δ-w/(Δ-w)^2+1/4κ^2].

Here, the second and third terms on the rights hand sides of both equations are corrections to the spring effect due to the higher-order interactions, resulting from the temperature-dependent expressions 

    [μ(w)]   =   w/Ω(m̅+1/2)+1/2, 
    [μ(w)]   =   Γ/2Ω(m̅+1/2).

The temperature-dependence of (<ref>) causes dependence of the spring effect on temperature as well.

The uncorrected standard expressions read <cit.>

    δΩ(w,Δ)   =   g_0^2n̅Ω/w[Δ+w/(Δ+w)^2+1/4κ^2+Δ-w/(Δ-w)^2+1/4κ^2],
    δΓ(w,Δ)   =   g_0^2n̅Ω/w[κ/(Δ+w)^2+1/4κ^2-κ/(Δ-w)^2+1/4κ^2],

from which we may observe

    δΩ(w,Δ)   =   -δΩ(w,-Δ),
    δΩ(-Δ,Δ)   =   g_0^2n̅Ω/2Δ/Δ^2+1/16κ^2,

which do hold for the standard spring effect at sufficiently high optical powers. 



 §.§ Examples


As an example, we first assume a side-band resolved cavity with optical resonance wavelength of λ=1μm and quality factor Q=10^6, a mechanical resonance frequency of Ω=2π× 1GHz and quality factor of Q_m=10^5, with the optomechanical interaction rate g_0=2π× 160kHz. Pumping resonantly at the rates α=10^10s^-1 and α=10^11s^-1 correspond to n̅=1.1× 10^2 and n̅=1.1× 10^4.

In Figs. <ref> and <ref>, optical spring effect has been illustrated at the two above different pumping rates, respectively, for a side-band resolved cavity. It has been supposed that a two-beam pump-probe experiment is undertaken, where the pump frequency is having the detuning Δ and the relative probe frequency is w. 

It can be seen that when the cavity photon number n̅ is increased while coherent phonon number m̅ is kept constant, the higher-order corrections to the spring effect are just negligible and can be ignored, as shown in Fig. <ref>. However, at lower pumping rates where n̅ is no longer much larger than m̅, the higher-order corrections become important as shown in Fig. <ref>. Ultimately, for a few or very low cavity photon number n̅<10, the higher-order corrections could be orders of magnitude stronger than the standard effects. 

In Fig. <ref> the same cavity parameters are employed with a lowered optical quality factor of Q=10^4. This choice puts the cavity in the Doppler regime. It is easy to see again that in the same manner of the previous examples, sufficiently high pump levels entirely masks the higher-order effects. Taking only the case of α=3× 10^11s^-1 here leads to n̅=1.01× 10^1, while at smaller pump levels only the higher-order optomechanical effects survive. At this cavity photon number around n̅≈ 10, differences between the standard and higher-order corrected responses are quite prominent and visible.










§ MINIMAL BASIS


Complete solution of optomechanical interaction ℍ_OM can be attained analytically using the minimal basis {A}^ T={n̂^2,n̂b̂,n̂b̂^†}={N̂,B̂,B̂^†}. Construction of Langevin equations leads to the system

    d/dt{[   N̂;   B̂; B̂^† ]}=
    [
    [     -2κ       0       0;    ig_0 -iΩ-γ/2       0;    ig_0       0  iΩ-γ/2 ]]
    {[   N̂;   B̂; B̂^† ]}-{[  √(4κ)N̂_ in;   √(γ)B̂_ in; √(γ)B̂_ in^† ]}.

Here, the multiplicative noise terms are defined as 

    √(4κ)N̂_in   =   2√(κ)(n̂â^†â_in+â_in^†ân̂),
    √(γ)B̂_in   =   √(2κ)b̂n̂_in+√(Γ)n̂b̂_in,
 
where n̂_in is already defined under (<ref>), and the spectral density of which can be estimated using the method described elsewhere <cit.>. A very effective method to deal with multiplicative noise is to be discussed in <ref>. This can be immediately noticed to be reducible as 

    d/dt{[ N̂; B̂ ]}=
    [
    [     -2κ       0;    ig_0 -iΩ-γ/2 ]]
    {[ N̂; B̂ ]}-{[ √(4κ)N̂_ in;  √(γ)B̂_ in ]}.

These will make the evaluation of spectral densities S_NN(ω) and S_BB(ω) possible. Interestingly, (<ref>) is actually decoupled, since the equation for N̂ is already independent of B̂, which admits the solution

    N̂(t)=N̂(0)e^-2κ t-2√(κ)e^-2κ t∫_0^tN̂_in(τ)e^2κτdτ.

We can be now plug (<ref>) in the second equation of (<ref>) to solve exactly for B̂. We define ϑ=iΩ+γ/2 and may write down

    B̂(t)=B̂(0)e^-ϑ t-e^-ϑ t∫_0^te^ϑτ[ig_0N̂(τ)+√(γ)B̂_in(τ)]dτ.


The treatment of multiplicative noise terms (<ref>) can be quite difficult in the most general form, especially that they demand prior knowledge of photonic and phononic ladder operators. However, assuming that the extra ladder operators can be replaced by their mean values, we can do the zeroth order approximations 

    √(4κ)N̂_in   ≈   √(κn̅)n̅(â_in+â_in^†)→ 2√(κn̅)n̅ǎ_in,
    √(γ)B̂_in   ≈   √(κn̅)b̅(â_in+â_in^†)+√(Γ)n̅b̂_in→ 2√(κn̅)b̅ǎ_in+√(Γ)n̅b̂_in.
 
Here, the real-valued Weiner process ǎ_in(t) with the symmetrized classical spectral density ǎ_in(ω) is obtained as

    ǎ_in(t)   =   â_in(t)+â_in^†(t)/2,
    ǎ_in(ω)   =   [â_in(ω)].


While this type of approximations in multiplicative noise could be useful for many cases, there are some phenomena which cannot be reproduced without correct treatment of multiplicative noise. This shall be discussed in details in <ref>. Nevertheless, it is also instructive the take the expectation values of (<ref>) to obtain the classical system

    d/dt{[ N(t); B(t) ]}=
    [
    [     -2κ       0;    ig_0 -iΩ-γ/2 ]]
    {[ N(t); B(t);      ]}+2√(n̅){[ n̅; b̅ ]}[α].

Together with (<ref>,<ref>), and setting the time-derivative on the left of the above to zero, makes the evaluation of steady-state values N(∞)=n̅^2 and B(∞)=nb≈n̅b̅ readily possible. Doing this and solving for n̅ and b̅ precisely gives back (<ref>). This not only is in agreement with the equilibrium equation (<ref>), but also confirms the general finding that the equilibrium intracavity photon population n̅(Δ) is independent of the coherent phonon population m̅(Δ). However, the opposite is not correct, and as it was shown in the previous sections, m̅(Δ) can actually be either determined from n̅(Δ) and fitting to the experimental data, or directly estimated using the expression (<ref>) in <ref>. 

Existence of such an exact transformation which puts the optomechanical interaction into exactly linear form should be connected to the polaron transformation <cit.> which leaves behind a Kerr nonlinear term as n̂^2 in the transformed optomechanical Hamiltonian. It furthermore highlights the fact that usage of higher-order operators ultimately can reach a fully linear system at which convergence of this method to the exact solution is evident.



§ HIGHER-ORDER SIDEBANDS


When the optical frequency is much larger than the mechanical frequency, apart from the mechanical sidebands which are roughly placed at Δ^(1)≈±Ω, there exist higher-order sidebands such as Δ^(2)≈±2Ω and so on. The occurrence of these higher-order sidebands, which are observable for sideband-resolved experiments, is obviously stringent on the existence of two- and multi-phonon processes. Normally, one would expect that these could be studied by constructing the Langevin equations for the operators b̂^2, b̂^† 2 and so on. But this guess turns out to be incorrect, since the corresponding Langevin equations would be totally independent of the one for â, implying that the second- and higher-order sidebands could not be reconstructed via the fully linearized Langevin equations. This has already been shown to be a nonlinear process which does not naturally appear in the Hamiltonian of the fully linearized optomechanics <cit.>. But this difficulty can be appropriately addressed by the method of Higher-order Operators, too.

In order to investigate this phenomenon, let us restrict the case only to the second-order sidebands roughly located at Δ^(2)≈±2Ω. In order to study these, it is sufficient to extend the basis {A}^T={â,âb̂,âb̂^†} to 

    {A}^ T={â,âb̂,âb̂^†,âb̂^2,âb̂^† 2},

where the third-rank higher-order operators âb̂^2 and âb̂^† 2 take care of the one-photon two-phonon processes, ultimately leading to formation of second-order sidebands at Δ^(2)≈±2Ω. The Langevin equations for this basis within the zeroth-order approximation of multiplicative noise reads

    [ 
    [         iΔ-κ/2           ig_0           ig_0              0              0;  ig_0(m̅+n̅+1)    -i(Ω-Δ)-γ/2              0           ig_0              0;    ig_0(m̅-n̅)              0     i(Ω+Δ)-γ/2              0           ig_0;              0 ig_0(m̅+2n̅+2)              0   -i(2Ω-Δ)-θ/2              0;              0              0 ig_0(m̅-2n̅-1)              0    i(2Ω+Δ)-θ/2 ]]
    ×{[       â;     âb̂;   âb̂^†;   âb̂^2; âb̂^† 2 ]}-[ 
    [      √(κ)         0         0; √(1/2κm̅)    √(Γn̅)         0; √(1/2κm̅)         0    √(Γn̅); 1/2√(κ)m̅  √(Γn̅m̅)         0; 1/2√(κ)m̅         0  √(Γn̅m̅) ]]
    {[    â_in;   b̂_in; b̂_in^† ]}+
    [ 
    [      1      0;     b̅      0;   b̅^∗      0;   b̅^2      0; b̅^∗ 2      0;        ]]
    {[   α; α^∗ ]} 
    =d/dt{[       â;     âb̂;   âb̂^†;   âb̂^2; âb̂^† 2 ]}.

Here, θ=κ+2Γ is the decay rate associated with the third-rank one-photon two-phonon processes âb̂^2 and âb̂^† 2. The approximation 2|b̅|^2≈m̅ is used following the discussions in <ref>. We do observe that this treatment of multiplicative noise causes non-negligible error in some cases, and is due to be discussed later in <ref>.



Vertical and horizontal partitions separate the single-phonon {âb̂,âb̂^†} and two-phonon processes {âb̂^2,âb̂^†^2}. So by retaining only the first 3× 3 blocks and first 3 rows what remains is nothing but the equations of first-order optomechanics in terms of the second-rank single-phonon operator basis {â,âb̂,âb̂^†}. 

The first-order δΔ^(1) and second-order δΔ^(2) sideband inequivalences take on similar expansions as

    δΔ^(1)/Ω≈-δΔ^(2)/2Ω≈(g_0/Ω)^2(n̅+1/2)-2(g_0/Ω)^4(n̅+1/2)(m̅+1/2) ≈𝒢^2-4𝒢^4𝒢_0^2.

Here, 𝒢_0=g_0/Ω and 𝒢=g/Ω are normalized interaction rates with respect to the mechanical frequency, where g=g_0√(n̅) is the enhanced optomechanical interaction rate. Furthermore, m̅ is approximated from (<ref>) in the above.


Results of noise spectrum calculations using the fully linearized and higher-order formulations of optomechanics is shown in Fig. <ref>. The single-photon optomechanical interaction rate g_0=1.68× 10^-3Ω and the enhanced optomechanical interaction on resonance satisfies g=0.31Ω, corresponding to the strong coupling regime and optical power of P_op=47.2μW at T=3K. Cavity is side-band resolved with the parameters given elsewhere <cit.>. The input power is high-enough to cause the cavity to exhibit asymmetric reflectivity, a very clear hallmark of bistability seen easily in experiments. 

This has been calculated and illustrated in Fig. <ref> for various linear and higher-order formulations resulting from simulating a scanning pump experiment. It can be seen that the fully-linearized optomechanics cannot reasonably reproduce the highly asymmetric and non-Lorentzian lineshape of cavity under strong pump. As a simple measure of reflectivity, one may use the Langevin equation for photons â with the semi-classical substitutions  â→a̅ and b̂→b̅, where b̅ is correspondingly given from (<ref>), and n̅ can be nonlinearly solved from (<ref>). For a side-coupled cavity where the reflectivity is not identity, and external coupling rate κ_ex is known, we have η=κ_ex/κ, leading to the approximation

    ℛ(ω,Δ)=1-iκ_ex/ω+Δ+2g_0^2n̅Ω/Ω^2+1/4Γ^2+i1/2κ,
 
where |ℛ|^2 is plotted as the dot-dashed curve in Fig. <ref>. More accurate solutions can be found by taking the scattering matrix element ℛ=Y_11, also shown in Fig. <ref>.

The existence of tiny second-order sidebands around ± 2Ω is illustrated by calculation of the reflectivity near the corresponding resonances. This has been shown in Fig. <ref> for both of the second-order red and blue sidebands at P_op=3mW. The depths of these resonances are rather small, being only around -0.13mdB.








§ MULTIPLICATIVE NOISE


Rewriting (<ref>) without the zeroth-order approximation for multplicative noise, gives the exact higher-order set of Langevin equations 

    d/dt{[       â;     âb̂;   âb̂^†;   âb̂^2; âb̂^† 2 ]}   =   [ 
    [         iΔ-κ/2           ig_0           ig_0              0              0;  ig_0(m̅+n̅+1)    -i(Ω-Δ)-γ/2              0           ig_0              0;    ig_0(m̅-n̅)              0     i(Ω+Δ)-γ/2              0           ig_0;              0 ig_0(m̅+2n̅+2)              0   -i(2Ω-Δ)-θ/2              0;              0              0 ig_0(m̅-2n̅-1)              0    i(2Ω+Δ)-θ/2 ]]
    {[       â;     âb̂;   âb̂^†;   âb̂^2; âb̂^† 2 ]}
       -   [ 
    [       √(κ)          0          0;     √(κ)b̂      √(Γ)â          0;   √(κ)b̂^†          0      √(Γ)â;   √(κ)b̂^2    √(Γ)âb̂          0; √(κ)b̂^† 2          0  √(Γ)âb̂^† ]]
    {[    â_in;   b̂_in; b̂_in^† ]}+
    [ 
    [      1      0;     b̅      0;   b̅^∗      0;   b̅^2      0; b̅^∗ 2      0;        ]]
    {[   α; α^∗ ]}.


To illustrate how the multiplicative noise terms on the second line are to be treated, let us assume that a simple equation is given as 

    d/dt𝒜(t)=(iΔ-1/2κ)𝒜(t)-√(κ)x̂(t)â_in(t),

where x̂ is some dimensionless and time-dependent operator and â_in corresponds to a white noise random process. The spectral density of the zero-average operator 𝒜 by definition is

    S_𝒜𝒜(w)=∫_-∞^∞dτ e^iwτ⟨𝒜^†(t)𝒜(t+τ)|.⟩

Therefore, the symmetrized spectral density via symmetrization operator 𝒮 which is the actual quantity measured in experiments is

    𝒮S_𝒜𝒜(w)   =   ∫_-∞^∞dτ e^iwτ⟨𝒮{𝒜^†(t)𝒜(t+τ)}|
    ⟩   =   ∫_-∞^∞dτ e^iwτ⟨𝒜^†(t)𝒜(t+τ)|_⟩S.

The equation (<ref>) admits a formal solution

    𝒜(t)=-√(κ)𝕃(t)x̂(t)â_in(t),

where 𝕃 is given as

    𝕃(t)=(d/dt-iΔ+1/2κ)^-1,

and is an operator which can be understood as an inverse Fourier transform such as

    𝕃(t)=ℱ^-1(1/iw-iΔ+1/2κ)(t)=ℱ^-1{L(w)}(t).

Here, we are not interested in an explicit form of 𝕃 although it is easy to be evaluated or looked up from table of Fourier transforms.

The formal solution (<ref>) gives rise to the symmetrized spectral density

    𝒮S_𝒜𝒜(w)=κ∫_-∞^∞dτ e^iwτ⟨â^†_in(t)x̂^†(t)𝕃^†(t)𝕃(t+τ)x̂(t+τ)â_in(t+τ)|_⟩S.

Now, we can employ the Isserlis-Wick theorem to decompose the expectation value as <cit.>

    ⟨â^†_in(t)ŷ^†(t)ŷ(t+τ)â_in(t+τ)|_⟩S   =   ⟨â^†_in(t)â_in(t+τ)|_⟩S⟨ŷ^†(t)ŷ(t+τ)|_⟩S
       +   ⟨â^†_in(t)ŷ^†(t)|_⟩S⟨ŷ(t+τ)â_in(t+τ)|_⟩S
       +   ⟨â^†_in(t)ŷ(t+τ)|_⟩S⟨ŷ^†(t)â_in(t+τ)|_⟩S.

where ŷ(t)=𝕃(t)x̂(t) is adopted for shorthand notation. Since â_in is a white noise Wiener random process, we may expect that to a very good approximation the second and third terms vanish and thus

    𝒮S_𝒜𝒜(w)=κ∫_-∞^∞dτ e^iwτ⟨â^†_in(t)â_in(t+τ)|_⟩S⟨ŷ^†(t)ŷ(t+τ)|_⟩S.

The random nature of a Weiner process requires that <cit.> 

    𝒮S_𝒜𝒜(w)   =   |κ∫_-∞^∞dτ e^iwτ⟨ŷ^†(t)ŷ(t+τ)|_⟩S|^2𝒮S_AA(w)
       =   |L(w)∗x̂(w)|^2𝒮S_AA(w).

Here, ∗ represents convolution in frequency and L(w) is defined in (<ref>) and actually represents the equivalent to the scattering matrix element. The expression (<ref>) presents a mathematically exact solution to the spectral density problem of multiplicative noise (<ref>).

In the context of higher-order quantum optomechanics and referring to (<ref>) the operator x̂ may represent either of the operators within the set {â,b̂,b̂^†,âb̂,âb̂^†,b̂^2,b̂^† 2}. However, not only these are not yet known, but also, they are influenced by random processes from the correspondingly lower-order interactions with the optical field and mechanical bath. The only approximation needed here is to replace these with corresponding non-operator functions which can be already obtained from the solution to lower-order equations. Doing this results in a set of equations for {a̅,b̅,b̅^∗,ab,ab^∗,b̅^2,b̅^∗ 2}, where solutions for {b̅,b̅^∗,b̅^2,b̅^∗ 2} can be obtained by having b̅. This is here calculated from the 3× 3 linearized optomechanics, giving rise to the expression 

    b̅(w)=α/i (w+Ω)+1/2Γ.

In a similar manner to <ref>, the next required expressions can be explicitly obtained  by Mathematica as

    a̅(w)   =   -α ^2 g_0 (w -Δ -Ω-i 1/2γ)/(w+Ω-1/2i Γ) {2 g_0^2[ (Δ-w+1/2 i γ)  (m̅+1/2)+Ω(n̅+1/2)]+(w-Δ -1/2i κ) [(w-Δ -1/2i γ)^2-Ω ^2]}, 
    ab(w)   =   α ^2 [g_0^2 (m̅-n̅)-(w-Δ -1/2i κ) (w -Δ -Ω -1/2i γ)]/(w+Ω-1/2i Γ) {2 g_0^2 [(Δ-w+1/2i γ)(m̅+1/2)+Ω(n̅+1/2)]+(w-Δ -1/2i κ) [(w-Δ -1/2i γ)^2-Ω ^2]}, 
    ab^∗(w)   =   -| α| ^2 [g_0^2 (m̅+n̅+1)-(w-Δ -1/2i κ) (w-Δ+Ω-1/2i γ)]/(w+Ω+1/2i Γ) {2 g_0^2 [(Δ-w+1/2i γ)(m̅+1/2)+Ω(n̅+1/2)]+(w-Δ -1/2i κ) [(w-Δ -1/2i γ)^2-Ω ^2]}.

In the above equations, it has to be noticed that α is a complex quantity which satisfies |α|=√(κ_exP_op/ħω), and also by (<ref>) we have

    α=√(n̅)[-κ/2+i(Δ+2g_0^2Ω/Ω^2+1/4Γ^2n̅)].


Now, we can rewrite the Langevin equations (<ref>) as

    d/dt{δ A}   =   [M]{δ A}-[Ĝ]{ A_in},
    
    [M]   =   [ 
    [         iΔ-κ/2           ig_0           ig_0              0              0;  ig_0(m̅+n̅+1)    -i(Ω-Δ)-γ/2              0           ig_0              0;    ig_0(m̅-n̅)              0     i(Ω+Δ)-γ/2              0           ig_0;              0 ig_0(m̅+2n̅+2)              0   -i(2Ω-Δ)-θ/2              0;              0              0 ig_0(m̅-2n̅-1)              0    i(2Ω+Δ)-θ/2 ]], 
    {δ A}^T   =   {δâ,δ(âb̂),δ(âb̂^†),δ(âb̂^2),δ(âb̂^† 2)},
    { A_in}^T   =   {â_in,b̂_in,b̂_in^†}, 
    
    [Ĝ]   =   [ 
    [       √(κ)          0          0;     √(κ)b̂      √(Γ)â          0;   √(κ)b̂^†          0      √(Γ)â;   √(κ)b̂^2    √(Γ)âb̂          0; √(κ)b̂^† 2          0  √(Γ)âb̂^† ]].

After defining the decay matrix

    [√(Γ)]=[[ √(κ)    0    0    0    0;    0 √(γ)    0    0    0;    0    0 √(γ)    0    0;    0    0    0 √(θ)    0;    0    0    0    0 √(θ) ]],

taking the Fourier transform, and using the input-output relation 

    {A_out(ω)}={A_in(ω)}+[√(Γ)]^T{δ A(ω)},

we arrive at the definition of the scattering matrix

    [Y(ω)]=I-[√(Γ)]^T([M]-iω[I])^-1[√(Γ)],

by which and (<ref>) we can evaluate the desired symmetrized spectral density of output optical field as

    S(ω)   =   |Y_11(ω)|^2 S_AA(ω)
       +   1/γ^2|[Y_12(ω)+Y_13(ω)]∗a̅(ω)|^2 S_BB(ω) 
       +   1/θ^2|[Y_14(ω)∗ab(ω)+Y_15(ω)∗ab^∗(ω)]|^2 S_BB(ω),

where spectral power densities S_AA and S_BB are already introduced in (<ref>) and convolutions ∗ take place over the entire frequency axis. 



§ ELEMENTS OF HIGHER-ORDER SCATTERING MATRICES

This section reports the explicit elements of the first row of scattering matrix [Y] in (<ref>), as needed for calculation of the spectral density according to (<ref>). These might be useful only when the method of residues are to be used for exact evaluation of complex convolution integrals, otherwise full numerical simulation of (<ref>) is much preferable.  



 §.§ Second-order 3× 3 Formalism


The elements of the scattering matrix are explicitly found using the supplementary Mathematica packages, and after some simplification they take the form

    Y_11(ω)   =   1-2iκ_ex[(ω-Δ -1/2i γ)^2-Ω ^2]/2 g_0^2[ (Δ-ω+1/2 i γ)  (m̅+1/2)+Ω(n̅+1/2)]+(ω-Δ -1/2i κ) [(ω-Δ -1/2i γ)^2-Ω ^2], 
    
    Y_12(ω)   =   -ig_0√(γκ_ex)(ω-Δ-Ω -1/2i γ)/2 g_0^2[ (Δ-ω+1/2 i γ)  (m̅+1/2)+Ω(n̅+1/2)]+(ω-Δ -1/2i κ) [(ω-Δ -1/2i γ)^2-Ω ^2], 
    
    Y_13(ω)   =   -ig_0√(γκ_ex)(ω-Δ+Ω -1/2i γ)/2 g_0^2[ (Δ-ω+1/2 i γ)  (m̅+1/2)+Ω(n̅+1/2)]+(ω-Δ -1/2i κ) [(ω-Δ -1/2i γ)^2-Ω ^2].

These expressions are useful in speed up of the code, as well as wherever the method of residues is to be used.



 §.§ Third-order 5× 5 Formalism


The convergence of 3× 3 is sufficiently good for most practical purposes, and also the explicit expressions for 5× 5 matrices decompose into products of fourth-order polynomials in terms of ω in their denominators, which severely limits the usefulness of applicability of the method of residues. For this reason, their explicit expressions are not included here. The interested reader may find them in the supplementary Mathematica packages instead.



§ MATHEMATICA PACKAGES

These are brief descriptions of supplied Mathematica packages along with this article, written by the author:

	
  * SuppleMath1.nb: Derivation of coherent phonon population m̅(Δ) in <ref>.
	
  * SuppleMath2.nb: Derivation of expressions corresponding to <ref> and <ref>.
	
  * SuppleMath3.nb: Code for generation of noise spectra in the main article including animated graphs.
	
  * SuppleMath4.nb: Code for generation of stability diagrams in the main article.
	
  * SuppleMath5.nb: Code for generation of higher-order optical spring effect in Fig. <ref>.
	
  * SuppleMath6.nb: Code for generation of higher-order optical spring effect in  Fig. <ref>.
	
  * SuppleMath7.nb: Code for generation of higher-order optical spring effect in  Fig. <ref>.
	
  * SuppleMath8.nb: Numerical calculation of side-band inequivalence.
	
  * SuppleMath9.nb: Code for generation of normalized side-band inequivalence in of the main article.
	
  * SuppleMath10.nb: Code for generation of reflection spectra in Fig. <ref>.
	
  * SuppleMath11.nb: Code for generation of second-order side-bands in Fig. <ref>.


99
	
	SBowen
	Bowen, W. P. & Milburn, G. J. Quantum Optomechanics (CRC Press: Boca Raton, 2016).
	SAspel1
	Aspelmeyer, M, Kippenberg, T. J. & Marquardt, F. Cavity optomechanics. Rev. Mod. Phys. 86, 1391 (2014).
	SPaper2
	Khorasani, S. Higher-order interactions in quantum optomechanics: Analytical solution of nonlinearity. Photonics 4, 48 (2017).
	SKip1
	Kippenberg, T. J. & Vahala, K. J. Cavity optomechanics. Science 321,  1172 (2008).
	SQuad9
	Doolin, C., Hauer, B. D., Kim, P. H., MacDonald, A. J. R., Ramp, H. & Davis, J. P. Nonlinear optomechanics in the stationary regime. Phys. Rev. A 89, 053838 (2014).
	SSpring1
	Teufel, J. D., Harlow, J. W., Regal, C. A. &  Lehnert, K. W. Dynamical backaction of microwave fields on a nanomechanical oscillator. Phys. Rev. Lett. 101, 197203 (2008).
	SSpring2
	Eichenfield, M., Camacho, R., Chan, J., Vahala, K. J. &  Painter, O. A picogram- and nanometre-scale photonic-crystal optomechanical cavity. Nature 459, 550 (2009).
	SSpring3
	Safavi-Naeini, A. H., Gröblacher, S., Hill, J. T., Chan, J., Aspelmeyer, M. &  Painter, O. Squeezed light from a silicon micromechanical resonator. Nature 500,  185 (2013).
	SSpring4
	Deotare, P. B., Bulu, I., Frank, I. W., Quan, Q., Zhang, Y., Ilic, R. &  Loncar, M. All optical reconfiuration of optomechanical fiters. Nat. Commun. 3, 846 (2011).
	SSpring5
	Sarabalis, C. J., Dahmani, Y. D., Patel, R. N., Hill, J. T. &  Safavi-Naeini, A. H.. Release-free silicon-on-insulator cavity optomechanics. Optica 4, 1147 (2017).
	SAspel3
	Hong, S., Riedinger, R., Marinkovic, I., Wallucks, A., Hofer, S. G., Norte, R. A., Aspelmeyer, M. & Gröblacher, S. Hanbury Brown and Twiss interferometry of single phonons from an optomechanical resonator. Science 358, 203 (2017).
	SAspel4
	Riedinger, R., Wallucks, A., Marinkovic, I., Löschnauer, C., Aspelmeyer, M., Hong, S. & Gröblacher, S. Remote quantum entanglement between two micromechanical oscillators. Nature 556, 473 (2018).
	
	
	SMacri
	Macrì, V., Ridolfo, A., Di Stefano, O., Kockum, A. F., Nori, F. & Savasta, S. Non-perturbative dynamical Casimir effect in optomechanical systems: Vacuum Casimir-Rabi splittings. Phys. Rev. X 8, 011031 (2018).
	STransduction
	Leijssen, R. & Verhagen, E. Strong optomechanical interactions in a sliced photonic crystal nanobeam. Sci. Rep. 5, 15974 (2015).
	SKip2
	Schliesser, A. & Kippenberg, T. J. Cavity optomechanics with whispering-gallery-mode optical micro-resonators. Adv. At., Mol. Opt. Phys. 58, 207 (2010). 
	SLoo
	Loo, S. G. Spectral density of random signals corrupted by multiplicative noise. Electron. Lett. 3, 238 (1967).
	SGirvin
	Børkje, K., Nunnenkamp, A., Teufel, J. D. & Girvin, S. M. Signatures of nonlinear cavity optomechanics in the weak coupling regime. Phys. Rev. Lett. 111, 053603 (2012).
	SKip3
	Schliesser, A., Rivière, R., Anetsberger, G., Arcizet, O. & Kippenberg, T. J. Resolved-sideband cooling of a micromechanical oscillator. Nature Phys. 4, 415 (2008).
	
	
	SWick1
	Hübener, R., Mari, A. & Eisert, J. Wick’s theorem for matrix product states. Phys. Rev. Lett. 110, 040401 (2013).
	
	SWick2
	Isserlis, L. On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables. Biometrika 12, 134 (1918).
	
	SWick3
	Wick, G.C. The evaluation of the collision matrix. Phys. Rev. 80, 268 (1950).
	



