Exploring corrections to the Optomechanical Hamiltonian

We compare two approaches for deriving corrections to the “linear model” of cavity optomechanics, in order to describe effects that are beyond first order in the radiation pressure coupling. In the regime where the mechanical frequency is much lower than the cavity one, we compare: (I) a widely used phenomenological Hamiltonian conserving the photon number; (II) a two-mode truncation of C. K. Law’s microscopic model, which we take as the “true” system Hamiltonian. While these approaches agree at first order, the latter model does not conserve the photon number, resulting in challenging computations. We find that approach (I) allows for several analytical predictions, and significantly outperforms the linear model in our numerical examples. Yet, we also find that the phenomenological Hamiltonian cannot fully capture all high-order corrections arising from the C. K. Law model.

Most experiments to date have explored the weak coupling regime, in which radiation pressure effects are only visible upon strong driving of the cavity, and the system can be described as a pair of coupled harmonic oscillators 1 . Several experimental platforms, however, are now approaching the single-photon strong coupling regime 2-4,16-18 (strong coupling, or SC, for brevity), in which the anharmonic nature of the radiation pressure interaction must be taken into account. Strong coupling occurs when the the coupling rate of a single photon is greater than the typical loss rates of the setup. Such regime features clear departures from classical behaviour 11,16,19 , and facilitates the production of nonclassical states in both light and mechanics [20][21][22] . A further appealing feature of strong coupling optomechanics is that the Hamiltonian can be diagonalized analytically upon linearising the cavity frequency around the origin, as per ω ω ω + ′  x x ( ) (0) ( 0) (equivalently, this may be seen as a first order expansion in the coupling constant). We shall call linear model the resulting Hamiltonian, arguably the most widely used modelling tool in SC optomechanics.
Despite the undeniable success of the linear model as a theoretical device, the need is arising to go beyond this approximation. For example, it was noted by Brunelli et al. 23 that the linearized optomechanical Hamiltonian is unbounded from below, with (unphysical) negative energies cropping up at very high photon numbers. While this pathology may be seen as somewhat mathematical (see Appendix A), it highlights that the model must be eventually refined, particularly in view of future experiments aiming to achieve ever larger coupling strengths. Even in systems that are far from the SC, the quest for ultra-precise measurements (e.g. those pertaining Planck scale physics) 14 , or for the detection of dynamical Casimir effects 24 , demand a more accurate Hamiltonian description of quantum optomechanics.
Armed with these motivations, in this paper we explore and compare two different starting points that may be used to go beyond the linear model: (I) a widely used phenomenological Hamiltonian, which conserves the cavity photon number (phenomenological approach); (II) a two-mode truncation of C. K. Law's microscopic Hamiltonian for an optomechanical Fabry-Perot cavity 15 (microscopic approach). We shall take this truncated Law Hamiltonian as the benchmark (or "true model") against which the various generalizations of the linear model should be judged. In this paper we focus on the realistic situation where the bare cavity frequency is much larger than the mechanical one. We find that approaches I and II agree at first order in the coupling, but at higher orders the microscopic model yields 'counter-rotating' terms that violate photon number conservation 24 and make computations challenging. We resort to numerical diagonalization of the Hamiltonian, truncated in Fock space, to deal with such situations (and also to tackle the exact phenomenological Hamiltonian). The numerical examples we explored suggest that the phenomenological approach does a good job in improving the linear model in the typical parameter regimes of optomechanics experiments, yet it does not fully capture the second-order corrections arising from the microscopic treatment. We thus conclude that the number-conserving Hamiltonians considered in this paper can provide meaningful refinements to existing optomechanical models, arguably good enough to describe experiments in the near future without adding an excessive computational burden. Yet, our results also highlight that photon-number conservation must be eventually given up to model cavity optomechanics with ultra-high accuracy, even in parameter regimes where dynamical Casimir physics would be usually discounted.
The paper is organised as follows. In Sec. 2 we review the linear model of optomechanics. In Sec. 3 we show how this model can be refined phenomenologically, while in Sec. 4 we discuss model corrections that are microscopically motivated. Sec. 5 features numerical examples comparing the various approaches. We draw our conclusions in Sec. 6, and in Appendix A we discuss in more detail the negative-energies pathology of the linear model 23 .

Linear Model: Brief Review
We begin by recalling the model that is commonly taken as a starting point in theoretical optomechanics. To lighten the notation we take  = 1, and consider a mechanical mode of unit effective mass. The linear model Hamiltonian reads where â is the annihilation operator for the cavity field mode, x and p are respectively the canonical position and momentum operators of the movable mirror, Ω is the bare mechanical frequency, ω 0 is the bare cavity frequency, and the coupling constant G is known as frequency pull parameter. The canonical commutators read , while all field operators commute with all mirror operators. In passing we note that the optomechanical interaction strength is often quantified via the vacuum optomechanical coupling strength g 0 = Gx zpf , where = Ω x 1/ 2 zpf is the position uncertainty in the (bare) ground state of the mechanical oscillator. Loosely speaking, the dimensionless ratio g 0 /Ω quantifies the influence a single photon can have on the mirror motion. That being said, many of our expressions will have a more compact form when written in terms of G. We assume that the hierarchy ω Ω   g 0 0 holds between the model parameters, which is the case in many experimental realizations. Hamiltonian (1) can be diagonalized exactly 16,20 . The properties of the corresponding eigensystem are well known, and can be obtained from the section: Approach 2: standard quadratic expansion by setting β = 0.

Phenomenological Approach
A first possibility in extending the linear model, often adopted in the literature, is to consider a phenomenological Hamiltonian similar to Eq. (1), but where the cavity frequency has a more general position dependence ω(x). Crucially, no change is assumed in the commutation rules for mirror and field operators, an assumption that drastically simplifies computations but may not be justified microscopically 15 -we shall elaborate on this point in the Section: Microscopic approach.
For the purposes of illustration, from now on we consider an optomechanical setup consisting of a Fabry-Perot resonator with one movable mirror (see Fig. 1), which amounts to setting ω(x) = ω 0 /(1 + x/L) with L the bare cavity length. However, the reader should keep in mind that the same formalism is applicable when ω(x) is a generic (positive) function. The phenomenological Hamiltonian is Figure 1. Schematic of an optomechanical Fabry-Perot resonator with one fixed (right) and one movable (left) mirror. Upon reflecting a photon, the movable mirror recoils due to radiation pressure. At the same time the total length of the cavity (L + x), hence the photon frequency ω(x), is modulated by the mirror coordinate x. The movable mirror is subject to a trapping potential (depicted as a spring), which is usually harmonic to a good approximation. If we now pick L = ω 0 /G = (ω 0 /g 0 )x zpf and take a first order expansion in G (or equivalently g 0 ), Eq. (1) is recovered. The natural next step is then to retain higher order corrections in G. In what follows, we discuss two inequivalent approaches that employ simple harmonic oscillator mathematics to approximate Eq. (2) beyond first order. Both methods are in principle amenable to analytical treatment.
n n phen 0 2 featuring the effective n-dependent mechanical potential The idea is now to perform a quadratic expansion of V x ( ) n around its minimum: n n n n n n 2 where the equilibrium position x n is implicitly determined by Note that the above procedure is well defined, since it is easily verified that V n has a unique global minimum at fixed n, and that for all n, . The outlined procedure returns the effective Hamiltonian n n n n n n where we have chosen the symbol Ĥn to signify that the photon number operator =ˆ † n a a controls which harmonic oscillator Hamiltonian H n is acting on the mirror. For each n, the eigenvalues and eigenvectors of H n may be found via straightforward harmonic oscillator techniques. For example the energy eigenvalues, labelled by the cavity photon number n and a mechanical 'phonon' number m (both semipositive integers), read These are all positive, and thus may be used to calculate the partition function without the need for artificial cutoffs 23 . A fully analytical treatment of Ĥn requires the explicit inversion of Eq. (6), which may be recast as This equation for x n is in principle solvable (it is equivalent to a third-degree polynomial equation) but leads to very cumbersome expressions. A number of observations, however, can be made without resorting to the explicit solution. First, it is evident that it must be > x 0 n , and a little more thought reveals that → ∞ x n when n → ∞. The latter limit is of course a purely mathematical result, indicating that the effective model becomes unreliable when n gets extremely large. Yet, at low photon numbers this model is very accurate at predicting mechanical uncertainties in the energy eigenstates (see Section: Numerical examples). As a second observation, notice that Eq. (10) describes the intersection between a strictly increasing function and a decreasing one, such that for any n there is a unique real solution for x n , as anticipated. Third, implicit differentiation in n shows that x n is a strictly increasing function of the photon number (as intution about radiation pressure would suggest). Finally, combining Eqs (8) and (10) the photon-dependent mechanical frequencies may be expressed as this result is unlikely to have any physical significance: the above model eventually becomes unreliable when too many photons occupy the cavity (as discussed, the associated limit → ∞ x n means that the average cavity length becomes infinite!). The eigenstates corresponding to E n,m take the form where the mechanical wavefunctions φ n, m (x) = 〈x|φ n,m 〉 are expressible in terms of Hermite functions centered at = x x n (details not shown). The following eigenstate properties are easily verified: nm n m n m n , , In light of the above discussion we gain the following intuition about the structure of the Hamiltonian eigenstates in Eq. (12): the mechanical wavefunctions φ n,m (x) are pushed further away from the fixed mirror as the photon number n is increased, while at the same time they become more and more squeezed in position. Such behaviour is schematised in Fig. 2. Approach 2: standard quadratic expansion. Our second approach consists in a more traditional second-order expansion of Eq. (2). Specifically, we expand where β 2 = 2ω 0 /L 2 = 2G 2 /ω 0 . For brevity, we shall refer to the resulting Hamiltonian as the quadratic model: While the 'quadratic coupling' has been studied before in optomechanics 25-28 , the focus has often been on the striking differences between a purely linear (β = 0) and a purely quadratic (G = 0) interaction. Here, on the other hand, our primary goal is to improve the accuracy of our Hamiltonian model: we want to explore the corrections to the linear model due to the higher order β term. This is in a similar spirit to ref. 14 (7), with parameters n n 2 so that the resulting energy eigenvalues read valid for all semipositive integers n, m. Also here we shall discuss the limit n→∞ as a matter of mathematical curiosity. Note that the quadratic model predicts a finite average mechanical displacement even in the limit of large n, while the mechanical frequencies are divergent with asymptotic scaling Ω ∝ n n , in turn implying that the mechanical wavefunctions become infinitely squeezed in position:

Microscopic Approach
In this section we take as a starting point the microscopic Hamiltonian derived by C. K. Law through canonical quantization of a Fabry-Perot cavity with a movable mirror 15 . For the sake of numerical tractability we shall assume that only two degrees of freedom, one optical and one mechanical, are sufficient to model our optomechanical system. More precisely: we shall retain only one optical degree of freedom in Eq. (3.5) from ref. 15 . Furthermore, we shall not perform the gauge-like transformation used in Section IV of the same paper: the latter is a useful theoretical tool at first order in the coupling, but does not easily generalise to include the higher order effects we are trying to capture. Under the described assumptions, the Law Hamiltonian simplifies to where ΠQ, are canonical quadrature operators for the cavity field mode and ωx ( ) retains the same form as before. The elementary canonical commutator for the cavity field is now Π =Q i [ , ] , and again all field operators commute with the mechanical variables x and p. Note that Hamiltonian (23) leads to the Heinseberg equation i.e., the instantaneous radiation pressure force is proportional to the intensity of a single field quadrature (rather than to the field energy). This makes sense: since the (transverse) electric field E must vanish on the mirror, the radiation pressure force will be proportional to the squared magnetic field |B| 2 alone (evaluated at the mirror surface). For completeness, the remaining equations of motion are where the last two equations may be seen as an approximate reformulation of the wave equation for a single-mode cavity field of frequency ωx ( ). Hence, within the two-mode approximation, we argue that Eq. (23) should be a more reliable description of radiation pressure physics than Hamiltonian (2). In the next section, Eq. (23) will play the role of benchmark, or "true model", against which our various approximations shall be assessed. In the present section, we instead discuss the corrections to the linear model that arise in this framework, highlighting the issues that were not present in the phenomenological approach. In passing we remark that a microscopic Hamiltonian, alternative to that of C. K. Law, has been recently proposed for optomechanics 29 . In future studies, our analysis could be generalised to deal with this alternative theory, and/or to quantitatively assess the robustness of the two-mode approximation that we are currently taking for granted. Note that Hamiltonian (25) does not yield photon number conservation in general. More precisely, the concept of photon itself relies on identifying a reference frequency, for which a natural choice is ω 0 ≡ ω(0). We then choose the following definition for the cavity annihilation operator: With the above definition, the annihilation operator â commutes with all mirror operators. This allows us to represent the pairs ˆx p ( , ) and ˆˆ † a a ( , ) in different Hilbert spaces, as is done for the phenomenological models of Sec 3. We found this approach convenient to set-up numerical calculations based on the truncation of Fock spaces. In passing we note that this was not the path taken in the seminal paper by C. K. Law 15 , where instead xdependent annihilation operators were introduced to describe the cavity field. Hamiltonian (23) may be recast as ) /2 , ( ) ( ( ) )/2 2 0 2 0 2 0 2 0 and the cavity photon number ˆ † a a is evidently not conserved due to the appearance of counter-rotating terms (i.e. +ˆ † a a 2 2  where the linear model (1) has been recovered through a rotating-wave approximation, justifiable in the parameter regime ω Ω   g 0 0 that we are considering here. More in detail, a back-of-the-envelope calculation shows that the counter-rotating terms only provide a contribution at order 0 -for example this can be proved by time-averaging the Hamiltonian on timescales ω 1/ 0 30 . In the parameter regimes of interest we can thus reach the reassuring conclusion that the microscopic and phenomenological approaches yield the same linear model. We should however point out that there can be situations where ω 0 and Ω have comparable magnitude, and in that case the +ˆ † a a ( ) 2 2 term cannot be neglected even at first order in G 24 . It is at second order in G that non-negligible discrepancies arise between the microscopic and phenomenological approaches. Following similar lines of reasoning as above, we find that the second-order microscopic Hamiltonian is where we recall β 2 = 2ω 0 /L 2 = 2G 2 /ω 0 , and we note two additional terms that were not present in the phenomenological model and cannot be neglected at the considered order. Due to the counter-rotating terms, Hamiltonian (27) presents the same computational challenges as the original model in Eq. (23) -albeit its numerical solution appears to be slightly faster in the examples we explored. The relevance of Eq. (27) to our discussion is more conceptual than practical: it highlights that the phenomenological approach of the Section: Phenomenological Approach, when taken beyond first order in G, may begin to miss some of the finer details of optomechanical physics. In particular, the appearance of counter-rotating terms provides a qualitative departure from the phenomenological model. Luckily, the examples to follow indicate that the analytical models derived via the phenomenological approach may still provide significant improvements to the linear model in the regime ω Ω   g 0 0 which is relevant to many experiments. Indeed, it can be shown that the corrections to H quad appearing in Eq. (27) cancel each other out on average over timescales ω 1/ 0 -see 30 for details on the time-averaging procedure. This may explain why the quadratic model (16) does improve the linear model, even though it does not fully capture all second order effects.

Numerical Examples
So far we have introduced several Hamiltonians that aim to improve the linear model. It is now interesting to compare these refined models in a concrete example: first by examining a few observables that characterize their low-energy eigensystems, then by looking at their predictions for the time evolution of mirror position and field amplitude. Since we are considering small deviations from the predictions of H lin , it is practically impossible to compare by eye the plots of absolute quantities. Instead, for any quantity A calculated via a given approximate model, we shall plot the absolute error |A − A mic |, where A mic is computed by solving H mic numerically. Note that numerical solutions are adopted also for H phen and H mic (2) . For all these cases we wrote a simple Mathematica script to diagonalize the Hamiltonian in a truncated space and hence calculate all quantities of interest. To confirm the solidity of these results we performed the time-evolution calculations (Fig. 4) also in Python, which resulted in visually identical plots. For our purposes truncations of n max = 20 photons and m max = 30 phonons were sufficient to obtain convergence: our plots did not show visible changes when considering nearby truncation numbers of n max + 1 and m max + 1. Yet, note that these numbers rely heavily on the chosen model parameters, and on the number of reliable eigenvalues one is seeking to obtain.
In the displayed examples we set Ω = 100g 0 , inspired by recent experiments featuring , while the bare cavity frequency is fixed at ω 0 = 100Ω. Qualitatively similar results were obtained by increasing ω 0 to 10 3 Ω. With our current numerical means we found it difficult to obtain reliable results when increasing ω 0 further: the handling of such cases required prohibitively high numerical precisions. It may be useful to overcome this limitation in the future, in order to address systems where many orders of magnitude separate mechanical and optical frequencies.
Low-energy eigensystems. To compare eigensystems, we focus on the following quantities that depend on the photon number n: (i) Lowest energy eigenvalue, E n,0 (ii) Average mechanical displacement, 〈 〉 x n,0 (iii) Mechanical uncertainty, Δx n,0 Note that the eigensystems of Hamiltonians , n lin p hen (2) are fully characterized by these quantities. In contrast, n is not a good quantum number for H mic and H mic (2) , since these Hamiltonians do not commute with n. Yet, in the parameter regimes here considered we were able define an effective photon number n as the integer SCiENTifiC REPORTS | (2018) 8:9157 | DOI:10.1038/s41598-018-26739-0 closest to = 〈 〉 † n a a , since − n n remained small in all the reported examples (while n is obviously an integer, note that n may not be so in a non-number-conserving Hamiltonian). In turn, this allowed us to identify the lowest eigenvalue at fixed n, E n,0 , as well as the associated eigenstate |Ψ n,0 〉. From these we were able to define, and calculate, quantities (i)-(iii) for the models H mic and H mic (2) . Figure 3 shows that all the approaches considered here provide a significant improvement to the linear model, but at this stage there is no clear winner: while H mic (2) is better at predicting energy eigenvalues, the phenomenological Hamiltonian H phen can predict mechanical equilibrium positions and (most) variances more accurately. This suggests that it may be worth keeping all these Hamiltonians in one's theoretical toolbox for the sake of future optomechanical studies. In passing, it is interesting to note that Ĥn can be seen as a sort of hybrid model emulating properties of both H quad (equilibrium positions) and H phen (mechanical variances).

Dynamics.
To compare dynamical predictions we take the initial state a b Figure 4. Performance of the different models in predicting low-energy dynamics in the system. As in Fig. 3, we consider parameters ω 0 = 10 2 Ω, Ω = 10 2 g 0 . The dynamics of 〈 〉 a and 〈 〉 x are calculated for the initial state (28)  where n n ( / !) a n n a /2 2 is a coherent state and |0〉 b is the bare ground state of the mechanics. It is apparent that the evolution of initial state (28) is going to involve several photon number sectors even with the number-conserving Hamiltonians. Moreover, according to the linear model, the evolution of this state can induce interesting nonclassical features in the system, such as Wigner-Function negativity and opto-mechanical entanglement 20 . Therefore, we think that the initial state (28) embodies an interesting test-bed for our various approximate Hamiltonians: a reliable description of the subsequent dynamics must capture well the evolution of self-correlations in the cavity field together with entanglement between the two subsystems. In our examples we picked α = 2, and looked at the accuracy of the various models in predicting the time-dependent averages 〈 〉 x t ( ) and 〈 〉 a t ( ) : the relevant results are displayed in Fig. 4. From these we infer that H mic (2) has an edge on the other models, as it may be expected from its construction as the 'correct' second-order expansion of H mic . As a corollary, we deduce that the inclusion of counter-rotating terms will eventually become necessary in trying to model optomechanical system with ever higher accuracy. On the other hand, the good news is that once again the phenomenological approach provides a significant improvement to the linear model (except for a few isolated points in time, where the linear model yields better estimates for 〈 〉 x t ( ) ). Furthermore, in the chosen example we do not observe significant differences between the dynamical predictions of H phen , Ĥn and H quad . Under these circumstances, our choice of model for further analysis would naturally fall on H quad due to its simple analytical solution.

Conclusions and Outlook
We have proposed and compared several Hamiltonians that go beyond the linear model of optomechanics, featuring corrections that are either phenomenologically or microscopically derived. Using a two-mode truncation of Law's Hamiltonian 15 as our benchmark, we were able to confirm that the phenomenological approach, common in the literature, is able to provide meaningful corrections to the linear model. A particularly reassuring observation is that, in the regime ω Ω   g 0 studied here, these corrections are well captured by the analytically solvable model H quad , i.e., a straightforward second-order expansion of the phenomenological Hamiltonian (2). However, we also pointed out that the phenomenological approach does not fully capture all second order corrections arising from a microscopic treatment of radiation pressure. In particular we have outlined an important effect, the breakdown of photon number conservation, which can provide non-negligible corrections even in parameter regimes where dynamical-Casimir effects are usually discounted.
Our work may benefit future experiments probing the single-photon strong coupling regime or aiming for high-precision optomechnical measurements. More importantly, it opens a discussion that may uncover new avenues for the theoretical modelling of optomechanics, beyond the domain of applicability of the linear model. There are several important questions that our work leaves open and that should be addressed in future studies. First, our analysis is confined to a two-mode scenario. In the quest for strong coupling and/or high precision modelling, a multi-mode treatment 31 (including several degrees of freedom for both light and mechanics) may eventually become necessary. Second, we have left system-environment interactions out of our discussion. Due to the presence of three very different energy scales (ω 0 , Ω, g 0 ), we anticipate that it will be challenging to develop an accurate open-system model for optomechanics that goes beyond the existing phenomenological treatments. Yet, identifying a reliable Hamiltonian is a crucial first step towards the rigorous modelling of an open quantum system 32 . Third, we have assumed the bare mechanical potential to be harmonic. Depending on the implementation, anharmonicity in the mirror potential 33 may provide new corrections that can be qualitatively different from those studied here. Finally, a model including relativistic corrections may be eventually needed. The theoretical challenge will increase significantly in trying to address any of these limitations of our work, yet it will be a worthwhile pursuit in the quest to model ever more sophisticated optomechanics experiments.
Availability of materials and data. The original Mathematica and Phython codes utilised to produce the plots in Fig. 3 are available upon request.