Time-evolution of nonlinear optomechanical systems: interplay of mechanical squeezing and non-Gaussianity

We solve the time evolution of a nonlinear optomechanical Hamiltonian with arbitrary time-dependent mechanical displacement, mechanical single-mode squeezing and a time-dependent optomechanical coupling up to the solution of two second-order differential equations. The solution is based on identifying a minimal and finite Lie algebra that generates the time-evolution of the system. This reduces the problem to considering a finite set of coupled ordinary differential equations of real functions. To demonstrate the applicability of our method, we compute the degree of non-Gaussianity of the time-evolved state of the system by means of a measure based on the relative entropy of the non-Gaussian state and its closest Gaussian reference state. We find that the addition of a constant mechanical squeezing term to the standard optomechanical Hamiltonian generally decreases the overall non-Gaussian character of the state. For sinusoidally modulated squeezing, the two second-order differential equations mentioned above take the form of the Mathieu equation. We derive perturbative solutions for a small squeezing amplitude at parametric resonance and show that they correspond to the rotating-wave approximation at times larger than the scale set by the mechanical frequency. We find that the non-Gaussianity of the state increases with both time and the squeezing parameter in this specific regime.

entropy of the non-Gaussian state and its closest Gaussian reference state. We find that the addition of a constant mechanical squeezing term to the standard optomechanical Hamiltonian generally decreases the overall non-Gaussian character of the state. For sinusoidally modulated squeezing, the two secondorder differential equations mentioned above take the form of the Mathieu equation. We derive perturbative solutions for a small squeezing amplitude at parametric resonance and show that they correspond to the rotating-wave approximation at times larger than the scale set by the mechanical frequency. We find that the non-Gaussianity of the state increases with both time and the squeezing parameter in this specific regime.
Keywords: optomechanics, nonlinear optomechanics, quantum optomechanics, non-Gaussianity, squeezing (Some figures may appear in colour only in the online journal)

Introduction
The mathematical understanding of optomechanical systems operating in the nonlinear quant um regime is a major topic of current interest. While most experiments effectively undergo linear dynamics, governed by quadratic Hamiltonians that emerge following a 'linearisation' procedure [1][2][3], many experiments now operate in the fully nonlinear regime [4][5][6] where this procedure fails. It is therefore highly desirable to provide a complete and analytic characterisation of the fully nonlinear system dynamics. Analytic solutions have previously been found for a constant light-matter coupling [7,8] and, more recently, the timedependent case was solved [9].
The inherently nonlinear interaction between the optical field and the mechanical element in an optomechanical system allows for the generation of non-Gaussian states. Starting from a broad class of initial states, including coherent states, the vacuum, and thermal states, this is only possible in the nonlinear regime; in contrast, quadratic Hamiltonians take input Gaussian states to output Gaussian states. As such, investigating the non-Gaussianity of optomechanical states can only be performed once the time-evolution in the nonlinear regime has been solved, which is the primary aim of this work. Interestingly, a number of non-classical and non-Gaussian states have been found to constitute an important resource for sensing. Schrödinger cat states [7,8], compass states [10,11] and hypercube states [12]-which are all non-Gaussian and highly non-classical states-have all been found to have applications for sensing. More generally, the detection and generation of non-Gaussianity in optomechanical systems has been extensively studied in theoretical proposals [13][14][15] as well as in experiments [4,5,16]. Beyond optomechanics, the presence of a nonlinear element is also key to a number of quantum information tasks, such as obtaining a universal gate set for quantum computing [17,18], teleportation [19], distillation of entanglement [20,21], error correction [22], and non-Gaussianity has been explored as the basis of an operational resource theory [23][24][25].
Optomechanical systems offer a natural nonlinear coupling which, if strong enough, may lead to substantial non-Gaussianity in the evolved state. It is therefore essential to better understand the dynamics of such systems, with special emphasis on the interplay between nonlinearities and other Hamiltonian terms in this dynamics. An important question to be answered is thus how do the different aspects of an optomechanical system affect the non-Gaussianity of the state at a given time? A preliminary study of non-Gaussianity in standard optomechanical systems provided the first tools to approach this question [9], however, optomechanical systems can exhibit additional, potentially more interesting, effects. An important non-classical effect that can be included into optomechanical systems is squeezing of the optical or mechanical modes. The addition of squeezing has been found to be beneficial for sensing since it increases the sensitivity in a specific field quadrature. For example, it has been shown that squeezed light injected into LIGO can be used to enhance the detection of gravitational waves [26]. Similarly, mechanical squeezing can aid the amplification and measurement of weak mechanical signals [27].
In this work, we study the non-Gaussianity of a quantum system of two bosonic modes characterised by an optomechanical Hamiltonian with the standard cubic light-matter interaction term, and with the addition of a mechanical displacement term and a mechanical squeezing term. We extend a recently developed solution of the time evolution operator induced by a plain optomechanical Hamiltonian [9] to include the additional terms of interest here. Interestingly, for time-dependent squeezing modulated at resonance, we find that the dynamics are governed by the well-studied Mathieu equation. We subsequently derive perturbaive solutions and show that these coincide with the physically intuitive rotating-wave approximation (RWA) for large times. The decoupling methods used in this work have a long tradition in quantum theory [28][29][30] and were recently applied to problems such as the one at hand [9,31,32]. We use the resulting analytic solutions to compute the amount of non-Gaussianity of the state using a measure of relative entropy [33,34] for both a constant and a time-dependent mechanical squeezing parameter.
Our results indicate that the non-Gaussian character of an initially coherent state decreases in general with an increasing squeezing parameter. However, when the squeezing is applied periodically at twice the mechanical resonance, the non-Gaussianity increases approximately linearly with time and the amplitude of the squeezing. The competition between the amount of squeezing and the strength of the nonlinear term is difficult to compute explicitly; instead, we provide asymptotic expressions in terms of upper and lower bounds to the non-Gaussianity in different regimes. A conclusive answer requires further investigation, potentially providing a concise expression where such competition can be easily understood.
The paper is structured as follows. In section 2, we introduce the nonlinear Hamiltonian with mechanical squeezing. This is followed by section 3 where we provide a short introduction to the methods used to solve the dynamics. The full derivation can be found in appendix B. Following this, we review the measure of non-Gaussianity and derive expressions for an asymptotic expression and a reduced measure in section 4. In section 5, we then specialise to two specific cases and compute the amount of non-Guassianity for constant squeezing (section 5.2), and for modulated squeezing (section 5.3). Finally, we conclude with a discussion in section 6 and some final remarks in section 7.

Dynamics
In this section we present the optomechanical Hamiltonian of interest to this work and explain the origin of the various terms. An extensive introduction to optomechanics can be found in the literature [1].

Hamiltonian
In this work we consider the two-mode Hamiltonian where Ĥ 0 := ω câ †â + ω mb †b is the free Hamiltonian, while ω c and ω m are the frequencies of the cavity mode and the mechanical resonator respectively.
The Hamiltonian (1) describes the dynamics of a number of different systems. For example, G(t) appears in optomechanical systems as a standard coupling term due to radiation pressure obtained for Fabry-Pérot cavities, where one end of the cavity is a mirror that can move freely [35]. Such coupling appears also within systems with a central translucent membrane in a rigid optical cavity [36], levitated nanodiamonds [37] or optomechanical crystals [38,39]. A depiction of a levitated nanosphere interacting with cavity modes can be found in figure 1.
The Hamiltonian (1) reduces to the standard optomechanical Hamiltonian when D 1 = D 2 = 0 9 . The term weighted by the coupling D 1 corresponds to an externally imposed displacement of the mechanical part, which can be induced by a piezoelectric element connected to its support or by an external acceleration, such as that caused by the gravitational force acting on the mechanical element [40,41]. The term governed by D 2 can be thought of as a modulation of the trap frequency and leads to squeezing of the mechanics, which can be externally imposed employing another strong optical field or an electrostatic force [42].

Dimensionless dynamics
To understand which dimensionless parameters are relevant to the dynamics of the system, we start by introducing dimensionless quantities and rescaling the Hamiltonian. Such a rescaling also serves to simplify the notation and any graphical representation of the system dynamics. We achieve this by dividing the functions in the Hamiltonian by the mechanical frequency ω m . The action corresponds to switching from the laboratory time t to τ = ω m t, where τ is the new, dimensionless time. The optical frequency becomes Ω c := ω c /ω m . In addition, the couplings in the Hamiltonian become We also rescale the Hamiltonian by , meaning that Ĥ becomeŝ which is the Hamiltonian that we will be working with. A levitated nanosphere in a cavity. The optical field is described by annihilation and creation operators â and â † , while the mechanics-in this case the mechanical motion of the nanosphere-is described by annihilation and creation operators b and b † . The system evolves under the optomechanical Hamiltonian (2). 9 We note here that non-Gaussianity for the case D 1 = D 2 = 0 has been already studied [9].

Solving the dynamics
Our aim is to provide the techniques to be used to understand the interplay of mechanical squeezing and non-Gaussianity in an optomechanical system, for which an analytic expression for the state evolution [9] is central. In this section we introduce the tools needed to solve the dynamics generated by (2). See appendix A for a more in-depth introduction to the underlying concepts, and appendix B for the full calculations.

Continuous variables and covariance-matrix formalism
When solving the dynamics, we employ methods from the continuous variable formalism [3,43]. Specifically, the methods are used to solve the time-evolution of the quadratic part of the system, and to describe its action on the nonlinear light-matter interaction term. We briefly review the continuous variable formalism here.
In recent years, thanks to progress in the mathematical framework provided by the covariance matrix formalism [3,43], it has become clear that Gaussian states constitute an extremely valuable toolkit to investigate quantum information processing in quantum setups, and in relativistic ones as well [44]. The main advantage is that the covariance matrix formalism provides a powerful set of mathematical tools to treat Gaussian states of bosonic fields that undergo linear transformations of the creation and annihilation operators fully analytically [43]. Ultimately, Gaussian states are the paramount resource for continuous variables quant um information processing and computation [17] and have become a standard feature in most quantum optics laboratories. However, it should also be pointed out that these methods can be used to describe the evolution of operators in the Heisenberg picture, even when the states considered are not Gaussian.
In quantum mechanics, the initial state ρ i of a system of N bosonic modes with operators {â n ,â † n } evolves to a final state ρ f through the standard Schrödinger equation ρ f =Ûρ iÛ † , where Û implements the transformation of interest, such as time evolution. If the state ρ is Gaussian and the Hamiltonian Ĥ is quadratic in the operators, it is convenient to introduce the vector X = (â 1 , . . . ,â N ,â † 1 , . . . ,â † N ) T , where T denotes the transpose of the vector. Similarly, the vector of first moments d := X and the covariance matrix σ are defined by σ nm := {X n ,X † m } − 2 X n X † m , where {·, ·} stands for anticommutator and all expectation values of an operator Â are defined by Â := Tr(Âρ).
In this language, the canonical commutation relations read [X n ,X † m ] = i Ω nm , where the 2N × 2N matrix Ω is known as the symplectic form [43]. We then notice that, while arbitrary states of bosonic modes are, in general, characterised by infinite real parameters, a Gaussian state is uniquely determined by its first and second moments, d n and σ nm respectively [43]. Furthermore, unitary transformations quadratic in the annihilation and creation operators, such as Bogoliubov transformations [45], preserve the Gaussian character of a Gaussian state and can always be represented by a 2N × 2N symplectic matrix S that preserves the symplectic form, i.e. S † Ω S = S Ω S † = Ω 10 . In a similar manner, the symplectic matrix S that encodes the evolution of a state is generated by the Hamiltonian matrix H, which is defined by Ĥ =X † HX/2 +X † d . The symplectic matrix becomes S = exp [ΩH]. 10 Note that σ is complex in our choice of basis, which implies taking the Hermitian conjugate of S.
The Schrödinger equation can be translated in this language to the simple equation σ f = S σ i S † for the second moments, and r f = S r i for the first moments, which shifts the problem of usually untreatable operator algebra to simple 2N × 2N matrix multiplication. Here, the indices i and f denote the initial and final state, respectively. Finally, Williamson's theorem guarantees that any 2N × 2N Hermitian matrix, such as the covariance matrix σ, can be decomposed as σ = S † ν ⊕ S, where S is an appropriate symplectic matrix. The diagonal matrix ν ⊕ = diag(ν 1 , . . . , ν N , ν 1 , . . . , ν N ) is known as the Williamson form of the state and ν n := coth( ωn 2 kB T ) 1 (where we have introduced normal frequencies ω n and a nominal temper ature T) are the symplectic eigenvalues of the state [46].
Williamson's form ν ⊕ contains information about the local and global mixedness of the state of the system [43]. The state is pure if ν n = 1 for all n and is mixed otherwise. As an example, the thermal state σ th of a N-mode bosonic system is simply given by its Williamson form, i.e. σ th = ν ⊕.

Decoupling of a time evolution operator
The time evolution of a system with time-dependent Hamiltonian Ĥ (t) iŝ where ← T is the time ordering operator. This expression simplifies dramatically when the Hamiltonian Ĥ is time independent, in which case one obtains Û (t) = exp[− i Ĥ t] as a solution to the time-dependent Schrödinger equation. We are, however, interested in Hamiltonians with time-dependent parameters. Any Hamiltonian can be cast in the form Ĥ = n g n (t)Ĝ n , where the Ĝ n are time independent, Hermitian operators and the g n (t) are time-dependent real functions. The choice of Ĝ n need not be unique, and if this is the case, a specific choice is motivated by the specific aims of the problem.
We say that the time evolution operator (3) has been decoupled if it can be written as [28,29]Û where the real functions F n (t) are in general time-dependent. It has been shown that these functions can be found as solutions to a set of differential equations and are determined solely by the parameters g n (t) of the Hamiltonian [28]. The order of the operators in (4) is not unique; a different order changes the form of the functions F n (t), but the not the expectation value of physical quantities. A more detailed outline of these decoupling techniques may be found in appendix A. It is possible to obtain an even more explicit decoupling (4) in the context of Gaussian states and linear (i.e. quadratic in the operators) interactions. Given a set of N bosonic modes, there are N (2 N + 1) independent quadratic Hermitian operators, which we can denote Ĝ n , that can be formed by arbitrary quadratic combinations of the creation and annihilation operators [47] 11 . We also recall that any unitary transformation induced by a quadratic operator, including the quadratic time evolution operator (3), can be represented by a 2N × 2N symplectic 11 For example, where the numbering and ordering of the generators Ĝ n is a matter of convenience. Work in this direction has also been done in [48]. matrix S. Combining all of this together, it can be shown [47] that the symplectic matrix S that represents the time evolution operator (3) takes the form where the symplectic matrices S n are given by S n := exp[F n (t) Ω G n ] and the matrices G n can be obtained through Ĝ n = 1 2 X † G n X, with the restriction that the generator matrix G n must be Hermitian. The techniques to obtain the real, time-dependent functions F n (t) are the same as in the more general case described above. More details can be found in appendix A.

Decoupling algebra of the nonlinear Hamiltonian
Decoupling of the Hamiltonian (1) can be done using different choices of the Hermitian operators Ĝ n . Here, we find it convenient to consider the closed finite 9-dimensional Lie algebra generated by the following set of Hermitian basis operatorŝ which form the smallest set of operators in the Lie algebra that generate the Hamiltonian (2) 12 .
A generic time evolution operator Û (τ ) induced by an arbitrary Hamiltonian cannot in general be written in the form (4) for finite number of operators Ĝ n . A finite decoupling (4) is however possible when the operators forms a finite Lie algebra that is closed under commutation. This is the case for the Hamiltonian in (1), since the commutator of any two elements in the algebra (6) yields a linear combination of the elements of the algebra. This allows us to make an informed ansatz for the evolution operator as we will see below.

Decoupling of a nonlinear time-dependent optomechanical Hamiltonian
In order to achieve the main aim of this work, we need an analytical expression of the decoupling (4) given our Hamiltonian (1). While we will show that we can always obtain a formal expression for the evolution, the coefficients that make up the evolution cannot always be computed analytically, as will be clear for certain choices of the mechanical squeezing function D 2 (τ ).
We find it convenient to proceed by collecting all quadratic terms-including the squeezing term with D 2 in (2)-as a separate operator which we call Û sq (τ ). This choice allows us to study the action of the quadratic and nonlinear parts separately, which can be solved through different means. Since we are interested in computing the first and second moments of the system for the purpose of computing the non-Gaussianity, it is straight-forward to apply Û sq (τ ) to the operator b as a symplectic transformation.
We now make an ansatz for the time-evolution operator Û (τ ) as a finite product of the operators in the algebra: where we have defined an evolution operator Û sq as a quadratic evolution operator of the mechanical degree of freedom: Here, we have effectively divided the Hamiltonian into a quadratic contribution Û sq (τ ) and a remaining nonlinear contribution with the addition of linear term proportional to B ±. The coefficients in the decoupling above can now be obtained in terms of definite integrals. The full calculations can be found in appendix B. We obtain where we have introduced the function and where P 11 (τ ) and P 22 (τ ) are defined below. The only problem that we encounter is computing a decoupled form of Û sq in (8). In fact, it has been shown that decoupling of the evolution operator does not yield analytical results except in very specific cases [49]. For our purposes, this is not problematic, because we can calculate the action of Û sq on the first and second moments analytically using the covariance matrix formalism.

Action of the single-mode squeezing component
Although it is not possible to obtain an analytical decoupling of (8), it is possible to obtain an expression for its action on the operators b and b † . First of all, we note that a Bogoliubov transformation of a single mode operator always has the general expression [49]. The challenge is to find an explicit expression for the Bogoliubov coefficients α(τ ) and β(τ ), which satisfy the only nontrivial Bogoliubov identity |α(τ )| 2 − |β(τ )| 2 = 1. In appendix B we show that the Bogoliubov coefficients α(τ ) and β(τ ) can be obtained through whose explicit form can be obtained once an explicit expression of the functions P 11 (τ ) and P 22 (τ ) is found. Given the previously defined function ξ(τ ) in (10), we also find The two functions P 11 and P 22 are determined by the two following uncoupled differential equations:P where the dot stands for a derivative with respect to τ and the initial conditions are P 11 (0) = P 22 (0) = 1 and Ṗ 11 (0) =Ṗ 22 (0) = 0. Furthermore, the second equation in (12) can be written asÏ which now has boundary conditions I P22 (0) = 0 and İ P22 = 1, and where The solutions to P 11 and P 22 (or I P22 ) can then be used in the expressions (9)-(11) to find the full dynamics of the state. While the solutions must in general be obtained numerically, we anticipate that there are scenarios, such as constant D 2 , where the equations above admit analytical solutions.

Initial state
In this work, we assume that both the optical and mechanical modes are initially in a coherent states, namely |µ c and |µ m respectively, defined as the eigenstates of the annihilation operators, i.e. â|µ c = µ c |µ c and b |µ m = µ m |µ m . For optical fields, this is generally a good assumption. On the other hand, within optomechanical systems the mechanical element is typically found initially in a thermal state. Our choice of initial coherent state can be generalised to that of a thermal state in a straightforward manner, that is, by integrating over the coherent state parameter with an appropriate kernel (as any thermal state may be written as Gaussian average of coherent states, as per its P-representation). Restricting ourselves hence to a single coherent state also for the mechanical oscillator, the initial state |Ψ(0) reads We now proceed to apply (7) to this state.

Full state evolution for general dynamics
For completeness, we present here the full state derived under the evolution with Û (τ ) for two initially coherent states (15): where we have defined K : Note that, in the above, we have kept the dependence on τ implicit but, in general, all exponentials will oscillate in time. We also note that the state (16) contains all terms that have been considered in the literature before, including the contributions from a constant nonlinear light-matter term [7], a timedependent light-matter term [9], and a linear, mechanical displacement term [40]. The main addition here is Û sq , which takes on the trivial form We note here that the expression of (16) allows us to compute the reduced state of the mechanics ρ Mech. (τ ) at any time τ , which readŝ We are now ready to consider the non-Gaussianity of the evolved state.

Measures of deviation from Gaussianity
The time evolution (7) is not linear. Therefore, an initial Gaussian state will evolve, in general, to a non-Gaussian state. Here we ask: is it possible to quantify the deviation from Gaussianity of the state evolving from an initial Gaussian state?
To answer this question we need to find one or more suitable measures of deviation from Gaussianity. In this work we choose to employ a measure based on the comparison between the entropy of the final state and that of a reference Gaussian state [33]. This measure can be understood simply as follows: let us assume that our initial state ρ(0) evolves into the state ρ(τ ) at time τ . We can analytically compute the first and second moments of ρ(τ ). We then consider a Gaussian state, which we call ρ G (τ ), with the same first and second moments as ρ(τ ). It has been shown that this reference state ρ G (τ ) is indeed the Gaussian state that is closest to ρ(τ ) [34]. In general, since ρ(τ ) is not determined uniquely by its first and second moments, as is the case for Gaussian states, the two states do not coincide, i.e. ρ(τ ) =ρ G (τ ).
One way to quantify the difference between two states ρ and ρ G is via the relative entropy S(ρ,ρ G ). It has been shown that the relative entropy S(ρ(τ ),ρ G (τ )) is equivalent to the difference between the local von-Neumann entropies of the states [33]. The measure of non-Gaussianity δ(τ ) can therefore be defined as where S(ρ) is the usual von Neumann entropy of a state ρ, defined by S(ρ) := −Tr(ρ lnρ).
Since the reference state ρ G (τ ) is Gaussian, it is fully characterised by its first and second moments. We therefore turn to the continuous variable formalism and consider the covariance matrix σ of ρ G (τ ). Furthermore, we can define the von Neumann entropy S(σ) of the state as given by S(σ) = j s V (ν j ), where j runs over all the modes, ν j are the symplectic eigenvalues of σ and s V (ν j ) is the binary entropy of the state defined by The symplectic eigenvalues are defined as ν j = |λ j |, where λ j are the eigenvalues of the matrix i Ωσ, where Ω is the 4 × 4 symplectic form. Note that, for all physical states, the eigenvalues satisfy 1 ν j . It follows from the above that a state is non-Gaussian at time τ if and only if δ(τ ) = 0.
An alternative interpretation of this measure is as a quantification of the impurity of ρ G (τ ). While the initial state ρ(τ ) remains pure throughout the evolution, such that S(ρ(0)) = S(ρ(τ )) = 0, the constructed Gaussian reference state ρ G (τ ) does not remain pure. This is not due to external noise, but occurs because we are, loosely speaking, 'approximating' the actual state with the Gaussian subset of states.
In this work, we consider unitary dynamics only. If the initial state ρ(0) is pure at τ = 0, it stays pure throughout its evolution, and the measure thus reduces to where ρ G (τ ) is the Gaussian reference state constructed form the first and second moments of ρ(τ ). Our challenge is therefore to compute the symplectic eigenvalues ν j in order to be able to find the expression of ρ G . Using the expression for the decoupled time evolution operator (8), we can obtain all of the elements of σ. These expressions are cumbersome and can be found in appendix C. The expression for the symplectic eigenvalues are too involved and we choose not to print them. Before we proceed, we also consider the effect of mechanical squeezing on the symplectic eigenvalues. In the continuous variable formalism, a squeezing operation can be represented as a symplectic transformation S acting on the covariance matrix σ through congruence: S σ S † . All symplectic transformations leave the symplectic eigenvalues ν j invariant when acted upon in this way. In this work, however, we consider the inclusion of mechanical squeezing as a term in the Hamiltonian, which acts on the fully non-Gaussian state ρ(τ ). The presence of the nonlinearity means that the squeezing term acts non-trivially on the full state and can actually affect the symplectic eigenvalues of the Gaussian reference state. The mechanical squeezing parameter D 2 (τ ) affects all F-coefficients, meaning that not only the mechanical subsystem but also the optical subsystem will be affected.

Application: non-Gaussianity for optomechanical systems
In this section, we demonstrate the applicability of our techniques by computing the non-Gaussianity of an optomechanical system. The solutions allow us to consider both constant and time-dependent light-matter couplings, however, in order to obtain explicit results we choose to set G (τ ) ≡g 0 constant throughout this work and refer the reader to [9] for a thorough analysis of the non-Gaussianity of the optomechanical state given a time-dependent light-matter coupling. Furthermore, we set D 1 = 0 throughout the remainder of this work. Since the second moments are not affected by a displacement term, the non-Gaussian character of a state remains unchanged [3].
We consider two cases in this section: one where we assume that the mechanical squeezing parameter is constant, and one where the mechanical squeezing is periodic. Our goal is to derive some general bounds on the non-Gaussianity of the state for each case. First, however, we provide some bounds on the total amount of non-Gaussianity.

Bounding the full measure
The exact expression for δ(τ ) is long and cumbersome due to the complex expressions of the covariance matrix elements (C.7). We therefore provide bounds to the measure that can be expressed as simple analytic functions. Since the full measure δ(τ ) is an entropy, it can be bounded from above and below by the means of the Araki-Lieb inequality [50], which reads where ρ AB is the full bipartite state and ρ A and ρ B are the traced-out subsystems. This inequality allows us to bound the behaviour of the full measure δ(τ ) in terms of the subsystem entropies.
We therefore proceed to define the lower and upper bounds as In our case, the subsystems are the traced out optical state ρ Op and the traced out mechanical state ρ Me . To quantify the entropy of the subsystems, we must find the symplectic eigenvalues of the optical and mechanical subsystems, which we call ν Op and ν Me respectively. Lengthy algebra (see appendix C), the use of the Bogoliubov identities |α| 2 = 1 + |β| 2 and α β * = α * β, and observing that |EB and its appearance in the first and second moments) allow us to find where we recall that KN a := FN aB− + i FN aB+ , and where we have defined θ(τ ) = 2 FN2 a + FN aB+ FN aB− . The optical symplectic eigenvalue (21) is bounded by which can be inferred by noting that KN a is generally given by an oscillating function multiplied by the strength of the optomechanical coupling g 0 . For specific τ which ensures that |KN a | 2 = 0, and then considering g 0 1, the exponentials in ν Op in (21) are suppressed, which means we are left with , the bipartite entropy of the Gaussian reference state S(ρ G ) is approximately equal to one of the subsystem entropies. To determine when this is the case, we consider the maximum values of ν Op and ν Me . In general, when |µ c | 2 1, and when |KN a | 2 1, which requires g 0 1 and specific values of τ , the eigenvalues ν Op and ν Me tend to their maximum values ν Op,max and ν Me,max , which are We note that there are three distinct scenarios which arise from the comparison of the coherent state parameter |µ c | 2 and the function KN a : Here, the non-Gaussianity is well-approximated by (ii) Secondly, we assume In this case, the Araki-Lieb bound is not very informative since the left-hand-side is zero and must evaluate the non-Gaussianity exactly.
Note that the first two cases might occur only for short periods of time τ since KN a is oscillating. Furthermore, we note that the squeezing parameter D 2 (τ ) affects the peak value of the non-Gaussianity because it enters into |KN a | through the F-coefficients (9). The dependence is non-trivial, but we will consider the analytic case for constant squeezing below. However, in general, when |µ c | |KN a |, we see from (24) that the non-Gaussianity is independent of D 2 (τ ) and can be accurately modelled by the standard optomechanical Hamiltonian without mechanical squeezing.
Let us now consider two specific cases where the squeezing term is either constant or modulated.

Applications: constant squeezing parameter
Here we assume that the rescaled squeezing parameter is constant, with D 2 (τ ) =d 2 . This case is equivalent to the case where the mechanical oscillation frequency ω m is shifted by a constant amount and where the initial state is a squeezed coherent state, see appendix D. We begin by deriving analytic expressions for the coefficients in (9) given this choice of parameters.

Decoupled dynamics.
We use the methods discussed in section 3 to start by solving the differential equation (12). We find the solutions P 11 = P 22 = cos ζτ, where we define ζ := 1 + 4d 2. This, in turn, yields the following Bogoliubov coefficients (defined in (11)): Furthermore, we find ξ(τ ) = cos ζτ − i ζ sin ζτ , which in turn can be integrated to obtain the coefficients (9), which now read where sinc(x) = sin(x)/x. Since D 1 = 0, all other coefficients are zero. The functions (27) now fully determine the time evolution through (7).

Quadratures.
To gain intuition about the evolution of the system, we include plots of the optical quadratures. These can be found in figure 2. The quadratures are the expectation values of x 1 = (â † +â)/ √ 2 and p 1 = i(â † −â)/ √ 2 and would correspond to classical trajectories in phase space. The full expression for the expectation values x 1 and p 1 can be found in (C.5) in appendix C.
In figures 2(a)-(d), we have plotted the quadratures for µ c = 1, µ m = 0, g 0 = 1 and increasing values of d 2 . While it is generally difficult to engineer a coupling of this magnitude, these To show the directionality of the evolution, the colour of the curve starts as light blue for τ = 0 and becomes increasingly darker as τ increases. We observe that the addition of mechanical squeezing causes the system to trace out highly complex trajectories, compared with the case when d 2 = 0.

Measure of non-Gaussianity.
We now proceed to compute the non-Gaussianity δ(τ ), defined in (18), of the state evolving at constant squeezing parameter. A fully analytic expression for δ(τ ) exists but is too cumbersome to include here. Instead, we plot the measure of non-Gaussianity in figure 3. In the first row of figure 3, we present a comparison between the full measure δ ( figure 3(d)) and the lower and upper bounds δ min and δ max provided by the Araki-Lieb inequality in figures 3(e) and (f).
We note that the non-Gaussianity increases for large light-matter coupling g 0 and large coherent state parameter µ c . This feature was also observed for standard optomechanical systems in [9]. However, the most striking feature here is that the larger d 2 is, the less non-Gaussian the system becomes. To understand why this is the case, we examine the dependence on For large d 2 , and therefore large ζ, the first term inside the brackets dominates and for ζτ = nπ with integer n, we are left with |KN a | ∼g 0 sin 2 (ζ τ)/ζ 2 . In general, we find limd 2 →∞ |KN a | 2 = 0. The consequences for the non-Gaussianity are difficult to predict given the complexity of the expressions, but we note that the mechanical symplectic eigenvalue ν Me decreases, while the optical symplectic eigenvalue ν Op increases.
Furthermore, the quantity θ(τ ) = 2 FN2 a + FN aB+ FN aB− is given by We find that limd 2 →∞ θ(τ ) = 0. We then look at the symplectic eigenvalues (21) in this limit. We find that ν Me → 1, and ν Op → 1, which means that both the upper and the lower bounds of the non-Gaussianity tend to zero, and hence δ(τ ) → 0 as d 2 increases. We conclude that increasing the amount of constant squeezing in the system reduces the overall non-Gaussianity of the state.

Applications: modulated squeezing parameter
In this section, we consider a modulated squeezing term. The dimensionless squeezing D 2 (τ ) = D 2 (t)/ω m is time-dependent and of the form where d 2 = d 2 /ω m is the amplitude of the squeezing and Ω 0 denotes the time-scale of squeezing 13 . The differential equations in (12) are not generally analytically solvable for arbitrary choices of D 2 (τ ). However, for the choice of D 2 (τ ) in (30), both equations have a known form. Consider the differential equation for P 11 , which we reprint here for convenience, Equation (31) is that of a parametric oscillator, which is used elsewhere in physics to describe, for example, a driven pendulum. As shown in appendix B, the equation for the integral of P 22 (B.19) takes the same form. The equation (31) is known as the Mathieu equation. In its most general form, and using conventional notation, it reads: where a, q, and x are real parameters. The general solutions to this equation are linear combinations of functions known as the Mathieu cosine C(a, q, x) and Mathieu sine S(a, q, x), the exact form of which will be determined by the boundary conditions for y .
To find which values the a, q and x parameters correspond to, we note that the cosine-term in D 2 (τ ) has the argument Ω 0 τ, which means that we must rescale time τ as τ = Ω 0 τ /2. Inserting our expression for D 2 (τ ) and using the chain-rule to change variables from τ to τ , we rewrite the equation for P 11 as where we identify the variables a = 4/Ω 2 0 , and q = −8d 2 /Ω 2 0 . The boundary conditions P 11 (0) = 1 and Ṗ 11 (0) = 0 will yield the Mathieu cosine C(a, q, x), and for I P22 as the solution, and the boundary conditions I P22 (0) = 0 and İ P22 (0) = 1 will yield the Mathieu sine S(a, q, x) as the solution. For our specific choice of D 2 (τ ) in (30), the system is resonant at Ω 0 = 2 (see appendix E), which means that a = 1 and q = −2d 2 .

Approximate solutions at resonance
The Mathieu equations are notoriously difficult to evaluate numerically. Instead, we use a twoscale method to derive perturbative solutions to P 11 and I P22 . The perturbative solutions are valid for d 2 1 and make use of specific resonance conditions to ensure that the solutions do not diverge. See appendix E for the full derivation, where we also show that these approximate solutions correspond exactly to the more physically intuitive RWA when τ 1. For smaller values of τ , the approximate solutions are still valid, but they cannot be interpreted as equivalent to the RWA.
The squeezing term is resonant when Ω 0 = 2. We find that the approximate solutions for P 11 and I P22 (the integral of P 22 ) are given by, respectively, We then compute ξ(τ ) in (E.15). We assume that d 2 τ 1 to find where we have expanded the hyperbolic functions to second order. By using the relations between ξ(τ ) and the Bogoliubov coefficients (B.30), we find that the Bogoliubov condition is approximately satisfied as: With this expression, we can now compute the non-zero F-coefficients (9): where we have discarded terms with d 3 2 . With these expressions, we are ready to compute the non-Gaussianity of the system when the squeezing is applied at mechanical resonance.

Measure of non-Gaussianity at resonance.
We first compute the full measure of non-Gaussianity δ(τ ) and plot the results in figure 4. Figure 4(a) shows the full measure δ(τ ), the lower bound δ min (τ ) and the upper bound δ max (τ ) as a function of τ for the parameter g 0 = 1, µ c = 1, d 2 = 0.1, and µ m = 0 as a function of time τ and the squeezing d 2 . The second plot in figure 4(c) also shows the full measure δ(τ ), the lower bound δ min (τ ) and the upper bound δ max (τ ) as a function of g 0 at τ = π, d 2 = 0.1, µ c = 1, and µ m = 0. We find that the non-Gaussianity increases with g 0 , as expected.
In figure 4, we considered d 2 = 0.1; a value consistent with the validity of the approximate solutions to the Mathieu equation. For this value, the non-Gaussianity is found to increase very slightly with d 2 . To demonstrate this, we consider the regime where 1 2|µ c | |KN a |, which occurs when 2|µ c | g 0 for specific values of τ . In this regime, the non-Gaussianity was approximately given by s V 2 |KN a ||µ c | (25). Given the functions (37), we find that |KN a | 2 = 4g 2 0 sin 2 (τ /2) +g 2 where we have again removed terms proportional to d 3 2 and d 4 2 . The behaviour of |KN a | 2 is markedly different compared with the constant case. Firstly, while |KN a | 2 still oscillates, it also increases with τ and d 2 . If we consider the leading term with τ 2 , we find that the non-Gaussianity scales with δ ∼ ln(τd 2g0 ), which confirms that in this specific regime, the non-Gaussianity increases logarithmically with τ , d 2 , and g 0 . We conclude that squeezing is not necessarily detrimental to the non-Gaussianity if the squeezing is modulated at resonance, although more work needs to be done to ascertain the full interplay between the two effects.

Discussion
Before presenting our conclusions, we discuss the applicability and scope of the techniques we developed. We also comment on the effect of squeezing on the non-Gaussian character of the system.

Advantages over direct numerical simulations
With our techniques, we have shown that it is possible to analytically solve the dynamics of a nonlinear optomechanical system even when the mechanical squeezing is time-dependent. To emphasise this point, we wish to compare our approach, which relies on numerically solving the differential equations in (12), with a general numerical method using a standard higherorder Runge-Kutta solver to evolve a state in a truncated Hilbert space, e.g. using the Python library QuTiP [51].
When the dynamics is solved with a Runge-Kutta method, the continuous variable (pure) states are represented as finite-dimensional vectors in a truncated Hilbert space. When the system is nonlinear, information about the state is quickly distributed across large sectors of the Hilbert space. If the computational Hilbert space is too small, numerical inaccuracies quickly enter into the evolution, as a result of truncating the vectors. It follows that the dimension of the Hilbert space must be large enough to prevent this, which requires significant amounts of computer memory. It is also very difficult to consider parameters of the magnitude g 0 = 10 and d 2 = 10, as done in this work, since these cause the system to evolve very rapidly and, consequently, require the evolution of the system to be calculated using smaller and smaller time intervals.
The methods developed here excel at treating systems numerically for large parameters g 0 ,d 2 , µ c and µ m . However, we note that it becomes increasingly difficult to numerically evaluate the dynamics at longer times τ when the system of differential equation (12) is numerically solved for arbitrary functions D 2 (τ ). The difficulty is primarily caused by the double integral that determines the coefficient FN2 a in (9), which must be evaluated numerically. For each value of τ , the integral will be evaluated from 0 to the final τ , and then from 0 to τ . As a result, the integrals take an increasingly long time to evaluate for large τ . We therefore conclude that the key strength in our method lies in evaluating the state of the system at early times τ ∈ (0, 2π) for large parameters µ c ,g 0 , and d 2 . We also emphasise that, the computation using our methods is numerically exact, which a naive computation using QuTiP or a similar library is not.
To conclude, our methods allow for the evaluation of the state of the system with large parameters, e.g. g 0 = 100 and d 2 = 10, which would be nearly impossible to perform with QuTiP or a similar library unless one had access to significantly more computational resources.

Competing behaviours of nonlinearity and squeezing
We concluded from figure 3 that the addition of a constant squeezing term has a detrimental effect on the non-Gaussianity of the system. We also noted that inclusion of a constant squeezing term is equivalent to changing the mechanical trapping frequency ω m to a specific value and starting the computation with a squeezed coherent state (see appendix D). With this interpretation, our results also show that an initially squeezed state evolving under the optomechanical Hamiltonian can be expected to exhibit less non-Gaussianity compared with coherent states. The reason for this overall behaviour can be found by simple inspection of the total Hamiltonian. If a strong squeezing term is included in the Hamiltonian (1), it dominates over the interaction term, leading to a decrease in the non-Gaussianity. However, such a process is not fully monotonic, since an increase of the squeezing parameter does not always decrease the non-Gaussianity. This is, however, reasonable, as it cannot be expected that only the relative weight of the two parts of the Hamiltonian matter; the precise dynamics is much more complex, and the non-Gaussianity depends on the entire state, which is driven by the full Hamiltonian.
The finding that the non-Gaussianity increases with both time τ and d 2 when modulated at mechanical resonance is interesting and warrants further investigation. We leave this to future work.

Conclusion
In this work, we solved the time-evolution of a nonlinear optomechanical system with a timedependent mechanical displacement term and a time-dependent mechanical single-mode squeezing term. We found analytic expressions for all first and second moments of the quadratures of the nonlinear system and used them to compute the amount of non-Gaussianity of the state. We considered both constant and modulated squeezing parameter, and found that a squeezing parameter modulated at twice the mechanical resonance results in the Mathieu equations, for which we provide approximate solutions equivalent to the RWA.
In general, we find that the relationship between the squeezing and non-Gaussianity is highly nontrivial. The inclusion of a mechanical squeezing term in the Hamiltonian, which is equivalent to starting with a coherent squeezed state evolving with the standard optomechanical Hamiltonian with a shifted mechanical frequency, decreases the overall non-Gaussianity of the state. If the squeezing term is modulated at twice the mechanical resonance, however, we found that the non-Gaussianity increases with both time and the squeezing parameter in specific regimes. These results hold interesting implications for quantum control of nonlinear optomechanical systems.
Our results also suggest that the combination of non-Gaussian resources and mechanical squeezing may not necessarily be beneficial if the application relies specifically on the non-Gaussian character of the state. However, more work is needed to conclude if this has a significant effect on, for example, sensing applications. More work is also necessary to properly study the instabilities of the full solutions to the Mathieu equations and how they affect the dynamics. The effect of squeezing the optical rather than mechanical mode is another question we defer to future work.
The decoupling methods demonstrated here constitute an important step towards fully characterising nonlinear systems with mechanical squeezing and can be used both to model experimental systems and to test numerical methods. The results presented here apply to any system with the same characteristic cubic Hamiltonian interaction term and single mode squeezing term. More broadly, the Lie algebra method can be used to solve any system where a finite set of Lie algebra operators has been identified. Our work can further be extended to more complicated quadratic Hamiltonians of bosonic modes, such as Dicke-like models [52], which would allow for applications in other areas of physics to be developed. SQ acknowledges support from the EPSRC Centre for Doctoral Training in Delivering Quantum Technologies and thanks the University of Vienna for its hospitality. DR would like to thank the Humboldt Foundation for supporting his work with their Feodor Lynen Research Fellowship. This work was supported by the European Union's Horizon 2020 research and innovation programme under Grant agreement No. 732894 (FET Proactive HOT). DEB acknowledges support from the CEITEC Nano RI.

Appendix A. Decoupling of techniques for time evolution
In this appendix, we outline the general decoupling techniques that we shall be using throughout this work to find a decoupled time-evolution operator generated by the Hamiltonian in (2).

A.1. Decoupling for arbitrary Hamiltonians
The time evolution operator Û (t) induced by a Hamiltonian Ĥ (t) readŝ Any Hamiltonian can be cast in the form Ĥ = n g n (t)Ĝ n , where the Ĝ n are time independent, Hermitian operators and the g n (t) are time-dependent functions. The choice of Ĝ n need not be unique. It has been shown [28,47] that it is always possible to obtain the decouplinĝ where we have defined Û n := exp[−i F n (t)Ĝ n ] and the real, time-dependent functions F n (t), and the ordering of the operators is Û 1Û2 . . .. The functions F n (t) are uniquely determined by the coupled, nonlinear, first order differential equations where Ĥ and all F n are taken at time t. This is the general method we will be employing in the following sections.

A.2. Decoupling for quadratic Hamiltonians
If the Hamiltonian is quadratic in the mode operators the techniques described in the previous subsection have a more powerful representation. As explained in the main text, the time evolution operator has a symplectic representation S(t) and the decoupled ansatz (A.2) has the form S(t) = n S n (t), where we have introduced S n := exp[F n (t) Ω G n ] and the real, timedependent functions F n (t) are the same as those obtained with the technique above. The real, time-dependent functions F n (t) can be obtained by solving the following system of coupled nonlinear first order differential equations where the matrix H can be obtained by Ĥ (t) = 1 2 X † H X and the summation is over N (2N + 1) terms [47]. This is the matrix version of the operator differential equation (A.3) for quadratic Hamiltonians, which reduces the problem of operator algebra to matrix multiplication.

Appendix B. Decoupling of the nonlinear Hamiltonian
We use the techniques presented in appendix A to decouple the Hamiltonian (2). The decoupling below is obtained in the same fashion of previous results in the decoupling of Hamiltonians [47,49].

B.1. Solving the quadratic time-evolution
We start from the dimensionless Hamiltonian (2), which we reprint herê This Hamiltonian can be conveniently re-written aŝ Our goal now is to separate the Hamiltonian into a linear and a quadratic contribution. We then study the effect of the quadratic Hamiltonian on the nonlinear part, which allows us to focus fully on decoupling the nonlinear dynamics. We therefore rewrite the time evolution operator (3) in the alternative form where we have introduced The transition to the expression (B.3) is similar to moving from the Heisenberg (or Schrödinger) picture to the interaction picture. Our first goal is to solve the time-evolution of the quadratic part Û sq . We define the basis vector X := (b,b † ) T . From standard symplectic (i.e. Bogoliubov) theory we know that where the 2 × 2 symplectic matrix S sq (τ ) is the symplectic representation of Ĥ sq (τ ) and satisfies S † sq (τ ) Ω S sq (τ ) = Ω.
Here Ω = diag(−i, i) is the symplectic form.
The matrix S sq (τ ) therefore has the expression S sq (τ ) = Here Ĥ sq = 1 2 X †H sq X and one has that Therefore, we have that It can be shown that the time independent orthogonal matrix M ort with expression puts the Hamiltonian matrix H s in diagonal form through (B.9) This allows us to manipulate (B.7) and find This means that we have We introduce the new vector Y := M ort X, which is just a rotation of the operators. Then we have

B.2. Solving the matrix time-ordered exponential
Here we seek a formal expression for (B.12). We write in terms of the matrix K defined as and the matrix P, which we will determine and which is diagonal. This follows from the fact that the matrix on the left-hand side of (B.13) is diagonal when squared, and therefore any even powers in the expansion will be diagonal. We use the fact that to find the equation Since K is invertible, we manipulate this equation and obtain, after some algebra, We can now solve the four differential equations contained in (B.17) which read The differential equation (B.18) must be supplemented by initial conditions. We note that since the left hand side of (B.13) is the identity matrix for τ = 0 we have that P(0) = 1 which implies P 11 (0) = P 22 (0) = 1. In addition, taking the time derivative of both sides of (B.13) and setting τ = 0 implies Ṗ 11 (0) =Ṗ 22 (0) = 0. By introducing the integral I P22 = τ 0 dτ P 22 (τ ), one can also rewrite the second equation as so that it becomes the same as that for P 11 . The boundary conditions are now I P22 (0) = 0 and İ P22 = 1. We were not able to find a general solution to the differential equation for P 11 (B.18) and the equation for I P22 (B.19), but they can be integrated numerically when an explicit form of D 2 (τ ) is given. For specific choices of D 2 (τ ), which we discuss in the main text of this work, the equations become the well-studied Mathieu equation.
Proceeding with the solution to the quadratic time-evolution, we have In turn, this allows us to get where we have introduced the Bogoliubov matrix S sq (τ ) with coefficients α(τ ) and β(τ ).
This gives us, after a little algebra, These quantities can also be written in terms of I P22 as This means that the basis vector X (B.5) transforms as As a side remark, the Bogoliubov (symplectic) identities |α(t)| 2 − |β(t)| 2 = 1 read

(B.30)
We now note that the operators N a ,N 2 a ,B + ,B − ,N aB+ ,N aB− form a closed Lie sub-algebra of the full algebra of our operators. Therefore, we can apply the decoupling techniques described above.
We proceed to decouple the remaining part of the operator (B.28), which readŝ (B.31) We make the ansatẑ (B.32) Taking the time derivative on both sides and arranging in a similar fashion to (A.3) we find (B.33) Therefore our main differential equations can be obtained by equating the coefficient of the different, linearly independent operators of the Lie algebra from the equation below (B.34) The solutions with operators proportional to B ± can be independently solved. However, the solutions for FN a , FN2 a , and FN aB± are less straight-forward. We identify the following four coupled differential equations: . We find the following key expression for this work to be complemented with the functions listed in (B.36).

Appendix C. First and second moments and covariance matrix elements
We can employ the results summarised in section 3 to compute all the relevant quantities needed in this work. They include the first and second moments of the state of the system at any moment in time, which we use to compute the measure of non-Gaussianity, see section 4.
To do this, we use the explicit form of the time-evolution operator, written out above in (B.37). We further assume that the initial state is a separable coherent state |µ c ⊗ |µ m , defined in (15). In order to compute the second moments that constitute the covariance matrix σ for the state, we must calculate the expectation values of â , b , â † a , b †b , â 2 , b 2 , âb , and âb † . To achieve this goal, we start by defining the following quantities for ease of notation: and where we also introduce which we use when deriving the state (16). The last two quantities can be computed using the expression for the Weyl displacement operator D α = e αb † −α * b and the combination D αDβ = e (αβ * −α * β)/2D α+β , and read 3) The time-evolution of the operators â and b in the Heisenberg picture is thereforê

(C.4)
These expressions allow us to compute the expectation values of the first and second moments given our initial state. We have here transformed into a frame rotating with e −i Ωc τ , which read The two-mode covariance matrix is fully determined by its elements σ nm , defined in this basis as where all the other elements follow from the symmetries of σ, imposed by the requirement that σ † = σ. Given all of the above, we can compute an explicit expression for elements (C.6) of the covariance matrix, which reads (C.7)

C.1. Symplectic eigenvalues of the optical and mechanical subsystems
In this appendix, we use the covariance matrix elements to compute the symplectic eigenvalues of the optical and mechanical subsystems.
( C.11) Putting them together, we find which can be simplified into Next, we compute the mechanical symplectic eigenvalue.

(D.9)
To cancel the term proportional to b †2 +b 2 , we have to fix D 2 e 2r = −ω m cosh(r) sinh(r) = −ω m (e 2r − e −2r )/4, and therefore, = ω câ †â + ω m (1 + 2 sinh 2 (r))b †b + sinh 2 (r) + D 2 e 2r 2b †b + 1 (D.10) In particular, Ĥ 0 →Ĥ 0 − (ω m − ω m )/2 under the replacement ω m → ω m . Furthermore, we find that b +b † transforms aŝ When applying the same transformation to the nonlinear part of the Hamiltonian, we find If G(t) ∝ 1/ √ ω m and D 1 (t) ∝ 1/ √ ω m , which is indeed fulfilled for the interesting cases, we find that Ĥ opt →Ĥ − (ω m − ω m )/2 under the replacement ω m → ω m , wherê (D.13) We define Ĥ ω m opt :=Ĥ opt [ω m → ω m ], and we obtain for the full time evolution (D.14) For the expectation values of quadrature of a state ρ, this leads to where ρ sq :=Û † sqρÛsq . For the initial separable coherent state |µ c |µ m , we find that the time evolution of the quadratures induced by the full Hamiltonian Ĥ with constant D 2 can be obtained by calculating the corresponding time evolution of the quadratures induced by Ĥ with vanishing D 2 by replacing ω m with ω m and considering the squeezed coherent initial state Û † sq |µ c |µ m = |µ c |µ m , r , where µ m = µ m cosh(r) + µ * m sinh(r). As a result, the techniques we have developed here can also be utilised to model all the expectation values for an optomechanical system for an initially squeezed states |D 2 (z) .

Appendix E. Solutions to the differential equations at resonance
In this appendix, we first derive perturbative solutions to the Mathieu equations and then compare them with dynamics that we obtain from performing the RWA. We present two approaches which both amount to the same approximation while starting from different assumptions.

E.1. Approximate two-scale solution
Our goal is to obtain approximate solutions to the differential equations (B.18) and (B.19). We will do so by following the derivation in [53] with some modifications.
The general form of Mathieu's equation is given by We will use the general notation in this appendix and then compare with the system in the main text. We begin by defining a slow time scale X = qx. We then assume that the solutions y depend on both scales, such that y(x, X). This means that we can treat x and X as independent variables and the absolute derivative d/dx in (E.1) can be split in two: Mathieu's equation (E.1) therefore becomes (∂ x + q ∂ X ) 2 y(x, X) + (a − 2q cos(2x)) y(x, X) = 0 .

(E.3)
We then expand the solution y(x, X) for small q as y(x, X) = y 0 (x, X) + q y 1 (x, X) + O(q 2 ) and insert this into the differential equation above. Our goal is to obtain a solution for y 0 which incorporates a number of restrictions from the differential equation for y 1 .
To zeroth order, we recover the regular harmonic differential equation for y 0 , which is the limiting case as q → 0: where we know that the solutions are sinusoidal, while the coefficients must depend on X. We choose the following trial solution: Our goal is now to find explicit solutions to the complex function A(X). We continue with the equation for y 1 . We discard all terms of order q 2 to find q ∂ 2 x y 1 + 2 q ∂ x ∂ X y 0 + a q y 1 − 2 q cos(2x)y 0 = 0 . (E.6) We divide by q and insert our solution for y 0 to find At this point, we specialise to a = 1, which corresponds to setting Ω 0 = 2 in the main text. We combine the exponentials to find (E. 8) In order for the solution to be stable, we require that secular terms such as resonant terms e ix vanish. If these do not vanish, the solution will grow exponentially [53]. We also neglect terms that oscillate much faster, such as e 3ix . This leaves us with the condition that 2i ∂A * (X) ∂X + A(X) = 0 , (E.9) which can be differentiated again and solved with the trial solution A(X) = (c 1 − i c 2 ) e X/2 + (c 3 − i c 4 ) e −X/2 for the parameters c 1 , c 2 , c 3 and c 4 . From the requirement in (E.9), it is now possible to fix two of the coefficients in (E.10). We differentiate A(X) and use (E.9) to find that the conditions c 1 = c 2 and c 3 = −c 4 must always be fulfilled. We then recall that X = qx and after combining some exponentials, we obtain the full trial solution for the zeroth order term y 0 : y 0 (x) = A(qx) e i x + A * (qx) e −ix = 2 g 0 c 1 e qx/2 + c 3 e −qx/2 cos(x) + 2 g 0 c 1 e qx/2 − c 3 e −qx/2 sin(x) .
(E.10) We now proceed to compare this solution with the parameters and initial conditions given for P 11 in (B.18) and I P22 in (B.19) in the main text, which are both solved by the Mathieu equation.
First, we note that q = −2d 2 and that x = τ. Then we consider the boundary conditions for P 11 , which are P 11 (0) = 1 and Ṗ 11 (0) = 0. From these conditions, we find that g 0 c 1 = c 3 = 1/4, and the the approximate solution to P 11 is given by P 11 (τ ) = cos(τ ) cosh(d 2 τ ) − sin(τ ) sinh(d 2 τ ) . (E.11) The equation for I P22 has the opposite initial conditions I P22 (0) = 0 and İ P22 = 1. For this case, we find that g 0 c 1 = −c 3 = 1/(4 (1 −d 2 )). The full solution to I P22 is therefore cos(τ ) sinh(d 2 τ ) − sin(τ ) cosh(d 2 τ ) , (E.12) and thus Both solutions reduce to the correct solutions for the zero-squeezing case as d 2 → 0. The validity of the pertubative approach can be determined as follows. By inserting the trial solution for P 11 (τ ) into the Mathieu equation, we are left with terms that are multiplied by d 2 . The leading term is in fact d 2 cosh(d 2 τ ). These terms must be approximately zero to solve the Mathieu equation, which means that we require d 2 cosh(d 2 τ ) 1. This means that while d 2 1, we can allow d 2 τ ∼ 1, which means that τ can be large provided that d 2 is sufficiently small.
From the expression for ξ(τ ) in (B.29) we then find (E.14) For very small d 2 1, which was the condition for deriving the approximate solutions in the first place, we can approximate the fraction as unity and we find the compact expression ξ(τ ) = e −i τ cosh(d 2 τ ) + i e i τ sinh(d 2 τ ) .
(E. 15) To better understand what this approximation entails physically, we compare it with the RWA, which has a well-known physical interpretation.