Enhanced nonlinearity through optomechanical modulation

We study the non-Gaussian character of quantum optomechanical systems evolving under the fully non-linear and time-dependent optomechanical Hamiltonian. By using a measure of non-Gaussianity based on the relative entropy of an initially Gaussian state, we quantify the amount of non-Gaussianity induced by the characteristic cubic light–matter coupling and study its general and asymptotic behaviour. We find analytical expressions for the measure of nonlinearity in specific cases of interest, such as for constant and modulated light– matter coupling. More importantly, we also show that it is possible to further increase the amount of non-Gassuianity of the state by driving the trilinear light– matter coupling at the frequency of mechanical resonance, suggesting a viable mechanism for increasing the non-Gaussianity of optomechanical systems. ar X iv :1 81 2. 08 87 4v 2 [ qu an tph ] 7 J an 2 01 9 Enhanced nonlinearity through optomechanical modulation 2


Introduction
Understanding non-linear, interacting physical systems is paramount across many areas in physics. Specifically, "non-linear" (or "anharmonic") dynamical systems include all those whose Hamiltonian cannot be expressed as a second-order polynomial in the quadrature operators. Therefore, these systems allow us to exit the realm of Gaussian states, which are the ground and thermal states of second-order Hamiltonians, and which are also preserved by such Hamiltonians.
In the context of quantum information and computation, there has been a drive towards the realisation of anharmonic Hamiltonians as well as more general methods and control schemes capable of generating and stabilising non-Gaussian states. On the one hand, this is motivated by the fact that, in order to obtain effective qubits from the truncation of infinite dimensional systems, one needs unevenly spaced energy levels, such that only the transition between the two selected energy levels may be targeted and driven. In turn, this requires a sufficiently anharmonic Hamiltonian. On the other hand, it has always been clear that protocols entirely restricted to Gaussian preparations, manipulations and read-outs, through quadratic Hamiltonians and general-dyne detection, are classically simulatable, as their Wigner functions may be mimicked by classical probability distributions. ‡ In recent years, the intrinsic value of non-linear systems, as opposed to the aforementioned limitations that linear systems face, has been made clearer and more rigorous. Non-linearities in the form of non-Gaussian states constitute an important resource for quantum teleportation protocols [1], universal quantum computation [2,3], quantum error correction [4], and entanglement distillation [5,6,7]. This view of non-Gaussianity as a resource for information-processing tasks has inspired recent work on developing a resource theory based on non-Gaussianity [8,9,10]. In addition, it has been found that non-Gaussianity provides a certain degree of robustness in the presence of noise [11,12].
In optomechanical systems [13], where electromagnetic radiation is coherently coupled to the motion of a mechanical oscillator, the light-matter interaction induced by radiation pressure is inherently non-linear [14,15,16]. Optomecahnical systems also have been studied in the context of force sensing [17] and gravimetry [18,19] due to their excellent sensing capabilities. Furthermore, their comparatively large mass make them ideal for testing fundamental phenomena, such as collapse theories [20] and, possibly, signatures of gravitational effects on quantum systems at low energies [21,22].
A number of different optomechanical systems have been experimentally implemented, including Fabry-Pérot cavities with a moving-end mirror [23], levitated nano-diamonds [24,25], membrane-in-the-middle configurations [26] and optomechanical crystals [27,28]. However, most experimental settings can be fully modelled with linear dynamics [29,30], and it is generally difficult to access the fully non-linear regime. While several experiments have demonstrated genuine non-linear behaviour (see for example [31,32,30]), significant effort has been devoted to the question of how the non-linearity can be further enhanced. Most approaches focus on the few-excitation regime, where increasing the inherent light-matter coupling allows for detection of the non-linearity. This enhancement can be achieved, for example, ‡ We should note here that, since uncertainties in quantum Gaussian systems are fundamentally bounded by the Heisenberg uncertainty relation, which in principle does not hold for classical systems, Gaussian operations are in fact sufficient to run some protocols requiring genuine quantum features, such as continuous variable quantum key distribution. by using a large-amplitude, strongly detuned mechanical parametric drive [33], or by modulating the optomechanical coupling [34]. Similar work has shown that the inclusion of a mechanical quartic anharmonic term can be nearly optimally detected with homodyne and heterodyne detection schemes, which are standard measurements implemented in the laboratory [35].
A natural question that arises considering the approaches above is: Are there additional methods by which the amount of non-Gaussianity in an optomechanical system can be increased further ? In this work, we address precisely the question of developing a new method to fully characterise the non-Gaussianity of the nonlinear optomechanical state for an ideal closed system, and to increase it. More precisely, given an initial Gaussian state evolving under the standard optomechanical Hamiltonian, we quantify how non-Gaussian the state becomes as a function of time and the parameters of the Hamiltonian in question. To do so, we make use of recently developed analytical techniques to study the time-evolution of timedependent systems [36], and employ a specific measure based on the relative entropy of the state [37]. This measure has previously been extensively used to compute the non-Gaussianity of various states [38], as well as in an experimental setting where single photons were gradually added to a coherent state to increase its non-Gaussian character [39]. Several additional measures for the quantification of non-Gaussianity have been proposed in the literature, linking it to the Hilbert-Schmidt distance [40] or to quantum correlations [41]. Furthermore, a connection has been put forward between the non-Gaussianity of a state and its Wigner function [42], and similarly there appears to be an intrinsic link between the quantum Fisher information and the lowest amount of non-Gaussianity of a state [43].
Our results include the fact that the non-Gaussianity of an optomechanical system initially in a coherent state scales with the inherent light-matter coupling as expected. We also analyse how it scales with the coherent state parameter of the optical system, and illustrate how this behaviour differs for small and large coherent states. We find that in the regime of small initial coherent state amplitude, the non-Gaussianity grows at least polynomially with the light-matter coupling, while for a large initial coherent state amplitude the non-Gaussianity grows at most logarithmically. This asymptotic behaviour does not depend on the specific form of the non-linearity and is a consequence of the measure we use being based on the von Neumann entropy. Therefore, the scaling itself is not an indication of a fundamental property of the non-Gaussianity. Most importantly, however, we find that the non-Gaussianity of the state can be further increased by modulating the light-matter coupling strength at the mechanical resonance frequency. Since such modulation might be achievable in certain set-ups, this realisation may lead to an additional method to enhance and observe non-linear behaviour in the laboratory. This paper is organised as follows. In Section 2 we present the Hamiltonian of interest and solve the resulting dynamics. In Section 3 we present the measure of non-Gaussianity from [37] and compute both an exact and a more intuitive approximate expression. We proceed then to examine the behaviour of the non-Gaussianity for two cases: a constant light-matter coupling in Section 5; and a time-dependent coupling in Section 6, where we also show that driving the coupling results in larger amounts of non-Gaussianity. Finally, in Section 7 we discuss our results and relate them to a number of applications taken from the literature. Section 8 concludes this work.

Dynamics
We start our work by presenting the necessary mathematical tools needed to solve the dynamics. We refer the reader to Appendix A for the full derivation.

Hamiltonian
We begin by considering two bosonic modes corresponding to an electromagnetic mode and the other to the mechanical motion of a harmonic oscillator. Without loss of generality, we shall henceforth refer to the electromagnetic mode as the optical mode. The operatorsâ andb of the cavity and mechanical modes respectively, obey the canonical commutation relations [â,â † ] = [b,b † ] = 1, while all other commutators vanish. The radiation pressure induces a non-linear interaction and the whole systems is modelled by the Hamiltonian where ω c and ω m are the frequencies of the two modes and g(t) drives the, potentially time-dependent, non-linear light-matter coupling. This light-matter coupling strength g(t) takes on different functional forms for different optomechanical systems, and we also note that this Hamiltonian governs the evolution of many similar systems, including electro-optical systems [44]. To simplify our notation and expressions, we rescale the laboratory time t by the frequency ω m , therefore introducing the dimensionless time τ := ω m t, the dimensionless frequency Ω := ω c /ω m , and the dimensionless couplingg(τ ) := g(tω m )/ω m . This choice will prove convenient throughout the rest of this work, and dimensions can be restored when necessary by multiplying by ω m . This renormalisation effectively is equivalent to the use of time τ and the Hamiltonian To determine the action of this Hamiltonian on an initial pure state, we now proceed to solve the dynamics induced by (2).

Time evolution of the system
We now need an expression for the time evolution operatorÛ (τ ) for a system evolving with (2). The unitary time evolution operator readŝ where ← T is the time-ordering operator [45]. This expression simplifies dramatically when the HamiltonianĤ is independent of time, in which case one simply haŝ U (τ ) = exp −iĤ τ / . As we will here consider stime-dependent light-matter couplingsg(τ ), we instead seek to solve the full dynamics of the time-dependent Hamiltonian.
To do so, we make the ansatz that the time-evolution operator can be written as a product of a number of operatorÛ i . This is possible if there exists a finite set of generatorsĜ j that form a closed Lie algebra under commutation [45]. We thus write the evolution operator (3) aŝ where F j (τ ) are generally time-dependent coefficients determining the influence of the generatorĜ j on the quantum state. Our task is to find these coefficients. This has been done in [36] for an analogous setting. We note that in [36], the operatorsâ andb have been swapped and that in our case the coupling g(t) is preceded by a minus-sign. For clarity, we therefore present the full derivation for the case considered here in Appendix A.
With these techniques, we find that the time-evolution operatorÛ (τ ) can be cast into the convenient form where the operators, given by, form a closed Lie algebra, and where the coefficients are given by Note that in (5), we have transformed into a frame rotating with ω cNa in order to neglect the phase-term e −iΩτ . Given an explicit form ofg(τ ), it is then possible to write down a full solution forÛ (τ ). The decoupling techniques necessary to obtain this compact solution have a long tradition in quantum optics [46] and were generalised and refined recently [45].

Recovering the standard dynamics
Let us show here that this method reproduces the standard evolution operator for the optomechanical Hamiltonian. For this specific case, the light-matter interaction is held constant withg(τ ) =g 0 . The functions (7) simplify to These coefficients allow us to write the time evolution operatorÛ (τ ) aŝ This expression matches that found in the literature (see e.g. equation 3 in [15], which can be obtained with some rearrangement of the terms in (9)).

Initial states of the system
We now proceed to define the initial state of the system that we shall be using throughout this work. We assume that both the optical and the mechanical modes are in a coherent state, which we denote |µ c and |µ m respectively. These states satisfy the relationsâ |µ c = µ c |µ c andb |µ m = µ m |µ m . For the optical field, this is a readily available resource, since coherent states model laser light quite well. The mechanical element in optomechanical systems is most often found in a thermal state or, assuming perfect preliminary cooling, in its ground state, with µ m = 0. An analysis for coherent states can be generalised to thermal states in a straightforward manner by integrating over the coherent state parameter with an appropriate kernel [47].
Throughout the main text of this work, we will assume that the mechanical state is in the ground state with µ m = 0. The initial state |Ψ(0) of the compound system will therefore be |Ψ(0) = |µ c ⊗ |0 .
By starting in a Gaussian state, we ensure that any non-Gaussianity revealed by our work is due to the non-linear coupling in Eq. (2). Indeed, the only way an initially Gaussian state may at any time be non-Gaussian is for the corresponding Hamiltonian to induce some non-linear evolution [48].
Having established the physical model and our initial state, we now proceed to define the measure of non-Gaussianity that we shall use for the remainder of this work.

Measures of deviation from Gaussianity
Given a HamiltonianĤ, and an initial Gaussian stateρ(0), we ask the following question: can we quantify how much the stateρ(τ ) deviates from a Gaussian state at time τ ? This question stems from the following observation. The dynamics of our system is non-linear. Therefore, we expect an initial Gaussian state, characterised by a Gaussian Wigner function, to become a non-Gaussian state at later times. In fact, the only way for a Gaussian state to preserve its Gaussian character would be to evolve through a linear transformation, which is induced by a Hamiltonian with at most quadratic terms in the quadrature operators [49].
To answer our question we first need to find a suitable measure of deviation from Gaussianity. In this work we choose to employ a measure for pure states, which we denote δ, that is based on the comparison between the entropy of the final state and that of a suitably chosen reference Gaussian state [37]. A similar measure has been used to compute features of mixed systems [50].
Let us detail here the construction of the non-Gaussianity quantifier δ(τ ) for our non-linear dynamics. First, our initial stateρ(0) evolves into the stateρ(τ ) at time τ . With our full solutions for the dynamics in Section 2.2, we can find analytic expressions for the first and second moments ofρ(τ ). Then, we construct a stateρ G (τ ), which is the Gaussian state defined by the first and second moments that coincide with those ofρ(τ ). Now, we recall that a Gaussian state is fully defined by its first and second moments. Therefore, if two Gaussian statesρ 1 andρ 2 have equal first and second moments they are the same state [48,49]. In general, our stateρ(τ ) at time τ will not be Gaussian and therefore cannot be specified fully by its first and second moments. This implies that we can introduce a measure δ(τ ) that quantifies howρ(τ ) deviates fromρ G (τ ): where S(ρ) is the von Neumann entropy of a stateρ, defined by S(ρ) := − Tr(ρ lnρ). This measure has been shown to capture the intrinsic non-Gaussianity of the system, and it vanishes if and only ifρ(τ ) is a Gaussian state [37]. In other words, if at all times δ(τ ) = 0 this implies that the state is Gaussian and the dynamics is fully linear. We now note that the time evolution is unitary. This means that S(ρ(τ )) = S(ρ(0)). If we start from a pure state, then S(ρ(0)) = 0 and δ(τ ) = S(ρ G (τ )). Sincê ρ G (τ ) is a Gaussian state, its entropy can be exactly computed using the covariance matrix formalism [48,49]. The covariance matrix consists of the second moments of a quantum state, and can be used to fully characterise a Gaussian state (along with its first moments). This is convenient, as the construction ofρ G involves finding the first of second moments ofρ anyway. While we could compute the entropy forρ G by finding a diagonal basis in the Hilbert space, there exists a straight-forward method within the covariance matrix formalism.
To compute the entropy, we introduce the 4 × 4 covariance matrix σ(τ ) of the stateρ G (τ ). This matrix contains the second moments of the stateρ G (τ ) which in our specific choice of basisX = •} is the anti-commutator, we have defined the expectation value of an operator • := Tr{• ρ G }, and where for the sake of simpler notation, we have chosen not to write out the time-dependence explicitly. To compute the entropy S(ρ G (τ )) we require the symplectic eigenvalues {±ν where the property ν ± (τ ) ≥ 1 holds for all states. The symplectic eigenvalues can be computed by finding the eigenvalues of the object i Ω σ(τ ), where Ω is the symplectic form given by Ω = diag(−i, −i, i, i) in this basis. The von Neumann entropy S(σ) is then given in this formalism by . In summary, the stateρ(τ ) is non-Gaussian at time τ if and only if δ(τ ) > 0.

Measures of deviation from Gaussianity: General behaviour
We shall now infer some general characteristics of the measure of non-Gaussianity δ(τ ). As mentioned above, the symplectic eigenvalues ν ± satisfy ν ± ≥ 1 [49]. Therefore, we can conveniently write ν ± = 1 + δν ± , where δν ± ≥ 0 captures any deviation from purity. If the state is pure at time τ = 0 then it follows that ν ± (0) = 1 [49]. If the evolution is linear, then it is also the case that ν ± (τ ) = 1 for all τ . For closed dynamics, the symplectic eigenvalues may only change if the evolution is non-linear. The statements made here can be generalised to initial mixed states as well. In this case, we would define ν ± = ν 0± + δν ± with ν 0± > 1. Then, we would have that ν ± (0) = ν 0± and, again, linear evolution would imply that ν ± (τ ) = ν 0± . The preceding statements imply that δν ± are functions of the non-linear contributions alone. Thus, when the non-linearity tends to vanish, then δν ± → 0. Among the possible asymptotic regimes we have that δν ± → +∞ or that it becomes constant.
These observations are important. To understand their implications we use the expression ν ± = 1 + δν ± to compute the general deviation from Gaussianity as Using this form, we see that in the nearly linear (Gaussian) regime with only small contributions from the non-linear dynamics, we will have δν ± 1 and therefore On the contrary, in the highly non-linear (non-Gaussian) regime we have δν ± 1 and therefore δ(τ ) ≈ ln δν+ 2 + ln δν− 2 . If the symplectic eigenvalues depend on a large parameter x 1, then one will in general find that they have the asymptotic form ν ± ∼ x N± N± n=0 ν (n) ± x −n for some appropriate real coefficients ν (n) ± , where N ± constitutes the upper limits of the sum [51].
A careful asymptotic expansion of the measure of nonlinearity in this regime gives These general results allow us to anticipate that these general behaviours will be confirmed by the explicit analytical and numerical computations below. A detailed computation can be found in Appendix B.

Results
We now proceed to evolve the initial state |Ψ(0) = |µ c ⊗ |0 with the evolution operatorÛ (τ ) in (5). The explicit form of |Ψ(τ ) for a constant coupling g 0 has previously been used in [15]. To compute the amount of non-Gaussianity of this state δ(τ ), we must first find the elements σ nm of the covariance matrix σ, which we construct form the first and second moments of |Ψ(τ ) . For the initial state defined in (10), the elements of the covariance matrix can be computed with some algebra (see Appendix A.2), where we have transformed into a frame rotating with ΩN a . The elements are given by where we have introduced F := FN aB− + iFN aB+ and θ a := 2(FN 2 a + FN aB+ FN aB− ) for simplicity. Due to the symmetries of the covariance matrix, which in our basis reads σ † = σ, all other matrix elements can be obtained from those listed in (14).
Our challenge now is to compute the symplectic eigenvalues ν ± , given the expressions (14). The process of computing the eigenvalues can be simplified by using the expression 2 ν 2 ± = ∆ ± ∆ 2 − 4 det(σ), which is based on the existence of symplectic invariants [49]. The definition of ∆ is given in Appendix A.3.
The full analytic expression for the symplectic eigenvalues ν ± of σ in (14) is too long and cumbersome to be printed here. It also does not yield any immediate insight into the behaviour of δ(τ ). Instead, we will proceed to derive two analytic expressions for the two different regimes we identified in Section 3.1; small and large coherent state parameters respectively.

Asymptotic behaviour for small coherent state parameter
We begin by looking at the case where |µ c | 2 1. Here, one can take the covariance matrix elements (14) and, after some algebra, show that the perturbative expansion of the symplectic eigenvalues gives This implies that the behaviour of δ(τ ) for small |µ c | goes as in perfect agreement with (12). This approximation suggests that δ(τ ) scales with ∼ |F | 2 |µ c | 2 ln |µ c | to leading order.

Asymptotic behaviour for large coherent state parameter
We now investigate the case where |µ c | 1. Before making any quantitative evaluation, we recall that the measure will have the form (13), where now x ≡ |µ c | 2 . Let us proceed to demonstrate this result analytically for our specific case.
For large µ c , it is clear that whenever θ a = 2πn for integer n, the matrix elements σ 31 , σ 21 and σ 41 vanish, due to the exponentials containing the |µ| 2 factor. Therefore, far (enough) from the times where θ a = 2πn we are left with the following covariance matrix elements and all other elements are zero. Note that we need to keep the next leading order in each element of (17), which is a constant in the case of σ 11 and σ 22 . Naively neglecting of this element would give an incorrect result when computing the entropy, as the neglected factor becomes significant in the logarithm.
With this simplified matrix, we are able to find a simple analytic expression for the symplectic eigenvalues, which reads We note that both eigenvalues grow with |µ c |, as expected. The amount of non-Gaussianity for large µ c is now given by the simple expression which scales asymptotically as δ(τ ) ∼δ(τ ) := 4 ln |µ c |, in perfect agreement with (13). We note that (19) is also valid for time-dependent light-matter couplingg(τ ). In all cases, the nonlinearity grows at with ln |µ c | at leading order.

Application: Constant coupling
Let us now move on to a quantitative analysis of the evolving non-Gaussianity in different contexts. We begin by considering the case where the non-linear lightmatter interaction is constant:g(τ ) =g 0 . To a large extent, this is the case for most experimental systems. The coefficients which determine the time-evolution are those found in (8), and we note that the function F which appears in the covariance matrix elements σ nm is now given by F =g 0 (1 − e −iτ ).
We now proceed to compute the exact measure of non-Gaussianity δ(τ ) for this case. The exact expression is again too long and cumbersome to be reprinted here, but we plot the results in Figure 2 and Figure 3. In Figure 2a we see the measure of non-Gaussianity δ(τ ) as a function of the time τ for different values of the coherent state parameter µ c , over the period 0 < τ < 2 π. It was shown already that the full non-linear dynamics is periodic (or recurrent) with period 2 π, see [18] and this is clearly reflected in our plot. Therefore, we expect, and indeed see, that δ(0) = δ(2 π). Since δ(τ ) is a continuous function over a closed interval, it must take a maximum value in between the two points 0 and 2 π. It turns out that the maximum is achieved at τ = π, which will be the time of reference when discussing the magnitude and scaling of the effects withg 0 and µ c . Increasing µ c yields a logarithmic increase in δ(τ ), which is evident from the approximation in Eq. (19). Figure 2a also implies that for closed dynamics, the non-linear system will almost immediately become maximally non-Gaussian. It will then retain approximately the same amount of non-Gaussianity until τ = 2π, meaning there will be a rapid decrease of non-Gaussianity before the system revives again. Such plateaus will be justified mathematically by analysing the asymptotic behaviour for large |µ c |. This shows that the maximum amount of non-Gaussianity available during one cycle can be accessed almost immediately without requiring the system to evolve for a long time. As a side remark, we note that the functional form of δ(τ ) in Figure 2a resembles the linear entropy of the traced-out subsystems as found in [15,18].
To better understand the behaviour of δ(τ ) for small times τ , we plot the behaviour of log 10 δ(τ ) for τ 1 for different values of µ c in Figure 2b. We note that δ(τ ) increases quickly at first, but soon tends to a near-constant value. This means that δ(τ ) grows linearly for an interval of small times, which can be seen as the increasing and decreasing parts in Figure 2a.
We now proceed to examine the scaling behaviour of δ(τ ) at constant time τ = π. Figure 3a shows a log 10 -log 10 plot of the measure δ(τ ) as a function of the non-linear couplingg 0 . Asg 0 increases, the amount of non-Gaussianity first grows linearly in the logarithm, then plateaus asg 0 increases further. The same behaviour occurs for larger µ c , only more rapidly. This suggests that it becomes increasingly difficult to add more non-Gaussianity to the system by increasingg 0 , meaning that focusing on increasing the couplingg 0 only gives marginal returns. Similarly, 3b shows log 10 δ(τ ) as a function of increasing µ c for various values ofg 0 .

Small coherent state parameters in the constant coupling regime
For a small amplitude coherent state, with |µ c | 2 1, we found in 16 that δ(τ ) scales with ∼ |F | 2 |µ c | 2 ln |µ c |. Given the explicit form of F , we see that it scales with F ∝g 0 . Since δ(τ ) in this regime is proportional to |F | 2 , it follows that δ(τ ) grows quadratically with the light-matter coupling in this regime.
To see how well the asymptotic form in (19) approximates the exact measure, we have plotted both the exact form of δ(τ ) (solid lines) with the asymptotic form (dashed lines) in Figure 4. We note that even for |µ c | ∼ 1, the asymptotic form closely approximates the exact value of δ(τ ) when far enough away from τ = 0 and τ = 2πn. In fact, it becomes even more accurate as the coherent state parameter µ c increases.
The asymptotic form also becomes even more accurate once we also increaseg 0 , as evident in Figure 4b. Forg 0 = 10 2 , the approximation is almost entirely accurate. This occurs because the function θ a increases withg 0 , which further suppresses the offdiagonal covariance matrix elements at the beginning and end of each cycle. Finally, we note that the approximation breaks down close to the recurrnt points, as we have discussed above.

Applications: Time-dependent coupling
In all physical systems, such as optomechanical cavities, the confining trap is not ideal. This means that, in general, the couplingg(τ ) is time-dependent as a consequence of, for example, trap instabilities. Here, time-dependent variations such as phase fluctuations in the laser beam used to traps a levitated bead will modulate the coupling. In this work, we want to exploit the possibility of controlling the coupling g(τ ) by considering its periodic modulation in time. In practice, such time-dependent control would be achievable for an optically trapped and levitated dielectric bead that interacts with a cavity field, by controlling the optical phase of the trapping laser field. In fact, such phase determines the bead's equilibrium position which, in turn, affects the cavity light-bead coupling through the varying overlap between the bead and the cavity mode function. Technically, the trapping laser's optical phase may be controlled through an acoustic-optical modulator. Alternately, control on a bead's equilibrium position may also be enacted by adopting Paul traps, which work for levitated nanospheres [52], and have been used to shuttle ions across large distances, typically for the purpose of quantum information processing [53,54].

Modelling the trap instability
We shall model this situation by assuming the coupling has the form g(τ ) =g 0 (1 + sin(Ω 0 τ )) .
Here,g 0 is the expected value of the coupling, is the amplitude of oscillation and Ω 0 := ω 0 /ω m is the dimensionless frequency that determines how the coupling oscillates in time. We can insert this ansatz in the general expressions (7) and obtain an explicit form for this case. The full expressions for the coefficients in (7) are again very long and cumbersome, and we do not print them here. They are listed in (C.1).
We can now compute δ(τ ) for this time-dependent coupling and we display the results in Figure 5. In 5a we plot δ(τ ) vs. τ for different values of the oscillation frequency Ω 0 . Note that we here include a larger range of τ to capture potentially recurring behaviour. In the limit Ω 0 → 0 we recover the time-independent solution, as expected. Interestingly, when Ω 0 = 0 we see that we can achieve higher values for the nonlinear measure δ(τ ). This is especially pronounced as Ω 0 → 1, where the trap oscillation frequency is equal to the mechanical frequency ω m , for which δ(τ ) ceases to oscillate periodically, but instead steadily increases. We discuss this case in detail in the following section.

Resonance
The functions (C.1) contain denominators of the form Ω 0 − 1. Therefore, among all possible values of Ω 0 , we can ask what happens on resonance, that is, when Ω 0 = 1. We have already evidence from Figure 5a that the system should behave markedly differently.
We have plotted in Figure 5 the exact measure of non-Gaussianity δ(τ ) in the resonant case for different values of µ c . As anticipated before, in this case we no longer have periodic dynamics. Instead, the non-linearity increases as ln τ . Formally, this growth can continue for arbitrarily large times τ , however, the maximum time τ that can be achieved in practice is limited by the coherence time of the experiment. Similarly, we plotted δ(τ ) for various values of in Figure 5c. We note that δ(τ ) oscillates increasingly rapidly with larger but with decreasing amplitude for increasing τ , as |F | 2 ∼g 2 0 2 τ 2 becomes the dominant term for τ 1. Finally, in Figure 6, we compare the exact measure δ(τ ) at resonance with the asymptotic form derived in (19). The solid lines represent the exact measure δ(τ ) and the dashed lines represent the asymptotic expression. In Figure 6a we compare them for different values of µ c . We note that, except for at very small τ , the asymptotic form is entirely accurate and gets even more precise for increasing values of µ c . This is a consequence, as we noted before, of the exponentials in (14) that suppress some elements for large µ c , unless θ a = n π. Similarly, in Figure 6b we have plotted δ(τ ) and its asymptotic form for different values of the oscillation amplitude . Again, the suppression of the exponentials with increasing means that larger values of yield a more accurate expression.

Discussion and practical implementations
We have employed a measure of non-Gaussianity δ(τ ) in order to quantify the deviation from linearity of an initial Gaussian state induced by the Hamiltonian (1). Our results show that, for a constant light-matter couplingg 0 , the non-Gaussianity δ scales differently in two contrasting regimes: (i) For a weak input coherent state |µ c |, the non-linear character of the state grows asg 2 0 |µ c | 2 ln |µ c |, (ii) conversely, for large |µ c |, the non-linear character of the state grows logarithmically with the quantitỹ g 0 |µ c |. The same general scaling with |µ c | occurs wheng(τ ) is time-dependent.
Crucially, we also find that the amount of non-Gaussianity can be enhanced by driving the light-matter coupling at mechanical resonance. We will now discuss these results in the context of concrete experimental setups, and specifically discuss how the modulated light-matter coupling can be engineered.

Modulating the light-matter coupling in physical systems
We saw in Section 6 that the amount of non-Gaussianity in the system increases when the light-matter couplingg(τ ) is modulated. How can this modulation be engineered in the laboratory? We shall explore several methods in this section. As a basis for this discussion, we present a derivation of a time-dependent lightmatter coupling for levitated nanobeads in Appendix D, which is based on the work in [55]. There are several practical ways in which one may envisage to increase the non-linearity by modulating the coupling, depending on the nature of the trap at hand: i) Optically-trapped levitating particles. The effect that we are looking for can be realised by modulating the phase of the trapping laser beam (which, in turn, can be achieved through an acousto-optical modulator). In our derivation in Appendix D, this phase is denoted by ϕ(τ ) and it affects the light-matter coupling strength by determining the particle's location with respect to the standing wave of the cavity field. Thus if we let ϕ(τ ) = π 2 (1 + sin Ω 0 τ ), with Ω 0 = ω 0 /ω m , and where ω 0 is the phase modulation frequency, we obtain the expression used in Section 6. If, then, the phase frequency is resonant with Ω 0 = 1, it should be possible to increase the non-Gaussianity even further.
ii) Paul traps. The shuttling of ions has been demonstrated [54,53] using Paul traps, which are customarily used for ions but which have also recently been used for trapping nanoparticles [52,56,57]. These works indicate that a modulation of the particle's position, and hence, a modulation of the coupling as per point i), can be obtained in a Paul trap as well.
iii) Micromotion in hybrid traps. Paul traps display three different kinds of particle motion. Firstly, we have thermal motion, whereby the particle moves around the trap. Secondly, and most importantly to our scheme, we have micromotion, which induces small movements around the potential minimum.
Finally, there is mechanical motion, which is the harmonic motion in the trap, here denoted by ω m . Since the micromotion moves the bead around the potential minimum with a frequency ω d , this already modulates the light-matter coupling, and is, in a way, an equivalent implementation to the "shaking" of the trap. If the micromotion can be engineered to occur with a frequency ω d equal to ω m , then one could, instead of averaging it out, adopt the micromotion's variables to increase the non-Gaussianity with the scheme we propose in Section 6.2.
There are potentially many more ways in which the light-matter coupling could be modulated, including with optomechanically induced transparecny [58,59] and by using the Kerr effect to change the refractive index of the oscillator.
We conclude that the enhancement of the non-linearity predicted by our work can be realised in experiments, given the capabilities mentioned above. There are, of course, many challenges to be overcome. In fact, to take advantage of the rather slow logarithmic scaling with time τ , one must keep the system coherent for longer, which is difficult. However, although our analytical results are restricted to Hamiltonian systems, we note that there is no reason to expect that this enhancement should disappear in a noisy setting.

Experimental regimes
There are two relevant experimental regimes for optomechanical systems. They are determined by the magnitude of the light-matter coupling g compared to the other frequencies in the system. In the weak single-photon optomechnical coupling regime, the light-matter coupling g is small compared to the resonant frequency ω m and the optical decoherence rate κ. Such experiments usually involve a strong laser drive, which tends to wash out the non-linearity. In the strong single-photon coupling regime, non-linear effects are in practice small but more significant. Under these conditions, a single photon displaces the mechanical oscillator by more than its zeropoint uncertainty and weak optical fields tend to be used [60]. In summary, most approaches fall into one of two categories: (i) small g and linearised dynamics and (ii) large g and low number of photons.
Ultimately, our work suggests that we can further increase the amount of non-Gaussianity by modulating the light-matter coupling. We emphasize that this scheme is applicable in both the weak and strong coupling regimes. This sets it apart from other schemes, which usually focus on enhancing the non-Gaussianity in one of the two categories mentioned above.

Conclusions
We have quantified the non-Gaussianity of initially Gaussian coherent states evolving under the standard, time-dependent optomechanical Hamiltonian. We used a measure of non-Gaussianity based on the relative entropy of a state and characterised the deviation from Gaussianity of the system. Our techniques allowed us to derive asymptotic expressions for small and large optical coherent-state amplitudes, see Equation (16) and Equation (19) respectively. We found that for coherent states with amplitude |µ c | ≥ 1, the amount of non-Gaussianity grows logarithmically with the input average number of excitations |µ c | and with the light-matter coupling. At resonance, this is further enhanced by a logarithmic scaling with time of interaction.
An important and promising aspect of our study consists in the prediction that the amount of non-Gaussianity in the system can be increased by devising a timemodulated optomechanical coupling. This is particularly pronounced at resonance, where the driving frequency matches the inherent mechanical frequency. Even though we have not discussed direct applications of the non-linear dynamics, the measure of non-Gaussianity we have adopted, based on the relative entropy, has strong operational implications [8]. As such, this points to a new, and practically accessible, mechanism to enhance the non-linear character of optomechanical dynamics at a given light-matter coupling strength. While we have not considered open systems dynamics, one should not expect significant qualitative changes to the non-Gaussianity in the presence of noise.
Our work can be extended to more complicated Hamiltonians of bosonic modes, and we can include modifications such as squeezing of the mechanical state. This setting will be explored in future work.

Appendix A. Derivations and tools
In this appendix, we will derive the coefficients in (7) that determine the time-evolution of the system. This follows the derivation in [36]. We will also show the explicit timedependence of the second moments and discuss some methods related to computing the symplectic eigenvalues of the covariance matrix.
Appendix A.1. Properties of of the non-linear Hamiltonian Finding a simple expression forÛ (t) is straight-forward when the light-matter coupling g is not time-dependent. If g(t) is time-dependent we require a more rigorous framework. This is what we present here.
We will here follow the derivation in Appendix A in [36]. We note that compared with [36], we have here swapped the definition ofâ andb, and we have a minus-sign in front of g(t). In this derivation, we will also use t rather than τ to show the full dependence of ω m .
For the time-dependent HamiltonianĤ in (1), the time-evolution operator is given byÛ where ← T is the time-ordering operator. The basis for decoupling the operator is finding a Lie algebra of generatorsĜ i that induce the time-evolution. This Lie algebra must be closed under commutation, that is, either [Ĝ j ,Ĝ k ] ∝Ĝ l , or [Ĝ j ,Ĝ k ] = c where c is a scalar. This will allow for the terms inÛ (t) to be moved with the Baker-Campbell-Hausdorff formula such that U (t) can be written in a simpler form.
We start with the ansatz that the evolution operator (3) can be written aŝ where F j are coefficients corresponding to each of the generatorsĜ j . Our task is now to find the coefficients F j . We begin by defining the operatorsĜ j in the algebra: It can be verified that the operators in (A.3) form a closed Lie algebra under commutation. With these operators, our ansatz can be written aŝ To find the coefficients, we note the following equivalence: Differentiating both sides brings down the Hamiltonian (1) on the left, which we here write in terms of the generators (6). We then multiply both sides by U † to obtain the following differential equation: whereḞ i = d dt F i . This is the equation that determines the coefficients. We can now commute all the operators through, where we find  Depending on the form of g(t), we can now use these equations to find a simplified form ofÛ (t).

Appendix A.2. Time-evolution of first and second moments
In the Heisenberg picture, the time evolution of the mode operatorsâ andb induced by the Hamiltonian considered here is simplyâ(t) :=Û † (t)âÛ (t) andb(t) :=Û † (t)bÛ (t).
In terms of the generators of the Lie algebra defined in (A.3), we explicitly havê This expression can also be rewritten in more compact notation aŝ where we have introduces the quantities For the initial coherent state |Ψ(t = 0) = |µ c ⊗ |µ m (note that in the main text we specialise to µ m = 0), and ignoring the global phases e −i ωc t , which can be done by transforming into a frame rotating with ω câ †â , we obtain where we have introduced F := FN aB− + iFN aB+ and θ a := 2(FN 2 a + FN aB+ FN aB− ). Note that, in the main text, we had set t = τ /ω m and µ m = 0, but we here show their explicit dependence.
These expectation values can now be used to compute the elements of the covariance matrix σ, which in our basis are given by where we have suppressed the time-dependence for notational convenience. The covariance matrix elements are, with ω m and µ m reinstated as compared with (14) given by In the first case, δν ± 1, and we have and it leads us to the expressions: 4 sin(τ ) cos(Ω 0 τ )(cos(τ ) cos(Ω 0 τ ) − 2) + 8 cos(τ ) sin(Ω 0 τ ) It can be seen from these expressions that there are some resonances expected, namely a drastic change in the behaviour of (some of) the functions in the limit Ω 0 → 1, which occurs when ω 0 = ω m . It is straight-forward to see how the terms FN aB+ and FN aB− simplify as Ω 0 → 0 by noting that lim Ω0→1 sin τ (1 − Ω 0 /(1 − Ω 2 0 ) = τ /2. The long expression for FN 2 a is more challenging. We note that the terms independent of remain unchanged with Ω 0 .
Thus, at resonance, these coefficients read: Appendix D. Deriving the modulated light-matter coupling In this appendix, we will show how the time-dependent term used in Section 6 can be derived for levitated nanobead systems. In [55], a fully general theory of light-matter coupling is presented. We will recount some of the derivation here and show how the cavity volume can be modulated in a manner such that it is useful to our scheme. Given a number of assumptions regarding the light-matter interaction (see [55] for a full description) the full Hamiltonian that describes the light-matter interaction for a homogeneous dielectric object is the following: The termĤ f m =p 2 /2M , where M is the total mass of the system, is the kinetic energy of the centre-of-mass position along the cavity axis.Ĥ F c = ω câ †â is the energy of the cavity mode.Ĥ f out andĤ f free are terms describing an open system, which we shall ignore in this work. We likewise ignoreĤ i cav−out which describes a coupling between the cavity input and the output mode.
The last termĤ i diel describes the light-matter coupling and can be written in the general formĤ where P (x) is the polarization of the levitated objects (which we assume to be a scalar quantity) andÊ(x) is the total electric field, which can be obtained from solving Maxwell's equations given a set of well-defined boundary conditions. The quantised modes of the electric field can thus be written as [49] where s is the spin-polarization index and m signifies the field-mode number, and E m = ωm 2 0 Vc is the field amplitude with V c being the cavity mode volume. The functions χ s,m must obey the spatial solutions to the wave-equations, where the full classical solutions separate into E(r, t) = χ(r) T (t).
If we assume that the polarization is given by P (x) = c 0 E(x), we obtain the simpler expressionĤ where c = 3 r −1 r +2 , and where r is the relative dielectric constant of the nanodiamond. We now assume that the electric field operators are displaced by a classical part: a → â 0 +â. The classical part â 0 will form the optical trapping field, while the quantum part describes the light-matter interaction.
Thus the classical contribution to the electrical field is given by where α is a complex prefactor and f (x) is a complex function which describe the standing waves inside the cavity. We now write our full electric field asÊ tot (x) = E(x) + E(x), whereÊ(x) is the quantum contribution containingâ andâ † , and E(x) is the classical part. The full Hamiltonian is noŵ The classical contribution, E(x) will yield a trapping frequency, while the operator termsÊ(x) will yield the light-matter interaction term for the levitated sphere. The cross-term,Ê(x)E 2 (x) will generate elastic scattering processes inside the cavity which converts cavity photons and tweezer photons into free modes [55]. We shall ignore them here and focus on the generation of the trapping frequency ω m and the coupling g(t).
We begin with the trapping frequency.
Appendix D.1. Mechanical trapping frequency We now assume that the classical field has a Gaussian profile which extends in the y-direction for a cylindrical geometry. The cavity extends along the z-direction. We here follow the derivation presented in [61]. If we denote the radius of the cylinder by r, we can write down the trapping field as where E 0 = Pt 0cπW 2 0 , P t is the trapping laser power and W 0 is the beam waist with the full beam as a funtion of y being W (y) = W 0 1 + y 2 λ 2 π 2 W 4 0 . It follows that the narrowest part of the beam W 0 occurs at y = 0, which is the minimum in the potential where the nanobead is trapped.
We can now expand [E(y, r)] 2 to second order in r and y around the origin y 0 = r 0 = 0. We start with the exponential, which we expand as . (D.8) Next, we expand the inverse beam width to second order in y: (D.9) Combining the two expressions give us [E(y, r)] 2 ≈ E 2 0 1 − (D.10) If we now assume that y W 0 , meaning that the beam waist is much larger than the region we consider, we can approximate the above as We then insert this now constant expression into the integral for the Hamiltonian and we drop all constant terms as they are just constant energy shifts. To perform this integral, we now assume that the radius R of the bead is much smaller than the wavelength of the light. This is often referred to as the 'point-particle approximation', or the Rayleigh approximation. Essentially, this means that the field inside the bead is constant (although the field still changes in space with x and y). Thus we can assume that wherever the sphere is located in the field, the integral just simplifies to the volume of the sphere times the field amplitude. For a derivation which includes arbitrary particle sizes, see [62]. This gives where V is the integration volume. The result is a harmonic trapping of the form where ρ = m V is the density of the levitated object and where we have used E 2 0 = 2I c 0 , where I is the intensity of the laser beam, and c = 3 r −1 r +2 .

Appendix D.2. The light-matter interaction term
We now come to the most important term, which is the light-matter interaction term denoted g in this work. We will continue to follow the derivation in [55] to show exactly where time-dependence could potentially be included. If the sphere is sufficiently small, we can choose a TEM 00 (transverse electromagnetic mode) as the cavity mode, which is aligned in the z-direction. In this mode, the cross-section in x and y is perfectly Gaussian, and it is one of the most commonly used modes in experiments. If the sphere is smaller than the laser waist and if it is placed close to the centre of the cavity, we can approximate the field at the centre of the beam by Here, the laser waist is given by W c = λL (2π) 2 , L is the cavity length. λ is the laser wavelength. We assume that the wave-vector k c points in the z-direction, along the axis of the cavity, and ϕ is a generic phase which determines the minimum of the potential seen by the bead. For laser-trapped nanobeads, this phase can be made time-dependent, whereas for a Paul trap, it is static. We will leave out the timedependence for now for notational simplicity. Finally,â andâ † are the annihilation and creation operators of the electromagnetic field.
To obtain the Hamiltonian term, we now integrate over the full energy within the volume of the nanobead. For a bead situated at r = (x, y, z) leads tô 16) In all traps, optical and Paul traps, the bead is trapped in a minimum of the potential. This occurs at ϕ = π 2 . In optical traps, we can now modulate ϕ → ϕ(t), to change the light-matter coupling. If we let ϕ(t) = π 2 (1 + sin ω 0 t), we obtain the scenario we investigate in Section 6. Finally, we note that there might be many additional ways in which the coupling can be modulated that we have not discussed in this work.