A Asymptotic behaviors

# Enhancing Quantum Effects via Periodic Modulations in Optomechanical Systems

## Abstract

Parametrically modulated optomechanical systems have been recently proposed as a simple and efficient setting for the quantum control of a micromechanical oscillator: relevant possibilities include the generation of squeezing in the oscillator position (or momentum) and the enhancement of entanglement between mechanical and radiation modes. In this paper we further investigate this new modulation regime, considering an optomechanical system with one or more parameters being modulated over time. We first apply a sinusoidal modulation of the mechanical frequency and characterize the optimal regime in which the visibility of purely quantum effects is maximal. We then introduce a second modulation on the input laser intensity and analyze the interplay between the two. We find that an interference pattern shows up, so that different choices of the relative phase between the two modulations can either enhance or cancel the desired quantum effects, opening new possibilities for optimal quantum control strategies.

## I Introduction

Theoretical studies and huge technological progresses over the last decades made it possible to reach a considerable level of control over quantum states of matter in a large variety of physical systems, ranging from photons, electrons and atoms to bigger solid state systems such as quantum dots and superconducting circuits. This opened the possibility for novel tests of quantum mechanics and allowed, among other things, to take important steps forward in investigating the quantum regime of macroscopic objects. In this perspective, one of the main goals in today quantum science is controlling nano- and micromechanical oscillators at the quantum level.

Quantum optomechanics (1); (2); (3); (4), i.e. studying and engineering the radiation pressure interaction of light with mechanical systems, comes as a powerful and well-developed tool to do so. First, radiation pressure interaction can be exploited to cool a (nano)micromechanical oscillator to its motional ground-state (5); this is a necessary step for quantum manipulation and could not be accomplished by direct means such as cryogenic cooling (at the typical mechanical frequencies involved this would require cooling the environment to a temperature of the order ). Backaction cooling has been experimentally demonstrated for a variety of physical implementations, including micromirrors in Fabry-Perot cavities (6), microtoroidal cavities (7) or optomechanical crystals (8). Second, there exists a strong analogy between quantum optomechanics and non-linear quantum optics, so that many (if not all) optomechanical effects can be mapped onto well-known optical effects. As a result, optomechanics becomes a natural way for controlling a mechanical resonator at the quantum level. Experimentally, the strong coupling regime needed to observe quantum behaviors has been demonstrated only very recently (7); (9), and detection of quantum effects is still awaiting. Nevertheless, a lot of theoretical studies on the subject has been carried out in the last decade and several proposals have been produced (10). These cover, among other things, the generation of entanglement between one oscillator and the radiation in a Fabry-Perot cavity (11), the generation of entanglement between two oscillators (12), or the generation of squeezed mechanical states (13); (14). In particular, references (14); (15); (16) introduced a new and effective way of enhancing the generation of quantum effects, which relies on applying a periodic modulation to some of the system parameters (a similar result has also been found in the analogous contest of nanoresonators and microwave cavities (17)).

In this paper we further investigate the properties of periodically modulated optomechanical systems and we address the following questions: which is the fundamental link between modulation and enhancement of quantum effects? is there an optimal choice of the modulation, for which the visibility of quantum effects is maximal? Is this optimal regime robust against parameter fluctuations? What happens when two independent modulations are applied simultaneously? To tackle these issues we analyze the paradigmatic case of a mechanical oscillator whose natural frequency is externally modulated when it evolves under the action of the noise and of the radiation pressure exerted by the photons of an externally driven optical cavity mode. While quantum optomechanics is nowadays extensively studied within a variety of experimental setups, the modulation of the mechanical frequency we analyze here is a very crucial aspect of our system and one that has not been implemented yet. However, very recent proposals for doing optomechanics with levitated dielectric spheres (18); (19); (20) can be a good answer. In these proposals the mechanical degree of freedom is represented by the center of mass motion of a nanodielectric sphere which is trapped and levitated by means of an optical trap. The sphere is then put inside an optical cavity, where it interacts with the intracavity radiation via the usual optomechanical Hamiltonian (2). The frequency of the center of mass motion depends on the shape of the trapping potential and can thus be modulated adjusting the intensity of the trapping laser, as shown in (18). Moreover, typical parameters attainable with such setups are comparable to those we have adopted in our simulations (see below), assuring the feasibility of the system under analysis.

In the above scenario we study the formation of squeezing, entanglement and discord (21), showing that in the steady state all these quantum effects are enhanced when the modulation frequency is twice the original value of . As we shall see, such resonance admits a simple interpretation in terms of an effective parametric phase-locking between the external driving forces and the natural evolution of the involved degree of freedom. Similar enhancements were also observed in Refs. (14); (15), where an harmonic modulation of was imposed on the amplitude of the cavity mode laser, and in Ref. (16), where a harmonic modulation of was imposed on the coupling rate between two generic bosonic modes. Since several mechanisms can lead independently to the same effect, an interesting question is how they can be best exploited to control specific quantum properties in the system. This goes in the direction of developing optimal quantum control protocols, a topic which is currently benefiting from many contributions (22). In the present case, to study the interplay of different mechanisms we add a second modulation in our model and we observe the arising of interference pattern in the system response. Specifically we notice that the ability in cooling and squeezing the mechanical oscillator strongly depends upon the relative phase of the two modulations, the relative variation being almost 50.

The material is organized as follows. In section II we present the system and solve its dynamical evolution under the action of a periodic modulation of the mechanical frequency. In section III we then characterize the asymptotic stationary state in terms of entanglement, squeezing, etc. In section IV we compare our findings to other recent proposals (14); (15) and we study what happens when a second independent modulation is applied to the system [specifically, in our case we introduce a modulation on the amplitude of the input laser]. Conclusions and general remarks follows in section V. Some technical derivations are finally reported in Appendix A.

## Ii The system

Our choice falls on the simplest optomechanical system of all, i.e. a Fabry-Perot cavity of lenght with a movable mirror at one end (see Fig. 1), which nevertheless captures all interesting physics. We can reasonably assume (10) that a single optical mode is interacting with a single mechanical mode, be it the center of mass oscillation. The mirror can thus be modeled as a mass attached to a spring of characteristic frequency and friction coefficient ; it is described by dimensionless position and momentum operators , which obey the canonical commutation relation . The optical mode has frequency and decay rate ; it is described by annihilation/creation operators and , which obey the canonical commutation relation .

In our analysis the cavity is assumed to be driven by an external laser which, to begin with, we take to have constant power and quasi-resonant frequency . In this context a periodic modulation is inserted at the level of the spring constant, which we express as the following time dependent parametric rescaling of the mirror frequency

 ω2(t)=ω2M[1+ϵcos(Ωt)], (1)

with . Accordingly the Hamiltonian of the system writes as (23)

 ^H= ℏωC^a†^a+ℏωM2^p2+ℏωM2[1+ϵcos(Ωt)]^q2 −ℏG0^a†^a^q+iℏE(e−iωLt^a†−eiωLt^a), (2)

where is the optomechanical coupling rate and is the driving rate. Including dissipation and decoherence effects the system dynamics can then be described with the following set of quantum Langevin equations (10)

 ⎧⎪⎨⎪⎩∂t^q=ωM^p,∂t^p=−ωM[1+ϵcos(Ωt)]^q−γM^p+G0^a†^a+^ξ,∂t^a=−(k+iΔ0)^a+iG0^a^q+E+√2k^ain, (3)

which we have written in a frame rotating at . Here is the unperturbed cavity laser detuning while is the radiation vacuum input noise with autocorrelation function (24)

 ⟨^ain(t)^a†in(t′)⟩=δ(t−t′). (4)

Similarly is the Brownian noise operator describing the dissipative friction forces acting on the mirror. Its autocorrelation function satisfies the relation (25)

 ⟨{^ξ(t),^ξ(t′)}⟩=2γMωM∫dω2πωcoth(ℏω2kBT)e−iω(t−t′). (5)

which for the specific case of an harmonic oscillator with a good quality factor , acquires the same Markov character of Eq. (4), i.e.

 ⟨{^ξ(t),^ξ(t′)}⟩≈2γMcoth(ℏωM2kBT)δ(t−t′), (6)

(this is a consequence of the fact that for only resonant noise components at frequency do sensibly affect the motion of the system). In the above expressions is the system temperature while is the anti-comutator (26).

### ii.1 Solving the dynamics

The evolution of the system is ruled by a set (3) of non-linear stochastic differential equations with periodic coefficients, whose solution is in general very difficult. In the following we will then introduce some useful approximations to simplify the calculations. First, we expand each operator as the sum of a c-number mean value and a fluctuation operator, i.e.

 ^a(t) = ⟨^a(t)⟩+(^a(t)−⟨^a(t)⟩)≡A(t)+δ^a(t), ^q(t) ≡ Q(t)+δ^q(t), ^p(t) ≡ P(t)+δ^p(t). (7)

We recall that the cavity is usually driven by a very strong laser in order to attain satisfactory levels of optomechanical interaction, so that the mean value will be much bigger than the fluctuations, which are due to the presence of random noise. This allows us to write (3) as two different sets of equations, one for the mean values (8), one for the fluctuations (9) and linearize the latter neglecting all terms which are second order small, obtaining

 ⎧⎪⎨⎪⎩∂tQ=ωMP,∂tP=−ωM[1+ϵcos(Ωt)]Q−γMP+G0|A|2,∂tA=−(k+iΔ0)A+iG0AQ+E, (8)
 ∂t⎛⎜ ⎜ ⎜ ⎜⎝δ^qδ^pδ^Xδ^Y⎞⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜⎝0ωM00−ωM[1+ϵcos(Ωt)]−γMG0Re[A]G0Im[A]−G0Im[A]0−kΔ0−G0QG0Re[A]0−Δ0+G0Q−k⎞⎟ ⎟ ⎟ ⎟⎠⋅⎛⎜ ⎜ ⎜ ⎜⎝cδ^qδ^pδ^Xδ^Y⎞⎟ ⎟ ⎟ ⎟⎠+⎛⎜ ⎜ ⎜ ⎜⎝0^ξ^Xin^Yin⎞⎟ ⎟ ⎟ ⎟⎠, (9)

where we have introduced the phase and amplitude quadratures for the cavity and the input noise fields, i.e. , , and . Equation (9) can be also expressed in a more compact form

 ∂t^u=S^u+^ζ, (10)

with being a time-dependent matrix, and with and being the column vectors of elements and , respectively. We stress that Eqs. (8) and (9) must be solved in the correct order, because the mean values , and play the role of coefficients in the equations for the fluctuations.

Equation (8) is nonlinear but can be solved numerically. Assuming that we are far from optomechanical instabilities and that we keep the modulation strength small enough to avoid additional instabilities due to parametric amplification, one finds that the mean values evolve toward an asymptotic periodic orbit with the same periodicity of the applied modulation. In this regime, an approximate analytic solution can also be derived, which we detail in Appendix A. Indeed since the modulation strength is not too strong, one can guess a perturbative expansion of the form

 Q(t)=∞∑j=0Q(j)(t), (11)

where does not depend on , is linear in , is quadratic in and so on. It turns out that each order is exactly solvable, as long as previous orders are known. This originates a chained set of equations and by keeping a finite number of orders , we can finally obtain the asymptotic solution up to the desired precision (e.g. see Fig. 2).

Equation (9) is stochastic and needs some more manipulation. Nonetheless since we have linearized the dynamics and the noises are zero-mean gaussian noises, fluctuations in the stable regime will also evolve to an asymptotic zero-mean gaussian state. The state of the system is then completely described by the correlation matrix of elements

 Cij(t)=Cji(t)=12⟨^ui(t)^uj(t)+^uj(t)^ui(t)⟩, (12)

whose evolution can be derived directly from equations (10) and (12):

 ∂tC=SC+CS⊤+N, (13)

where is transpose of , and where is the diagonal noise correlation matrix with diagonal entries , defined by

 12⟨^ζi(t)^ζj(t′)+^ζj(t′)^ζi(t)⟩≡Nijδ(t−t′). (14)

Equation (13) is now an ordinary linear differential equation. We know that its solution evolves toward a unique asymptotic configuration (independently of the initial state), proven that the eigenvalues of the matrix have negative real part for all times , which can be verified by applying the Routh-Hurwitz criterion (27). Again, we can either solve Eq. (13) numerically or obtain an approximate analytic solution with a perturbative expansion in (see Appendix A for the latter).

### ii.2 Quantum properties of the system

As already mentioned, thanks to gaussianity of the asymptotic solution all relevant informations about the system can be extracted directly from the correlation matrix . In particular we will focus on the following quantities: the number of phonons in the mirror, the squeezing in the mirror and in the radiation quadratures, and the nonclassical correlation between the mirror and the radiation degrees of freedom.

The number of phonons can be expressed using the approximate relation

 ℏωM(n+12) ≈ ℏωM/2⟨δq2+δp2⟩ (15) = ℏωM/2(C11+C12),

which holds if the modulation of the mechanical frequency is not too strong. This tells how far the system is from the ground state. Since both and are periodic in time, we will identify the number of phonons with the maximum over one period of the modulation, i.e.

 nMAX=maxτ{n(t)}, (16)

(here and in the following represents an optimization with respect to a time interval with being sufficiently larger than to guarantee that the system has reached the asymptotic steady state).

Squeezing of the generalized mirror quadratures is also easily found:

 ⟨δq2θ⟩=C11cos2θ+C22sin2θ+(C12+C21)cosθsinθ. (17)

Again we construct a time independent quantity to deal with. First, for each time we select the parameter for which is minimum. In terms of covariance matrix, this is just the smaller eigenvalue of the block matrix . We then minimize this quantity with respect to time over a period . This tells how much squeezing can be produced at most.

 Δ2qMIN=minτ{minθ⟨δq2θ⟩}. (18)

 Δ2XMIN=minτ{minθ⟨δX2θ⟩}. (19)

Non-classical correlations in the system can be described using quantum discord  (21), which includes entanglement as well as more general quantum correlations that are shown also by separable states (28). For a gaussian state, is easily constructed from the correlation matrix as demonstrated in (29). Time dependance is then eliminated by considering

 DMAX=maxτ{D(ρ(t))}. (20)

Entanglement alone will be specifically described using logarithmic negativity (30), which is also easily constructed from the correlation matrix as demonstrated in (31). Again, time dependance is eliminated by considering

 ENMAX=maxτ{EN(ρ(t))}. (21)

## Iii Results

We now present the results obtained by solving the dynamics of the system as detailed in the previous section. The parameters used in our analysis are ng, MHz, Hz, K, , mm, MHz, nm and mW: this choice is compatible with values attained in state of the art experiments and is also consistent with the stability requirement of section II.1 (furthermore, under the condition the optical and the mechanical variables are brought at resonance). The strength and the frequency of the modulation are left as variable parameters instead, since we want to characterize the optimal modulation regime, e.g. which and maximize the visibility of quantum effects.

In Fig. 2, we temporarily fix , (this particular choice will be justified in the following) and we report the solution of Eq (8) for the mean values and of the mirror position and momentum. We see that the evolution tends indeed to an asymptotic periodic orbit, which is very well approximated by the analytic solution.

We then focus on the solution of equation (13) and we plot the quantities described in section II.2, for multiple values of , . In particular: Fig. 3 shows the maximum number of phonons in the mirror, computed via Eq. (16); the maximum of the logarithmic negativity, computed via Eq. (21); the maximum of the quantum discord, computed via Eq. (20); and the minimum variance of all the mirror generalized quadratures, computed via Eq. (18).

As evident from the plots, the level of squeezing and entanglement is maximum when the modulation frequency is and increases monotonically with respect to the strength , until the system eventually reaches an instability point for too strong modulations (in the above figures, this instability is represented by a blank region around the point , ). It is also clear that the optimal modulation, the one that most enhances quantum effects, is also responsible for heating the system far from its ground state.

We can understand this behavior if we interpret Eq. (13) as describing the dynamics of a set of (classical) parametric oscillators with canonical coordinates defined by the correlations functions (12), which evolve under the action of damping and constant external driving forces. Indeed, by a close inspection of the matrix one notices that such oscillators possess natural frequencies which are periodically modulated through functions (i.e. , , and the direct term ) that, in first approximation, evolve sinusoidally with the same frequency – see Eq. (41) in Appendix A for details. Moreover, in the stability region we are sure that parametric modulation pumps energy into the system at a lower rate with respect to losses, since the system evolves toward a stationary orbit: we call this regime “below-threshold” to distinguish it from the exponential amplification usually associated with parametric oscillators. For this model phase-locking is expected to occur when matches the zero-order eigenfrequencies defined by the constant part of (and not twice this frequencies as in the case of parametric instability), resulting in an enhancement of the oscillations of the effective coordinates (12) and hence of the associated quantum effects defined in Sec. II.2 (32) (more details are found in Appendix A). It turns out that, at least for the figure of merit we are concerned here (i.e. , , , etc) the relevant frequency is indeed .

To see this, we can procede by steps. First of all notice that from the numerical solution, we can guarantee that the system is not unstable (see Fig. 3), i.e. that it is indeed in the below-threshold regime. Next, consider the case of no coupling () and no modulation (): we stress out three relevant aspects. First, the mechanical part and the radiation part are independent, so there is no entanglement. Second, each subsystem evolves with the Hamiltonian of a quantum harmonic oscillator, so the quadrature mean value evolves with a phase and the variance with a phase (we remind that we fixed ). This tells us that, at least in this regime, the frequencies which govern the quantities of interest are degenerate at the value . Third, each subsystem is also coupled to its own environment and will eventually relax to a thermal state characterized by , so there is no squeezing. Now turn on the coupling : this has three main effects. First it introduces entanglement in the system (11) (). Second the eigenfrequencies are brought out of degeneracy and shifted by a term (33), which is quite small with respect to for our choice of values (confirming that indeed the latter is the resonant value at which the modulation should provide an enhancement). Third, backaction cooling (5) is now active and the oscillator approaches the ground state . Squeezing is still absent at this level. Finally, turn on the modulation (1). Thanks to the phase-locking mechanism we have anticipated previously and detailed in Appendix A this will yield an enhancement of the correlations when matches the natural frequency . For instance, for the negative entropy and for mirror variance , we get

 EN ∼ E0+ϵK1(Ω)cos(Ωt+φ1), ⟨δq2θ⟩0 ∼ 1/2+ϵK2(Ω)cos(Ωt+φ2), (22)

where and are associated response functions analogous to the Lorentzian response of a simple harmonic oscillator (though an exact expression is rather cumbersome in our specific case) and are peaked around . We see that the quadrature gets periodically squeezed over time and entanglement is periodically increased to higher values with respect to the unmodulated case. In addition, these effects increase monotonically with , up to the instability threshold. A similar enhancement of the entanglement is also described in Ref. (16), where two harmonic oscillators are coupled via linear interaction and the coupling constant is a periodic function of time. This time dependance produces an effective modulation on the normal frequencies of the system: as a result, entanglement is shown to increase and become much more robust against temperature. This agrees very well with what we found here.

## Iv Interplay between two different modulations

Results analogous to those presented in the previous section has been found very recently by Mari and Eisert (14); (15), for an optomechanical system driven with an amplitude modulated input laser. For clarity, we rewrite their Hamiltonian

 ^H= ℏωC^a†^a+ℏωM2(^p2+^q2)−ℏG0^a†^a^q +iℏ(E+E1cos(Ωt))(e−iωLt^a†−eiωLt^a). (23)

At first sight, the situation appears to be somewhat different from our initial problem. In Eq. (23), internal parameters of the system are left unchanged; it is instead the external driving that undergoes an oscillatory behavior. Nevertheless the effects are strikingly similar: high levels of squeezing can be attained when the frequency of modulation is (14), and the same regime is also optimal to enhance entanglement between mechanical and radiation modes (15) The authors themselves comment that “…this dynamics reminds of the effect of parametric amplification, as if the spring constant of the mechanical motion was varied in time with just twice the frequency of the mechanical motion, leading to the squeezing of the mechanical mode…” (14).

In fact there is a strong analogy between the two cases. Independently of what Hamiltonian (2) or (23) one chooses, far from instability regions the mean values , , will be characterized by an asymptotic periodic orbit with the same periodicity of the applied modulation . This assures that in both cases the equation (13) for the covariance matrix has the same linear form, with being a periodic function of time (in the limit ) and being a constant driving. The conclusions we derived in section III, must therefore hold, at least qualitatively, also for the system studied in (14); (15).

An interesting question now rises. What if the two modulations are applied together? Can they interfere, either constructively or destructively, and sensibly alter the one-modulation picture?

To get an answer, we consider a new composite system, described by the Hamiltonian

 ^H = ℏωC^a†^a+ℏωM2^p2+ℏωM2(1+ϵcos(Ω1t))^q2 (24) −ℏG0^a†^a^q +iℏE(1+ηcos(Ω2t+ϕ))(e−iωLt^a†−eiωLt^a).

Note that we explicitly introduced a relative phase between the two applied modulation: if we expect any interference, the properties of the system should indeed depend on this new variable.

The analysis presented in the previous sections is straightforwardly generalized to the present case, so we will skip directly to the results [details can be found however in the Appendix]. Taking the same parameters as in Sec. III, we choose the optimal modulation frequencies and fix , (this is the same value used in (14)). These modulation strengths give comparable squeezing performances when considered singularly, and also assure that we are reasonably far from the instability region. To present the results, we plot the quantities introduced in Sec. II.2 against the relative phase in Fig. 4. An interference pattern is indeed evident and each of the above quantities oscillates between a minimum and a maximum as varies in the range . However, entanglement and quantum discord are affected very weakly and do not differ much from our initial one-modulation case. Besides we see that in order to generate quantum correlations, a modulation of the mechanical frequency is more suitable than a modulation of the driving laser amplitude. Adding the second modulation to the first is of little effect.

Squeezing generation instead, presents very interesting features. First, as we said, we choose two modulations that give comparable levels of squeezing when applied individually. Moreover, when applied together, they can strongly interfere. For example we see in Fig. 4(d) that for a phase , rises toward the threshold value and squeezing becomes weaker. Each modulation taken alone would generate more squeezing than the two combined: this is an unambiguous sign of a disadvantageous interplay. For a phase we find instead a great advantage in applying two modulations: is lowered to a value , a considerable performance if compared to our initial one-modulation case where instabilities prevent us from reaching . In fact, not only we attain the same high levels of squeezing, but we are also well inside the stability region, so that we could increase both and to perform even better.

We also note that the optimal(worst) phase choice for squeezing generation also corresponds to maximum heating(cooling) of the mirror, as can be seen in Fig. 4(a). This is another confirmation that parametric oscillation is indeed the main underlying mechanism: in fact not only does a stronger modulation enhance the generation of quantum effects, as inferred from equations (22), it also pumps more energy into the system.

We then see how the interplay between two independent modulations can be carefully exploited to increase levels of squeezing in an optomechanical system.

## V Conclusions

We have studied in great detail the effect of periodic modulations on optomechanical systems and we have characterized several ways in which such modulations can be exploited to enhance relevant quantum properties including squeezing, entanglement and quantum discord. While the idea that modulations can help accessing the quantum regime was already known from previous works (14); (15); (16), we have proposed a new interpretation of this enhancement mechanism in terms of a resonance between the modulation frequency and the natural frequencies of the system. This simple model allowed us to prove the existence of an optimal modulation regime and to understand the arising of instability thresholds. Finally, we have analyzed the interplay of different modulations and we have found that constructive (destructive) interference effects may arise when they are applied simultaneously, causing a further enhancement (a suppression) of quantum effects. We believe that these results could lead further on toward the development of optimal control strategies.

###### Acknowledgements.
We thank Rosario Fazio for useful discussions and comments. This work was supported by MIUR through FIRB-IDEAS Project No. RBID08B3FM.

## Appendix A Asymptotic behaviors

This section deals with some technical aspects related with the asymptotic solutions of Eqs. (8) and (13), which define the quantum properties of the system. Here we discuss the resonant mechanisms which is responsible for the enhancement of quantum effects at , as well as the role of the relative phase in the interplay between different modulations. We start in Sec. A.1 by presenting a simple paradigmatic case which captures the main aspects of the resonance. Then, in section A.2, we introduce the analytic framework which will be used to describe the dynamics of the system. Finally, sections A.3 and A.4 are devoted to analyze in details the asymptotic behavior of the system, in the one- and two-modulation scenario respectively.

### a.1 Single oscillator model

As anticipated in the main text the evolution of the correlations matrix describes the dynamics of a multi-dimensional (classical) oscillator which evolves in presence of damping and external constant driving (defined by the matrix ) and which possesses characteristic frequencies (determined by ) that are externally modulated at frequency [these statements are explicitly verified in Secs. A.2A.3 and A.4]. To enlighten the role of the modulation in the evolution of the correlation functions it is hence worth focusing on the simplest example of this sort. This is provided by a single parametric oscillator whose position evolves according to the equation

 ¨x(t)=−ω20[1+αcos(νt)]x(t)−γ˙x(t)+F, (25)

with and being the characteristic and the modulation frequency, being the amplitude of the modulation, being the damping rate and the strength of a constant driving. For this simple scenario, two cases are possible. If (above-threshold condition), parametric modulation pumps energy into the system at a faster rate with respect to dissipation; the system increases its energy exponentially and is therefore unstable. If (below-threshold condition), the system reaches a stationary regime, given by the balance of pumping and dissipation. We can then look for a stable solution of Eq. (25), assuming that is small and treating the solution perturbatively, i.e. . To order zero in the system is just a damped driven harmonic oscillator, which relaxes toward its equilibrium position . To first order in , the long time solution is then given by

 ¨x(1)(t) =−ω20x(1)(t)−ω20cos(νt)¯x(0)−γ˙x(1)(t). (26)

Therefore we see that the parametric modulation, for the below-threshold regime, can be mapped onto an effective external driving and the solution is easily found to be

 x(t)≃Fω20+αf(ν)Fcos(νt+ϕ), (27)

with being the Lorentzian response function of a classical harmonic oscillator. Clearly the superimposed oscillation, which we remind is an effect of the parametric modulation, will be much greater near resonance with the natural frequency and for just below the instability threshold. Going to second order in yields small deviation from this picture and we can stop our qualitative analysis here. In summary, parametric modulation can controllably enhance oscillations of the system coordinates if two main conditions are satisfied: the modulation must not be too strong, otherwise the system becomes unstable, and an external (constant) driving must also be applied, otherwise the system relaxes to (as from Eq. 27 with ). We also stress out that, in the below-threshold regime, the resonance condition is given by (i.e. the modulation frequency should be the same as the natural frequency of the system) and not by , as is the usual case of exponential parametric amplification.

### a.2 General treatment of the modulated optomechanical system

Turning back to Eq. (13), we will see that all conditions are indeed satisfied: the coefficient matrix is periodically modulated over time, stability can be verified with a numeric solution and external driving is provided by the noise correlation function . The above result implies that in the case of our multi-dimensional parametric oscillator, maximum enhancement of the oscillations is expected when matches the characteristic frequencies that govern the dynamics of the correlation functions in absence of the modulation. The latter are defined by the matrix of Eq. (13) when (and in the two-modulation scenario). As mentioned in the text, at least for the figure of merit we are concerned about in the paper (i.e. , , , etc), the relevant frequency is indeed .

We can thus generalize the simple model of Sec. A.1 to the present case and reproduce the numerical results we found in the main text with a semi-analytic solution of Eqs. (8) and (13), which we briefly sketch here. In doing so, we will also identify and comment the relevant points which are responsible for the behaviours observed in Figs. 3 and 4.

#### Classical solution.

Let us start with Eq. (8) for the mean values which, for the sake of completeness, we report here for the general scenario defined by the Hamiltonian (24) where both the frequency modulation (2) and the amplitude modulation (23) are activated, i.e.

 ⎧⎪⎨⎪⎩∂tQ=ωMP,∂tP=−ωM[1+ϵcos(Ωt)]Q−γMP+G0|A|2,∂tA=−(k+iΔ0)A+iG0AQ+E[1+ηcos(Ωt+ϕ)], (28)

having only assumed their frequencies to be identical, i.e. . As anticipated in the text – see Eq. (11) – we look for a perturbative solution in the modulations strengths and , i.e

 ⎧⎪⎨⎪⎩Q=Q(0)+Q(1)+Q(2)+…,P=P(0)+P(1)+P(2)+…,A=A(0)+A(1)+A(2)+…, (29)

where, for instance, is linear in and , is quadratic in and and so on (note that we can revert to the single modulation scenario simply by imposing ). At order zero we get

 ⎧⎪⎨⎪⎩∂tQ(0)=ωMP(0),∂tP(0)=−ωMQ(0)−γMP(0)+G0|A(0)|2,∂tA(0)=−(k+iΔ0)A(0)+iG0A(0)Q(0)+E. (30)

From the numeric simulation we know that this non-linear equation evolves toward a stable point (, , ) and by setting the derivatives to zero, we can find these asymptotic values. Next, at first order we get

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂tQ(1)=ωMP(1),∂tP(1)=−ωMQ(1)−ωMϵcos(Ωt)¯Q(0)−γMP(1)+G0(¯A(0))∗A(1)+G0(A(1))∗¯A(0),∂tA(1)=−(k+iΔ0)A(1)+iG0¯A(0)Q(1)+iG0A(1)¯Q(0)+Eηcos(Ωt+ϕ). (31)

These are the equations of three coupled and forced harmonic oscillators, with forcing terms and that are purely oscillating. In addition, damping makes sure that the system is stable. The solutions are easily obtained in the form

 ⎧⎪⎨⎪⎩Q(1)(t)=q1eiΩt+q∗1e−iΩt,P(1)(t)=p1eiΩt+p∗1e−iΩt,A(1)(t)=a1eiΩt+a2e−iΩt, (32)

with , , , and being complex parameters which can be computed by replacing (32) into Eq. (31). Hence, to first order in and , the effect of the modulation on classical values is to add an oscillating term with frequency and mean value . The amplitude of this oscillation clearly depends on the forcing term, i.e on the amplitudes , , on their relative phase and on the frequency (via the oscillator response function). Finally, the equations for second order are

 ⎧⎪⎨⎪⎩∂tQ(2)=ωMP(2),∂tP(2)=−ωMQ(2)−ωMϵcos(Ωt)Q(1)−γMP(2)+G0A∗(0)A(2)+G0A∗(1)A(1)+G0A∗(2)A(0),∂tA(2)=−(k+iΔ0)A(2)+iG0A(2)Q(0)+iG0A(1)Q(1)+iG0A(0)Q(2). (33)

These are again the equations of three coupled and forced harmonic oscillators, but this time the forcing terms , and have also a constant part. As for the first order corrections the solutions are easily obtained in the form

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Q(2)(t)=¯Q(2)+q3ei2Ωt+q∗3e−i2Ωt,P(2)(t)=¯P(2)+p3ei2Ωt+p∗3e−i2Ωt,A(2)(t)=¯A(2)+a3ei2Ωt+a4e−i2Ωt. (34)

To second order in and , the effect of the modulation on classical values is thus to add a constant shift and an additional oscillating term with frequency and mean value . Higher orders can be processed in the same way but for the parameter region we have selected in the main text, one can limit the analysis to second order since already at this point we get the correct result within a good degree of accuracy (see Fig 5). Full convergence of the approximation when higher orders are included can be seen from Fig 2 in the main text, where we plot the numeric evolution of classical values and and the analytic counterpart, computed up to order six.

#### Linearized quantum solution.

We now turn to Eq. (9) for the quantum fluctuation, which we rewrite below

 ∂tC=SC+CS⊤+N.

Recall that the matrix depends on and via the classical values and an additional explicit term . If we make use of the approximate solution found before, we can thus identify a matrix independent of the perturbation, a matrix linear in and and a matrix quadratic in and . Again, we look for a perturbative solution for the matrix , i.e

 C=C(0)+C(1)+C(2)+…. (35)

where is linear in and , is quadratic in and and so on. The calculations simply follow what we have done for the classical part. At order zero we get

 ∂tC(0)=S(0)⋅C(0)+C(0)⋅S(0)⊤+N. (36)

This equation is linear and evolves toward a stable point , which we can find by setting the derivatives to zero. Next, at first order we get

 ∂tC(1)=S(0)⋅C(1)+C(1)⋅S(0)⊤+S(1)⋅¯C(0)+¯C(0)⋅S(1)⊤. (37)

These are the equations of sixteen coupled and forced harmonic oscillators, with forcing terms that are purely oscillating. Since is real, the solutions are easily obtained in the form

 C(1)(t)=c1eiΩt+c∗1e−iΩt. (38)

Hence, to first order in and , the effect of the modulation on the correlations is to add an oscillating term with frequency and mean value . As in the case of an unidimensional resonator, the amplitude of this oscillation will be greater when the modulation frequency is chosen in resonance with the eigenfrequencies of the normal modes. Finally, the equations for second order are

 ∂tC(2)=S(0)⋅C(2)+C(2)⋅S(0)⊤+S(1)⋅C(1)+C(1)⋅S(1)⊤+S(2)⋅¯C(0)+¯C(0)⋅S(2)⊤. (39)

These are again the equations of sixteen coupled and forced harmonic oscillators, but this time the forcing terms and have also a constant part. The solutions are easily obtained in the form

 C(2)(t)=¯C(2)+c3ei2Ωt+c∗3e−i2Ωt. (40)

Hence, as for the linear solutions, to second order in and , the effect of the modulation on the correlations is to add a constant shift and an additional oscillating term with frequency and mean value .

### a.3 Asymptotic behavior in the single modulation regime

#### Classical solution.

We can come back to the single modulation scenario by putting in the above analysis. By doing so, we can get approximate analytic expressions for the asymptotic mean values , and . However the complete formulas are too long to be reported here and we must limit ourselves to a semi-numeric expression, where we substitute all values as in Sec. III except for the interesting parameter . For example we report the expression of the mirror position

 Q(t) =14684.7−ϵ22784.43 Missing or unrecognized delimiter for \Big +ϵ2(164.97cos(2Ωt)−0.50sin(2Ωt)). (41)

As anticipated in the main text, to first order in the mean values have an asymptotic oscillatory behavior, which well describes the exact asymptotic solution. To be precise however, we cannot neglect the second order contributions: indeed, while second harmonic oscillations are one order of magnitude smaller, the constant shift is comparable to first order effects and must be taken in account.

#### Linearized quantum solution.

We can also look at the quantum properties of the system in the asymptotic regime, as a function of the modulation strength . Fixing all other parameters to values in the text, we find for example the following expression for the number of phonons in the mirror

 nphon(t) =0.08+ϵ24.14 +ϵ(0.14cos(Ωt)−0.01sin(Ωt)) −ϵ2(0.02cos(2Ωt)−0.21sin(2Ωt)). (42)

For completeness we also report the expression for the single correlation

 C11(t) =0.56+ϵ24.01 +ϵ(0.28cos(Ωt)−1.63sin(Ωt)) +ϵ2(0.03cos(2Ωt)−0.20sin(2Ωt)). (43)

We see that (and similarly ) has strong oscillations in time proportional to . Hence we can say, at least qualitatively, that squeezing will be dominated by first order effects in the range of values considered. On the contrary the number of phonons can be considered time-independent, with oscillations that are negligible if compared to the constant term. Moreover, we know that the mirror is cooled close to its ground state when the system is unmodulated. Therefore the number of phonons is strongly dependent on , and the constant shift due to second order effects becomes quickly the dominant effect. These results agree very well with the numerical simulation summarized in Fig. 3.

We conclude this section with one last comment on why the number of phonons is constant in time. We know that the position and momentum of the mirror oscillate with frequency and a relative phase shift of (slight modifications being induced by the interaction with the optical subsystem). In the same way and oscillate with twice this frequency, i.e. , and with twice this relative phase shift , i.e. . In turn, the two oscillations cancel each other out when summing and , thus giving a time-independent number of phonons. In addition we see that the modulation is most effective on the mirror correlations when , as stated before.

### a.4 Asymptotic behavior in the two modulation regime

#### Classical solution.

We now reintroduce the second modulation and study the interplay between the two. Again we would like to fix all parameters except , and to the values found in the main text. However, already at the classical level, expressions for , and tend to become rather long and complex since we have now 3 free parameters. Therefore we will substitute also the numerical values of and (values are found in Sec. IV). This is not so bad: indeed recall that we are particularly interested in the dependence of quantum properties on the relative phase . For example we report the expression of the mirror position

 Q(t) =17523.4−357.13cos(ϕ)+315.98sin(ϕ) +1484.13cos(Ωt)−4.43sin(Ωt) +2201.22cos(Ωt+ϕ)−1870.83sin(Ωt+ϕ) +14.84cos(2Ωt)−0.04sin(2Ωt) +22.62cos(2Ωt+ϕ)−18.78sin(2Ωt+ϕ) +113.53cos(2Ωt+2ϕ)−32.43sin(2Ωt+2ϕ). (44)

Without losing much time on the cumbersome formula above, we only point out that again first order effects are dominating (as we can see by comparing oscillations at and oscillations at ). However we also see that, already at the classical level, the the phase has a strong influence on the amplitude of oscillations. This is a clear sign that the phase plays indeed an important role in the system dynamics.

#### Linearized quantum solution

We turn now to quantum properties of the system. Again we fix all parameters at the values used above, except for the relative phase . To second order in the perturbation, we get expressions for the number of phonons and the correlation that depends on the phase as

 nphon(t)=12( 1.167+0.087cos(ϕ)+0.753sin(ϕ)+0.086cos(Ωt)−0.005sin(Ωt) −0.004cos(2Ωt)+0.037sin(2Ωt)+0.008cos(Ωt+ϕ)−0.006sin(Ωt+ϕ) −0.024cos(2Ωt+ϕ)+0.005sin(2Ωt+ϕ)+O(10−5)), (45)
 C11(t)= 1.09+0.02cos(ϕ)+0.39sin(ϕ)+0.08cos(Ωt)−0.48sin(Ωt) +0.002cos(2Ωt)−0.02sin(2Ωt)+0.39cos(Ωt+ϕ)−0.05sin(Ωt+ϕ) +0.001cos(2Ωt+ϕ)−0.002sin(2Ωt+ϕ)+O(10