Higher-order interactions in quantum optomechanics: Revisiting theoretical foundations

The theory of quantum optomechanics is reconstructed from first principles by finding a Lagrangian from light's equation of motion and then proceeding to the Hamiltonian. The nonlinear terms, including the quadratic and higher-order interactions, do not vanish under any possible choice of canonical parameters, and lead to coupling of momentum and field. The existence of quadratic mechanical parametric interaction is then demonstrated rigorously, which has been so far assumed phenomenologically in previous studies. Corrections to the quadratic terms are particularly significant when the mechanical frequency is of the same order or larger than the electromagnetic frequency. Further discussions on the squeezing as well as relativistic corrections are presented.

Usually, the analysis in all these above works is ultimately done within the linearized approximation of ladder operators, thus being limited in accuracy where ℍ OM interaction is non-existent or vanishingly small. For applications where quadratic or even quartic effects are primarily pursued, 0 may be designed to be identically zero [55][56][57][58][59][60][61], which urges need for accurate knowledge of higher-order interaction terms. Similar situation also could arise in trapped ultracold atomic gases [62], where linear interactions identically vanish. Relevant Physical phenomena in optomechanics such as four-wave mixing, also is suitably described by higher order interaction terms [63]. Moreover, significance and prominent role of such nonlinear interactions has been observed in few recent experiments [64,65].
A careful review of the theory of this subject [5], however, reveals that there are a number of physical approximations in formulation of the problem such as the non-relativistic limit [66,67], which makes the unified description of relativistic photon momenta and non-relativistic mirror motion inaccurate. As it is being shown here, a full treatment of the latter will yield higher-order multi-particle interactions. Such quadratic interactions have been recently used phenomenologically [58] without a theoretical basis. Similar yet weaker interactions may be also drawn from relativistic corrections [68] as discussed here. In that sense, the quadratic interactions are shown to receive contributions from both non-relativistic and relativistic terms, which become quite significant when the mechanical frequency is comparable or larger than the electromagnetic frequency.

Classical Hamiltonian
The basic theory to be discussed is based on two important and basic assumptions. We first assume that the cavity mode decomposition is valid independent of the mirror motion. Second, the electric fields vanish at the mirror surface in the frame that the Lagrangian is constructed. These assumptions seem reasonable in the usual discussions with much lower mechanical frequencies when the cavity mode change can be treated adiabatically. However, when the end mirror oscillates at a very high frequency, it dresses the cavity modes so that the mode frequency might become undefined invalidating the assumption of mode decomposition. At higher frequencies, it has also been known and demonstrated that the mirror undergoing relativistic motion could produce photons from vacuum, known as the dynamical Casimir effect [69,70]. This is of course beyond the regime being explored in this article.
The focus of the first two subsections 2.1 and 2.2 is to assert the claim, and demonstrate what term is missing and why it happens. As it will be shown and rigorously proven, even for the simplest case of interaction with a single-optical mode, a new term of the type ̇2 2 representing quadratic momentum ̇ and optical field interactions is found, the origin of which is also identified. For the more general case of multi-mode optical fields, the situation is even much more complex and there are a few more missing terms to consider. Once the Lagrangian is known, the Hamiltonian is subsequently constructed in the subsection 2.3.

The Equations of Motion
The one-dimensional wave equation for transverse component of the magnetic potential ( , ) in the dimensional form is expressed as [5] 2 ( , ) = ( , ), (1) where and are respectively the position and time coordinates, and is the speed of light in free space. Suppose the Fourier series relations for the magnetic vector potential are defined as [5] where is the cross-sectional area, 0 is the permeability of vacuum. This arrangement ensures that the definition of canonical variables can be used later, so that ̇2 simply takes on the dimension of energy.
One may furthermore define the functions = √2⁄ sin( ) and = √2⁄ cos( ), and hence ( , ) = ∑ where = √ 0 ⁄ . Here, the inner product is also defined as ( | ) = ∫ 0 such that the following relations may be found After straightforward calculations one obtains Here, the anti-symmetric coefficients and ℎ are Therefore, the wave equation reads Now, the inner product relationships (3) and some further simplification (Appendix A) help us to obtain the equation of motion as where all summations are nonzero only for ≠ , and The related equation in Law's paper [5] is completely equivalent, but apparently different in presentation. Similarly, one may directly deduce from the Newton's equation of motion

Lagrangian
The associated Lagrangian which leads to the above set of Euler equations is given by It should be again noticed that ℎ = 0 by (5), which together with the antisymmetry of coefficients from (5), allows further simplification to reach where = 1 2 (ℎ + ℎ ) + is related to the symmetric part of ℎ , given by This is to be compared with the quite different expression by Law [5] given by which indeed consistently satisfies the Euler's pair of equations. Now, it is a matter of speculation whether (9) is equivalent to (11) or not. Firstly, the diagonal terms are equal if Now, for a single-mode system with only one radiation mode, we may easily notice that (9) and (11) become readily identical. This becomes more clear by noticing that = ℎ = 0, and hence for a system with only one electromagnetic radiation mode, (9) becomes ̈= − 2 + ( −2̇2 ) , with = 1 . This is while Law's expression (11) is ̈= − 2 + −2̇2 (∑ 1 1 ) . But (17) requires that = ∑ 1 1 , which confirms the equivalency for single-mode systems.
As for a multimode system, and by comparing (13) and (16) one would need in order to (9) and (11) be identical. This can be put to numerical tests (here done by the author by coding a simple Mathematica program), and is in fact accurately satisfied. Hence, the single-mode Lagrangian in the non-relativistic limit can be written in the form The last term has been usually ignored so far in the literature, and will result in momentum-field coupling. In the remainder of the paper, we focus on nonlinear terms arising from this interaction, and then also add up the relativistic corrections in the end.

Hamiltonian
The definition of canonical momenta taken here is with and being some transformation coefficients to be determined later. One here may take advantage of the degree of freedom in choice of and to get rid of unwanted summation terms in the Hamiltonian. It has to be here noticed again that the existence of the last term of (20), being completely new even under the single-mode operation, has nothing to do with the choice of canonical momenta. That implies the final resulting single-mode (and therefore multi-mode) Hamiltonian will be inevitably different, incorporating a few new terms. The Hamiltonian may be now derived from the Lagrangian by iterated use of (21) and through the relationship ℋ =̇+ ∑̇− ℒ as Law [5] arbitrates the choice = = . But by going further with this choice for a single-mode optical field, it will be evident that =̇ and =̇, hence, still resulting in an extra nonlinear term proportional to 2 2 2 ⁄ in the Hamiltonian, leading to a fourth-order momentum-field interaction. This will be discussed shortly in the subsection below.
Hence, the Lagrangian ℒ found in the above yields the Hamiltonian below after some algebra Here is readily evident now by (19), that the last term can be made identically zero, only if . This not only cannot be satisfied by Law's choice = = , but also the second summation term nonlinear in ̇̇ will also survive, further complicating the Hamiltonian formulation. Now, further elimination of ̇ and ̇ from (21) gives Here, Dealing directly with such an intractable and long expression is without doubt too tough. Instead, we may tweak Law's choice slightly as = = 1 2 , which allows (23) be greatly simplified as This can be further eventually expanded and simplified as The last two terms of this Hamiltonian can be expanded to obtain multiple orders of interactions. These include higher-order tripartite phonon/two-photon, and quadpartite two-phonon/two-photon interactions, which does not exist in the Law's Hamiltonian [5], given by As it appears, (27) is missing two very different types of momentum field interaction as the last two summation terms of (26). This fact becomes evident in below.

Single Optical Mode
The interesting difference between these two Hamiltonians becomes quite clear with consideration of only one optical mode in the cavity. This simplifies our derived Hamiltonian to where ( ) = ⁄ and ≅ 3.8, while the Law's Hamiltonian [5] gives rise to the significantly different form As it will be shown below, (28) and (29) agree only to the first order, and hence up to the standard optomechanical Hamiltonian.

Field Quantization
When the obtained Hamiltonian is moved to the realm of quantum mechanics, it is first needed to define where (̂) = /̂ is defined according to (2). Also for a mechanical resonant frequency Ω and a spring restoring potential the displacement operator may be defined as ̂= +̂, where is the reference position of mirror, and hence the phonon ladder operators as with [̂,̂ †] = 1 and [ , ] = .
Now, it is necessary first to symmetrize [67] the classical Hamiltonian prior to insertion of operators, to ensure correct quantization of parameters. The process of symmetrization is done according to [71][72][73] etc. Therefore, after symmetrization the final form of the Hamiltonian is given by For a single-optical mode, (34) greatly simplifies and one gets This has to be applied to the last interacting term, which involves But symmetrization of a term which contains n non-commuting terms, results in n! terms, which for this case sum up to a total of 4! = 24 different expressions. The direct way to get around this situation is to first make an estimate of which terms are the strongest in the limit of linearized interaction and ignore the rest. It is possible furthermore to use the approximate replacement 1 ≅ to obtain where further substitutions should be taken as This can be decomposed to the terms ℍ = ℍ 012 + ℍ 3 + ℍ 4 + ℍ 5 + ⋯ , ℍ 012 = 1 2 ℏΩ 2 + ( ) + 1 2 ℏ (ℙ 2 + ℚ 2 ).
Hence, there are several distinct types of nonlinear optomechanical multi-phonon/multi-photon interactions, which by defining = /4 ≈ 0.95 are respectively given by and so on for higher order interactions. Here, the expansion of symmetrized terms, for instance, gives Now, it is noted that since usually ≫ Ω, it may observed that the first terms are much weaker than the second terms. Hence, using the identity ℙ 2 + ℚ 2 = where ℍ 3 ≡ ℍ OM is the simple optomechanical interaction, and ℍ 4 is known as the quadratic interaction. It has to be emphasized that while ℍ 3 ≡ ℍ OM is actually nonlinear in the exact mathematical sense, it is the quadratic interaction ℍ 4 which is mostly referred to as the nonlinear interaction in the literature [55][56][57]. Since it is possible to make 0 and therefore ℍ 3 identically vanish by appropriate optomechanical design [55][56][57][58] in which the overlap integral of optical and mechanical modes sums up to zero, hence the quadratic interactions ℍ 4 can then find physical significance. The quadratic interaction has been a subject of growing importance in the recent years in optomechanical systems [59][60][61] and beyond [62]. In [59] the photon statistics and blockade under ℍ 4 interactions has been studied and analytical expressions were derived. The quantum dissipative master function has been numerically solved and the corresponding correlation functions were obtained. Interestingly, quadratic optomechanical interactions can arise at the single-photon level, too, where rigorous analytical solutions have been devised [60]. Such type of interactions can be also well described using equivalent nonlinear electrical circuits, where a Josephson junction brings in the desired nonlinearity of quadratic interactions and terminate a pair of lumped transmission lines [61]. Finally, ultracold atoms also can exhibit interactions of a comparable type which is mathematically equivalent to the quadratic interaction [62].
The single-photon multi-particle rates are given by This summarizes the Hamiltonian as in which ℍ 0 and ℍ are respectively the non-interacting and interacting Hamiltonians when ≫ Ω. Now, the dimensionless constant is defined as with being the r.m.s. value of zero-point fluctuations, by which the following is deduced This implies that every kind of higher-order interaction is typically times weaker than the interaction of the preceding-order. It should be noted that while such interactions are normally expected to rapidly vanish with the order increasing, is a well-known fact that certain physical phenomena such as magnetism in solid 3He cannot be understood without inclusion of four-particle interaction terms [74,75]. It is worth here to mention that a detailed theory of optomechanics in superfluid 4He has been developed [76], but no expression for the nonlinear terms has been reported. In general, the interaction of mechanical and optical modes is not strictly one-dimensional, implying that the overlap integral of normalized modes should also be taken into account. For instance, odd mechanical modes with even optical modes have zero interaction. In that sense, tuning the interaction to an odd mode and then shining an even optical mode, or vice versa, makes the optomechanical interaction identically zero by setting ≡ √2 0 = 0. Then the lowest order surviving interaction would be the ℍ 4 term. This method has been used in [55][56][57] to highlight the quadratic interaction and make its measurement much easier. It has been shown that these quadratic terms may be exploited for direct observation of mechanical eigenmode jumps [55,56], as well as two-phonon cooling and squeezing [57], while the coupling strength could be increased by three orders of magnitude [56].
Moreover, the origin of mechanical parametric coupling which has recently been phenomenologically hypothesized [58] for the associated physical interactions cannot be understood without the presented analysis, although based on some earlier experimental evidence [77].
It must be added that the condition ≫ Ω may be violated in carefully designed superconducting microwave circuits and also the recently demonstrated molecular optomechanics [78], which signifies the importance of the 2 ℚ 2 term in ℍ 4 . It is furthermore worthwhile to point out that the regime = Ω can be indeed be accessed and investigated, as it has been shown experimentally for superconducting circuit optomechanics [79]. The proposal of light propagation in a cylinder with rotating walls [80] also requires accessing regimes where and Ω fall within the same order of magnitude. Alternatively, in situations where ≪ Ω, the scaling will be then given as which shows a significant enhancement in this type of interactions.

Conditions for Observation of Momentum-Field Quadratic Interactions
In summary, two general criteria should be satisfied in a carefully designed experiment to allow investigation of momentum-field quadratic interactions:  The optomechanical interaction ℍ 3 must vanish to allow easier study of quadratic interaction ℍ 4 . This is quite possible by design as extensively has been discussed in the above and literature [55][56][57][58][59][60][61][62].
 The mechanical frequency Ω must be of the same order of magnitude or exceeding the electromagnetic frequency . This is also possible and at least one experiment using superconducting optomechanics [79] has accessed this regime. Other possibilities are molecular optomechanics [78] as well as a rotating cylinder [80]. Evidently, such momentum-field quadratic interactions might be more difficult to observe under normal experimental conditions compared to the regular optomechanical setups. However, progressive developments in the precision and accuracy of optomechanics experiments, such as what happened for the case of Laser Interferometric Gravitational Observatory (LIGO) [81], could make it eventually possible to realize and probe such unexplored domains.

Optical Field
The standard method to linearize the interaction Hamiltonian can be now used by making the substitutions ̂→ ̅ +̂ where the new ̂ operator from now on stands for the non-classical perturbations and 〈̂〉 is a measure of optical field amplitude. Then ignoring higher-order terms and retaining only the lowest-order interacting terms, we get as well as Here, = ∡ ̅ and the coupling frequencies are defined as

Mechanical Field
Following the same method to linearize the mechanical motions, with the replacement ̂→ ̅ +̂ where the new ̂ operator denotes the perturbations, gives rise to the expressions ℍ 4 = +ℏ 4 + (̂ † +̂)(̂ † + −̂) , ≫ Ω, where = ∡ ̅ is set to zero without loss of generality, 4 + = 2| ̅| 4 + cos , and 4 − = 2| ̅| 4 − sin . In general, when ≫ Ω is violated, one would expect the momentum of mirror be coupled to the first quadrature of the radiation field. This type of interaction can be compared to the normal optomechanical interaction (50), in which the position is coupled to the first quadrature of the field.

Squeezing Hamiltonian
Without When the optical and mechanical frequencies do not differ by orders of magnitude so that neither ≫ Ω nor ≪ Ω hold, then the linearized Hamiltonian could be recast as This can be written as It is here to be noticed that ̂ and ̂ are in the standard form of Bogoliubov squeezing operator [54,55]. It may be noted that the equation is actually a function of by definition of 4 + and 4 − . Simplifying the above gives the expression for squeeze ratio as This shows that quadratic interactions give rise to squeezed mechanical or optical states unless = √ Ω and of course = 0.

Special Case
As discussed above, the Hamiltonian ℍ 3 can be made identically zero [56,57,82,83] to access the quadratic interaction terms ℍ 4 directly. There is an interesting condition on the ratio of optical to mechanical frequencies, which could be sought here. Let = √ Ω, in which is a constant to be determined later. This allows the ℍ 4 to be written as Then, for the choice of = 1 2 , that is ≅ 0.69Ω, one may reach the desired interaction quadratic Hamiltonian, linearized in the electromagnetic operators where the interaction rate is = = 2 | ̅|. When expanded in its four terms and after the replacement ̂= 1 2̂2 (Appendix C), one may immediately recognize the Hamiltonian of the type The first parenthesis represents the Hopping or Beam-Splitter term, while the second is normally referred to as the dissipation. Interestingly, the above could have been further linearized in mechanical operators to obtain where = 2 | ̅ |. This latter form, may find application in non-reciprocal optomechanics [84].

Relativistic Considerations
As a final remark, the approximate nature of the Lagranian formulation by Law [5] has not been left unnoticed. It could be attributed first to the non-relativistic description of mirror's motion which ultimately ignores higher-order interactions, and then to the relativistic nature of radiation friction force and the associated Doppler shift [66]. As a result, in a subsequent paper by Cheung and Law [67], it has been made clear that the non-relativistic optomechanical Hamiltonian is correct only to the first-order in .
The nature of the relativistic corrections can be quite different, as follows:  relativistic Doppler shift [66], which causes corrections in ⁄ ,  relativistic correction in radiation pressure term [67], the lowest-order of which being proportional to ⁄ ,  length contraction [68], due to the moving mirror boundary, resulting again in corrections as ⁄ .
Not surprisingly, all these relativistic terms vanish in the limit of infinite light speed . These altogether could be taken into account in a fully relativistic formulation of the Lagrangian and equations of motions for the mirror and optical field [85], which has been recently carried out in an extensive research by Castaños & Weder [68].
As shown in Appendix D, the total relativistic correction terms added to the Hamiltonian takes the form For the single-mode cavity, = 0 ℏ Ω 4 2 ⁄ , to be compared with = ℏ Ω 2 ⁄ in (48). Hence, the relativistic correction to the quadratic Hamiltonian ℍ 4 is Again, it is seen that when ≫ Ω is violated, the relativistic corrections might be quite significant. In any case, there is no relativistic correction to ℍ 3 .

Conclusions & Future Work
The derivation of the optomechanical Hamiltonian has been carefully examined from the modal expansions, equations of motion, all the way to the Lagrangian, and ultimately the Hamiltonian and relativistic considerations. A set of correction terms to the nonlinear terms have been identified, which do not eliminate under any choice of canonical momenta. With the careful system design which allows 0 = 0, these type of interactions are particularly interesting and now being actively pursued. It was shown that under these conditions one may expect coupling of mechanical momentum to the field position. Other sorts of interactions emerge under various conditions. In general, when the optical frequency is not much larger than the mechanical frequency, novel nonlinear interactions may appear.
A future work of the author [86] discusses a semi-analytical method based on the Langevin equations, which will enable easier study of quadratic interactions in quantum mechanical systems. This modified method of Langevin equations could in principle greatly simplify the study of quadratic and higher-order interactions, which would otherwise need either a full numerical solution using the master equation approach or full expansion unto the infinite set of orthogonal basis number kets.

Appendix A. Equations of Motion
In this section, we present the step-by-step details of the derivation of (9) from the previous equations, as it constitutes the most critical part of this article.
Starting from (8), one has first to rename the dummy index from to , multiply both sides by , and then take the inner product. This will yield the expression Using ( which by plugging in the definition for from (10) takes the form This is exactly the equation (9).

Appendix B. Special Case
Expansion of (60) results in In the limit ≫ Ω with → ∞, ℍ 4 as in (53) is recovered. With linearization of the electromagnetic field operators, and removing the non-interacting terms the following is found )̂] × (̂ † +̂). (B5)

Appendix C. Squared Annihilator
The operator ̂ has clearly a simple solution for its eigenkets, which is the same as coherent states such as | ⟩ where | ⟩ = 1 2 2 | ⟩. Hence, the eigenvalue is simply the complex number

Appendix D. Relativistic Correction
The relativistic Lagrangian density for a light field with normal incidence to a fully reflective and noncompressible moving mirror, correct to the first-order in ̇ and , reads [68] and is a dimensionless shape function independent of , being zero outside mirror and relative susceptibility of the mirror's dielectric 0 inside, and 0 is the permittivity of vacuum. By expanding in the powers of Β, this Lagrangian gives the first-and second-order corrections to the quadratic Hamiltonian density as It is appropriate to assume the approximation of conducting interface [87][88][89] for the mirror, such as the thickness is let to approach zero, while it susceptibility increases proportionally. In that limit, one may set where is the mirror's thickness. This is similar to the assumption of the locality of interaction by Gardiner & Zoller [90], too. Hence, one gets It is quite remarkable that (D8) is purely relativistic, and vanish in the limit of infinite , as shown below. Here, the dependence on is hidden for convenience. This term translates after symmetrization into This one after insertion of operators gives ∆ℋ (2) = − 1 2 ∆ℋ (1) and thus the total relativistic correction is found as ∆ℋ (1) + ∆ℋ (2) = −ℏ(̂ † −̂) 2 ∑ (̂ † +̂)(̂ † +̂). (D12)