Optomechanics with a position-modulated Kerr-type nonlinear coupling

Cavity optomechanics has proven to be a field of research rich with possibilities for studying motional cooling, squeezing, quantum entanglement and metrology in solid state systems. While to date most studies have focused on the modulation of the cavity frequency by the moving element, the emergence of new materials will soon allow to explore the influences of nonlinear optical effects. We therefore study in this work the effects due to a nonlinear position-modulated self-Kerr interaction and find that this leads to an effective coupling that scales with the square of the photon number, meaning that significant effects appear even for very small nonlinearities. This strong effective coupling can lead to lower powers required for motional cooling and the appearance of multi-stability in certain regimes.

Nonlinear effects in such systems are currently a frontier research topic, and it is well understood that optomechanical systems with a linear optomechanical interaction,Ĥ int ∼â †âq , can yield mechanical Kerr-type nonlinear optical responses [21]. Furthermore, optomechanical systems whose interactions are higher order in the position of the mechanical element, e.g.Ĥ int ∼â †âq2 have been studied [22][23][24] and shown to allow for new ways of coherent control of matter. Additionally a cross-Kerr interaction between the optical and mechanical modes, i.e.Ĥ int ∼â †âb †b , was considered in [25] and was shown to stabilize the bistable solutions as well as leading to tristable solutions in certain parameter regimes. In this work we will extend the range of nonlinear phenomena in cavity optomechanics by considering a situation where a nonlinear medium is present in the cavity. This situation has recently attracted some attention where a χ (3) material in the cavity was considered [26,27], and mechanical resonators with large optical nonlinearities are currently under development [28,29]. Another route to generate large optical nonlinearities is via coherent processes in atoms and the N -system of Schmidt and Immamoglu [30] is a well-known example. In our work, rather than considering a stationary optical nonlinearity, we will focus on the situation where the cavity's optical nonlinearity is directly modulated by the mechanical position, e.g. H int ∼â †2â2q . Such an interaction can be engineered * mathias.mikkelsen@oist.jp in superconducting systems, cavity polaritonic systems and atom-optical systems. For example, a giant-self-Kerr microwave optical nonlinearity was recently shown to be possible in a superconducting coplanar resonator via a capacitive coupling between two Cooper-pair boxes [31]. By utilizing an electromechanical capacitor (such as in Ref. [32]), a giant-Kerr mechanical modulation can then be straightforwardly engineered. Cavity polaritons, where light strongly couples to excitons in a quantum well or dot semiconductor hetrostructure, have also been studied for their use in optomechanics [33], and ultrastrong optomechanical couplings have been described [34]. The large interactions present in cavity polaritons are good candidates for optomechanical Kerr modulations to be expected [35]. Strong coupling between light and matter on the scale of individual atoms has recently been reported for atoms trapped near fiber tapers [36,37], and the generation of giant-Kerr optical nonlinearities using electromagnetically induced transparency techniques of N -systems (similar to Ref. [38]), would also lead to a position dependent self-Kerr interaction.
In this paper we will describe the general situation of nonlinear optical cavity optomechanics and include both, the usual optomechanical coupling and the Kerrtype coupling. For this we will first derive the relevant Langevin equations and calculate the classical steadystate solutions. We will then study the quantum fluctuations about these and the spectrum of the fluctuations. Using a numerical analysis we will discuss how the inclusion of the Kerr nonlinear interaction alters the cooling rates and other key figures of merit and end by analyzing the appearance of optical multistability and the situation of position dependent absorption.

II. THEORETICAL FRAMEWORK AND ANALYTIC RESULTS
In the following we will first outline the model and the theoretical framework we use. From this we will derive the analytic equations that govern the physical properties of the system. The system we consider can be described by an extension of standard optomechanical setups, such as the membrane-in-the-middle model [39][40][41], where the membrane is modeled as a harmonic oscillator with a single mechanical frequency Ω m and nondimensional position and momentum operators (q,p). The cavity field is modeled as a harmonic oscillator and described by the bosonic creation and annihilation operators,â † andâ, which create photons with the resonant frequency ω 0 . The membrane and the cavity interact via radiation pressure, which pushes the membrane and therefore changes the resonance frequency, making it position dependent ω(q). Since the harmonic oscillator is described byhω(q)â †â , this gives rise to a position-dependent interaction which is linear in the optical-field operators. While it is possible to perform the analytic treatment with no assumptions of the form of ω(q) [41], we will consider the case of from the beginning as this clarifies the relation between the two types of interaction we will consider. Within the linearity assumption the optomechanical interaction is then described asĤ L,int = −hâ †â g Lq . The second type of interaction we consider is nonlinear in the optical field operators, which corresponds to a simultaneous two-photon process that can be facilitated by a χ (3) material and is described by the term H Kerr = ηhâ †â †ââ . In order to obtain a positiondependent interaction with a moving membrane the nonlinear coefficient has to become dependent on its position η(q). If this can be engineered, then the nonlinear coefficient will be given similarly to the linear coefficient as Throughout the paper we will generally use the dimensionless position-coordinateq, but in order to derive the values of these coefficients for a specific physical system, it is easier to use the dimensional position coordinatê x, which is related to the nondimensional coordinate by some characteristic length scale x 0 asq =x/x 0 . Correspondingly, one can relate coupling strengths for the nondimensional coordinates to those of the dimensional coordinates by g NL = x 0 G NL , g L = x 0 G L . The characteristic length scale in optomechanics is generally the zero-point motion of the mechanical element x zp [1].

Physical implementation in an optical cavity
A model where the entire space between the two mirrors (where one of them can move) of a cavity is filled by a χ (3) medium giving rise to a termĤ Kerr was recently considered in [26]. We note that this model leads to a position dependence of η as and both, the resonant frequency ω(x) = ω 0 − G Lx and the cavity volume V c (x) = (L 0 +x)V c,0 /L 0 , depend on the position coordinate of the end mirror. Here 0 is the vacuum permittivity, while L 0 is the length of the cavity and V c,0 is the cavity volume, both in the absence of coupling. If we evaluate the first-order derivative at x = 0 and use the form of the linear coupling in this setup G L = ω 0 /L 0 , then we find that G NL = −3 η0 L0 . This means that a relatively big nonlinear coupling can be achieved only when the cavity length L 0 is small. Since x zp L 0 always holds, however, the interaction termĤ NL,int = −hâ †â †ââ g NLq will be small compared toĤ Kerr . The photon blockade effects arising from thê H Kerr term which was the topic of the investigation in Ref. [26] would therefore obscure the physics of interest in this work.

Physical implementation in a microwave cavity
In order to overcome this issue and to obtain a variable strong nonlinear interaction we turn to a different setup in which η is more tunable. It has been shown that a giant-self-Kerr microwave optical nonlinearity is possible in a superconducting coplanar resonator via a capacitive coupling between two Cooper-pair boxes [31]. While the details are too involved to present here (see Appendix), such a setup allows for a nonlinear coupling that depends on the mutual capacitance between the two Cooper-pair boxes. As the mutual capacitance C m depends on the distance between the plates, C m can be made dependent onx by coupling the motion of one of the Cooper-pair boxes to the physical motion of a membrane. By manipulating the parameters within physically realistic constraints, nonlinear couplings g NL of similar size to the typical linear couplings g L achievable in microwave cavities are obtainable (g L , g NL ∼ a few kHz). The same caveat as in the optical case is still present, but due to the tunability of the artificial molecules it is possible to place a second molecule inside the cavity which does not couple to the position and has the same magnitude as η 0 but the opposite sign. See Appendix for more details and a discussion of the relevant physical parameters. As it is possible to engineer situations where η 0 = 0, we will ignore the constant Kerr-term in our Hamiltonian, as it will generally lead to photon blockade which diminishes the effects of the nonlinear interaction that we want to investigate. Instead we consider just the position-modulated nonlinear term, described by the interaction Hamiltonian H NL,int .
Finally the light field is produced by an input laser which has a frequency ω L and an energy E. The full Hamiltonian in the frame rotating at the driving laser frequency ω L is then given bŷ Here the first term corresponds to the optical harmonic oscillator with detuning ∆ = ω 0 − ω L , while the second term corresponds to the single-frequency Ω m mechanical oscillator and the third term corresponds to the the input laser field of strength E. The last two terms account for the linear and nonlinear interactions as described above.

B. Quantum Langevin formalism
The optomechanical setup we are considering is an open system and interactions with the environment in the form of photon losses, mechanical dissipation, and noise have to be taken into account. We do this by utilizing the quantum Langevin formalism and consider photon loss associated with the end mirrors of the cavity κ 0 and photon loss associated with a moving element, such as a membrane κ 1 (q). While experimental evidence suggest that the dependence of κ 1 (q) on the position is usually small [42], we take it into account to get the most general picture. The losses have associated noise operatorsâ in 0 andâ in 1 , which have the correlation relation [41] â in In addition to the losses associated with the optical field, mechanical dissipation of the membrane γ m must be considered. For this the associated noise operator has the correlation relation where ν is the Fourier space frequency and is the mean phonon number at temperature T 0 and δ (t− t ) is the derivative of the Dirac-delta function. Adding these to the Heisenberg equations for the Hamiltonian (4) we arrive at the quantum Langevin equations (QLE) The ∂ q κ 1 (q) term is an effective noise term arising from the absorption in the membrane, as long as it is not assumed constant [41]. Solving this full set of quantummechanical equations is not a trivial task and we therefore employ a semiclassical approximation where we assume the physical operators can be expressed as a classical average {q, p, α} plus some small quantum fluctuation {δq, δp, δâ} This assumption only holds when the classical part is large, which is the situations experiments are commonly in. In the next section we will look at various aspects of the solutions for the system.

C. Classical steady-state solutions
To determine the classical steady-state solutions we consider the quantum fluctuations to be small compared to the classical c-numbers, so that any terms containing them can be neglected. Setting all derivatives in the QLE to zero then leads to where κ(q s ) = κ 0 + κ 1 (q s ). Assuming that κ 1 (q s ) = κ L q s these two coupled equations can be expressed as a single seventh-order polynomial for the mean photon number n s = |α s | 2 of the form Even though this seventh-order polynomial equation has in principle seven roots, all complex solutions can be dis-carded as the mean photon number has to be real. The steady-state position of the membrane can then be found by inserting the resulting solution for n s into expression (12). We note that doing a similar analysis for the linear system leads to a third-order polynomial only.

D. Linearised quantum Langevin equations: Effect of quantum fluctuations
To determine the stability of the steady-state solutions and to calculate physical values that depend solely on the quantum fluctuations, such as the temperature of the membrane, we will in this section look at the effects of the fluctuation terms. For this we use the quadratures of the 2i and insert the steady-state solution plus the quantumfluctuations into the QLEs. We are still considering the quantum fluctuations and the noise operators to be small and therefore only keep terms up to first order in these operators. This leads to the linearized quantum Langevin equations (LQLE) of the forṁ where we have defined the effective loss rate due to the position-dependent κ 1 (q) as Γ = √ 2∂ q κ 1 (q s )α s and the overall effective coupling in the system as G = √ 2α s (g L + 2g NL |α s | 2 ). Additionally we have defined the effective detuning ∆(q s ) = ∆ − g L q s . It is worth noting that the nonlinear coupling enters the effective coupling scaled with the mean photon number in the cavity and it can therefore be expected to have a much stronger effect than the linear coupling.
The LQLE can be rewritten as the matrix equation of the formu where and The drift matrix is given by and can be seen to reduce to the linear drift matrix for g NL = 0 [41]. The four eigenvalues appear as complex conjugate pairs and give information about the quantum fluctuations in the system. Their real parts corresponds to the cooling (or heating rate) of the membrane and the cavity, which means that the system is only stable if both of these are negative and the system relaxes towards a steady state. The imaginary parts describe the dressed eigenfrequencies of the membrane and the optical field. The stationary state is characterized by the covariance matrix V, which is determined by the Lyapunov equation and where D, known as the diffusion matrix, is related to the noise vector c(t) by [41] It is explicitly given by [41] The mean occupation number for the phonons n can be obtained from the mean energy of the mechanical resonator which is given by where V 11 and V 22 are matrix elements from the covariance matrix.
The spectral function for the system is defined as [26] S q (ν) = 1 4π dΩe −i(ν+Ω)t δq(ν)δq(Ω) + δq(Ω)δq(ν) Fourier transformation of the LQLE one finds where with δ (q s ) 2 = ∆(q s ) 2 + 12g 2 NL q 2 s − 8∆(q s )g NL q s and δ(q s ) = ∆(q s ) − 2g NL q s . Using the correlation relations for the noise operators in Fourier space then gives the spectrum where the thermal and radiation pressure spectra are given by and we have an extra noise spectrum associated with membrane absorption Here |χ eff | 2 us the effective susceptibility and is given by Again we note that the spectrum reduces to the linear one discussed in Ref. [41] if we assume that g NL = 0.

III. NUMERICAL RESULTS
To stress the effects stemming from the nonlinear nature of the coupling, we focus in the following on results that show significant differences to the linear situation. Furthermore, for the clearest comparison we initially make the approximation that κ 1 (q) = κ 1 , which reduces the problem to the standard optomechanical setup. This approximation is also consistent with existing experimental data for a linear membrane [42]. In general the behavior of our system can be distinguished into two categories. In the first one, g L and g NL have the same sign and thus enhance each other. In this regime the nonlinear coupling leads to a very strong effective coupling but no qualitative differences from the system where only a linear coupling is present. In the second category g L and g NL have opposite signs, i.e., one of them is attractive and the other is repulsive. Here a parameter regime exists for which additional steady-state solutions, that are not present for the linear system, can be found. These arise from the seventh-order polynomial describing the steady state of the nonlinear system. In our calculations we use a mechanical frequency of Ωm 2π = 356.6 kHz and a loss rate of κ/2π = 77 kHz corresponding to the experimental values in Ref. [42]. We generally employ a mechanical dissipation given by γ m = 0.01κ and choose a temperature which gives an initial phonon number of n 0 = 1.

A. Nonlinear enhancement
The phonon number as a function of the input energy of the laser and the linear and nonlinear couplings, respectively, for ∆ = Ω m is shown in Figs. 1(a) 1(b). This corresponds to the resolved sideband regime where optimal cooling can be expected. One can see that for any value of g L or g NL there is a value of E for which the same maximal cooling (n ≈ 0.01) is obtained and this value reduces as the coupling strength is increased. In fact, the nonlinear coupling simply enhances this behavior that is already present in the linear case, reducing the value of E required to achieve the same amount of cooling as in the linear case. Additionally one can see a rise in temperature as E increases beyond the value for maximal cooling, until the solutions become unstable (white areas in the figures). The photon numbers for the same parameter ranges are plotted in Figs. 1(c) and 1(d). One can see that they generally increase with energy as one would expect; however, it is more relevant to confirm whether the photon number stays above 1 in the regime of relevance, i.e., the regime of maximal cooling, since our model is otherwise inapplicable. In fact, for very small values of E it drops below 1 and to get maximal cooling for large values of g NL , values of E approaching this inapplicable regime must be chosen. For the most part, however, this limitation does not pose a problem. Situations where both (g L and g NL ) are finite behave similarly to the two limiting cases discussed here.
To investigate the existence of a strong-coupling regime at small laser intensities, we will look next at the spectrum S(ν) for (g L , g NL ) = (0.1κ, 0) and (g L , g NL ) = (0.1κ, 0.01κ). In both cases we have a fairly weak linear coupling and therefore do not expect the first case to be in the strong-coupling regime. The spectrum as a function of the detuning ∆ at E = 4Ω m is plotted in Figs. 2(a) and 2(b), and one can see that without the nonlinear coupling the system is indeed in the weak-coupling regime (no splitting is visible). However, adding a comparatively small amount of nonlinear coupling, g NL = 0.1g L , brings the system into the strong-coupling regime, signified by the avoided crossing of the dressed eigenfrequencies. This is quite remarkable as it allows for strong coupling at much smaller laser energies by adding just a small nonlinearity to the system. The spectrum as a function of the energy E at ∆ = Ω m , is plotted in Figs. 2(c) and 2(d), where the normal-mode splitting for the nonlinear coupling g NL = 0.1g L is distinctly visible and clearly absent for the linear case at these energies. The figures also shows that our results are self-consistent as the imaginary part of the eigenvalues of the drift matrix A (dotted white lines) correspond to the two peak positions in frequency space.
Another situation where the dramatic effect of a small nonlinearity becomes visible is shown in Fig. 3(a). Here the photon number is given as a function of E for (g L , g NL ) = (0.1κ, 0) and (g L , g NL ) = (0.1κ, 0.01κ) at ∆ = 50Ω m . The linear system can be seen to have only one solution in the displayed energy range, but adding the small nonlinearity leads to optical bistability as evidenced by the characteristic S-shaped curve. Starting at E = 0 and slowly increasing the energy the system moves along the first stable solution with the smallest number of photons. When this solution becomes unstable there is a first order transition to the second stable solution which has the largest number of photons. If the energy is then decreased the system stays in this second steady state until it becomes unstable at which point a first-order transition to the first steady state takes place. This is the characteristic signature of optical hysteresis. In Fig. 3(b) the corresponding phonon numbers are plot- ted. The temperature (phonon number) of the solutions are fairly constant in the stable regimes but rise rapidly as a solution becomes unstable.
It is worth stressing again that this optical bistability appears at quite low energies, because the nonlinearity leads to strong effective coupling. To see bistability with linear coupling only, much larger energies would be needed.

B. Optical multistability
While we have shown above that choosing the linear and the nonlinear coupling strength to have the same sign leads to mostly an enhancement of the effects already present in the linear case, qualitatively new effects can appear when they have opposite signs. For this we consider in the following (g L , g NL ) = (10κ, −10 −4 κ) and ∆ = 50Ω m . The reason for the large difference between the linear and the nonlinear coupling strength is that if the couplings were chosen to be of comparable size, the nonlinear coupling will dominate the system as it enters the effective coupling scaled with the photon number, i.e., G = √ 2α s (g L + 2g NL |α s | 2 ). Therefore, in order to engineer competition between the two forces, values for which the contribution of the nonlinear coupling g NL |α s | 2 and the linear coupling g L are comparable must be chosen. In Fig. 4 we show the behavior of |α s | 2 , n, q s , and g NL |α 2 s . Five steady-state solutions can be found, but only three of them are stable. To ease the discussion we will denote the three stable branches as 1, 2, and 3, where 1 corresponds to the stable solution with the smallest number of photons, 2 corresponds to the middle solution, while 3 corresponds to the solution with largest photon number [see Fig. 4(a)]. For increasing E, starting from E = 0 the system will be in a state on the first branch, jump to the second branch at energies where the first branch is no longer stable, before finally moving to the third branch. Looking at the position of the membrane for these states, one can see from Fig. 4(c) that it moves in the positive direction along the first branch, then changes direction and moves in the negative direction along the second branch, and, finally, jumps to a negative value moving in the negative direction when the system enters the third branch. This means that as the system jumps from the second to the third branch, the position of the membrane should jump from being displaced in the positive direction to being displaced in the negative one. Since q s = (g L + g NL |α s | 2 ) |αs| 2 Ωm one can understand this behavior by looking at the size of g L and g NL |α s | 2 , which is plotted in Fig. 4(d). Along the first branch g L , which is positive, is dominant and so the membrane is pushed in the positive direction. Along the second branch g NL |α s | 2 , which depends on the photon number, becomes comparable with g L and though q s stays positive the effective coupling G becomes negative. Increasing values of E mean an effective decrease in the (positive) position q s . On the third branch g NL |α s | 2 is always much larger than g L and therefore the membrane has a negative effective coupling, while also being at a negative value in position space. This means that the system moves between attractive and repulsive regimes for the interaction between the membrane and the cavity field as it jumps between stable solutions. Finally, looking at the phonon number in Fig. 4(b), one can see that the temperature (phonon number) of the second branch is considerably higher than that of the other two, with the third branch having the lowest temperature.

C. Effect of position-dependent absorption
Finally we consider the case of letting the membrane absorption depend on the position, i.e., κ 1 (q) = κ Lq with κ L = 0.0130κ 0 and κ 0 = 77kHz, and study how this affects the cooling in the resolved sideband regime, where ∆ = Ω m . For this, the phonon number n and the photon number |α s | 2 as a function of g NL and g L are shown in Fig. 5, where the energy at each point has been chosen so as to minimize n. One can see that the position dependent absorption coefficient leads to overall worse cooling [compare with n ≈ 0.01 from Figs. 1 (a) and 1 (b)], which is to be expected. However, by increasing the coupling this trend can be counteracted, leading to more efficient cooling for a more strongly coupled system. As can be seen from Fig. 5 this favors the nonlinear coupling, which is always stronger than the linear coupling. Therefore in this case the nonlinear coupling can lead to more efficient cooling. Looking at the photon number however, one must be wary as the energy for which optimal cooling is obtained corresponds to photon numbers on the order of one for very strongly coupled systems. While this means that the basic assumption for our model breaks down, we find enhanced cooling for the majority of the investigated ranges without encountering this problem.

IV. CONCLUSION
Applying the Langevin formalism we have derived the classical steady-state solutions and the spectrum of quantum fluctuations for a generic optomechanical system including a nonlinear position-modulated self-Kerr optomechanical interaction. We find that the the effective nonlinear coupling scales with the square of the photon number, which implies that even for a small nonlinear coupling the system enters the strong-coupling regime. By analyzing the obtained solutions numerically this is confirmed as a small nonlinear coupling leads to normalmode splitting and an avoided crossing in the spectrum, which is associated with the strong-coupling regime. This also leads to lower powers being required for motional cooling and bistability compared to systems where only a linear coupling is present. Furthermore, we find that the addition of a weak nonlinear coupling with the opposite sign of the linear coupling leads to three stable solutions, who can all be occupied as a function of the energy. Finally, we find that in the case of positiondependent absorption the nonlinear coupling can help counteract the degradation of cooling previously predicted in this regime. All these effects are consistent with a strong effective nonlinear coupling, even at small coupling strengths g NL . This means that engineering an effective nonlinear system experimentally is easier than what would initially be assumed, as only a small nonlinearity is required to observe dramatic effects.
Since it is known that nonlinear systems can support self-sustained oscillations and researchers have investi-gated these in the case of optomechanical systems possessing motional Kerr optical nonlinearities [43,44], we expect that with nonlinear optical coupling the behavior of such self-sustained oscillations may be even richer and an interesting topic for future investigations. Some possible avenues of experimental realization for such systems are cavity polaritons [35] and atoms trapped near fibre tapers [36,37], but the most promising one is utilizing artificial multilevel cooper pair box molecules in a microwave cavity [31], which we discuss in the appendix of this paper.
where 0 is the vacuum permittivity, r is the permittivity of the material between the plates, A is the area of both plates ,and d 0 is the initial distance between the plates. This means that any quantity, including η which depends on the capacitance C m , becomes position-dependent as well. In Ref. [31] it is shown that the nonlinear coefficient η in the regime where ( g1 Ω C ) 2 1 is given by where ∆ and δ are the detunings between the cavity frequency ω c and the fourth and second energy states |4 and |2 , respectively. The γ ij are the spontaneous decay rates from |i to |j , g 1 and g 2 are the coupling strengths between |1 , |2 and |3 , |4 , respectively, and, finally, Ω C is the coupling between |2 and |4 (see Fig. 1 in Ref. [31]). The position-dependence is contained in the two detunings δ and ∆ as these depend on C m . Analytic expressions for the these detunings can be found at the coresonance point where they become equal, δ = ∆, and are given by [31] ∆ whereh with k 1 = C Σ1 +C Σ2 and k 2 = C Σ1 C Σ2 . Here e is the elementary charge, while C Σ1 , C Σ2 are capacitances coming from different parts of the circuit (see. Fig. 2 in Ref. [31] for a diagram of the circuit). Depending on the setup, the resonant frequency ω c (x) can also be dependent on the the position coordinate (and in fact has to be if we wish to engineer a linear coupling). In order to simplify our estimate of the nonlinear coefficient, we will also assume that g 1 = g 2 = g and γ 21 = γ 23 = γ 43 = γ (similarly to Ref. [31]), which leads to To find the nonlinear coefficient we then need to evaluate the derivative of η atx = 0 (see Sec. II.A), which we can obtain analytically by simply taking the derivative of Eq. (40). In order to evaluate atx = 0 we require and To estimate the obtainable sizes of G NL and η 0 we need to input physically realistic values for the different parameters. We use estimates similar to those in Ref. [31], g/2π = 300κ, κ = 1 MHz, γ = 10 MHz, Ω C /2π = 0.9478 GHz. Additionally, we model the mutual capacitance, using values corresponding to those in the experimental setup from Ref. [2] with C 0 = 940 fF, d 0 = 50 nm and use a similar cavity frequency ω c /2π = 7.54 GHz as well as assuming a linear coupling equivalent to what they obtain, G L /2π = 49 MHz/nm. Finally, we use the values C Σ1 = C Σ2 = 4 fF [45]. It turns out that ω x is a useful knob for changing G NL , while keeping the other values constant, and one can get a large range of values for the nonlinear coefficients: For example, for ω x /2π = 3.3595 GHz we get G NL /2π = 95.3 MHz/nm, η 0 /2π = −0.751 MHz, while ω x /2π = 3.34 GHz gives G NL /2π = 0.637 kHz/nm, η 0 /2π = 0.1001 kHz. The nonlinear cou-pling can be tuned over several orders of magnitude and while much smaller values than the ones depicted above are easily obtainable, the main takeaway is that the nonlinear coupling can be made comparable in size to the linear coupling. To get the nondimensional coupling coefficients we multiply by the zero-point position from the experimental setup in Ref. [2] x zp = 4.1 fm which for the first case above corresponds to g L = 1.26 kHz, g NL = 2.46 kHz. Finally, we note that it is possible to place a second artificial molecule within the cavity, which does not couple to the position element, in such a way that the term η 0 is canceled out. This is important, as it allows the nonlinear interaction term to be engineered without the photon blockade effect associated with thê H kerr term. In order to engineer a value η of the same size with the opposite sign, the easiest way is to simply engineer a detuning with the opposite sign, i.e., ∆ stationary = −∆ moving (0), keeping the other parameters constant. By allowing the parameters δ, ∆, γ 21 , γ 23 , γ 43 to vary as well there are, however, many more ways to engineer an artificial molecule which cancels η 0 .