Enhanced continuous generation of non-Gaussianity through optomechanical modulation

We study the non-Gaussian character of quantum optomechanical systems evolving under the fully nonlinear 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 both a constant and time-dependent cubic light–matter coupling and study its general and asymptotic behaviour. We find analytical approximate expressions for the measure of non-Gaussianity and show that initial thermal phonon occupation of the mechanical element does not significantly impact the non-Gaussianity. More importantly, we also show that it is possible to continuously increase the amount of non-Gassuianity of the state by driving the light–matter coupling at the frequency of mechanical resonance, suggesting a viable mechanism for increasing the non-Gaussianity of optomechanical systems even in the presence of noise.


Introduction
Understanding nonlinear, interacting physical systems is paramount across many areas in physics. Specifically, 'nonlinear' (or 'anharmonic') dynamical systems include all those whose Hamiltonian cannot be expressed as a second-order polynomial in the quadrature operators. Crucially, these systems allow us to generate non-Gaussian states, which cannot be done, given only quadratic couplings. One family of system where this is possible are optomechancial systems, where light interacts with a mechanical element through a cubic interaction term.
In recent years, the intrinsic value of nonlinear systems, as opposed to the aforementioned limitations that linear systems face, has been made clearer and more rigorous. It has been shown that nonlinearities 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 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 [13][14][15][16][17][18]. 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 that the system does not have to be cooled to the ground state in order to access significant amounts of the non-Gaussianity. Finally, we discuss several methods by which this modulation can be realised for optomechanical systems 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 introduce the measure of non-Gaussianity from [50], and in section 4 we derive some generic results based on the measure which apply in different regimes. We then proceed to examine the behaviour of the non-Gaussianity in optomechanical systems 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 continuously generated non-Gaussianity. Both preceding sections also include an analysis of the open system dynamics. Finally, in section 7 we discuss our results and propose various methods by which the modulation of the optomechanical coupling can be achieved. Section 8 concludes this work, and many of the derivations used in this work can be found in appendices A, B, C, and E.

Dynamics
We start our work by presenting the necessary mathematical tools needed to solve the dynamics. We refer the reader to appendix B for a more extensive introduction on the techniques presented below.

Hamiltonian
We begin by considering two bosonic modes corresponding to an electromagnetic mode and a mechanical oscillator. Without loss of generality, we shall henceforth refer to the electromagnetic mode as the optical mode. The operators â and b of the cavity and mechanical modes respectively, obey the canonical commutation relations a a b b , , while all other commutators vanish. The radiation pressure induces a nonlinear interaction between the light and mechanics, and the whole systems is modelled by the following Hamiltonian: where ω c and ω m are the frequencies of the cavity mode and the mechanical mode respectively, and g(t) drives the, potentially time-dependent, nonlinear light-matter coupling. The 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 [51].
To simplify our notation and expressions, we rescale the laboratory time t by the frequency ω m , therefore introducing the dimensionless time t To determine the action of this Hamiltonian on initial states, 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 U t ( ) for a system evolving with (2). The unitary time-evolution operator reads where  ¬ is the time-ordering operator [52]. This expression simplifies dramatically when the Hamiltonian Ĥ is independent of time, in which case one simply has U H exp i t t = -( ) (ˆ). As we will here consider timedependent light-matter couplings g t( ), 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 operators U ĵ . This is possible if there exists a finite set of generators G ĵ that form a closed Lie algebra under commutation [52]. We thus write the evolution operator (3) as here F j (τ) are generally time-dependent coefficients determining the influence of the generator G ĵ on the quantum state. Our task is to find these coefficients, which has been done in [49] for an analogous setting. We note that compared with [49], the operators â and b 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 B.
With these techniques, we find that the time-evolution operator U t ( ) can be cast into the convenient form U t e e e e , 5 here the operators, given by, ote that in (5), we have transformed into a frame rotating with N â W in order to neglect the phase-term e i t -W .
Given an explicit form of g t( ), it is then possible to write down a full solution for U t ( ). The decoupling techniques necessary to obtain this compact solution have a long tradition in quantum optics [53] and were generalised and refined recently [52]. Finally, before we proceed, we also define the following two parameters: hese quantities will be useful when discussing features of the non-Gaussianity.

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 with g g 0 t =( )˜. The functions (7) simplify to his expression matches that found in the literature (see e.g. equation (3) in [22], which can be obtained with some rearrangement of the terms in (10)).

Initial states of the system
In this work, we will examine the non-Gaussianity of the evolved state given two initial states: a coherent state and a thermal coherent state. . 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. The initial state 0 Y ñ | ( ) of the compound system will therefore be The assumption that the mechanics is in the ground state is not always justified, and therefore we shall also consider the non-Gaussianity of cases where the mechanics is in a thermal coherent state [54]. Such a state is obtained simply by integrating over the coherent state parameter with an appropriate kernel [55]. We define the thermal state th r as here n is the average thermal phonon occupation of the state, bñ | is a coherent state, and the integration occurs over the full complex space. Assuming that the optical mode in the coherent state c m ñ | , the full initial state is therefore given by By starting in an initial Gaussian state, we ensure that any non-Gaussianity revealed by our work is due to the nonlinear coupling in equation (2). Indeed, the only way an initially Gaussian state may at any time be non-Gaussian is for the corresponding Hamiltonian to induce some nonlinear evolution [56]. We do however note that while the measure of non-Gaussianity that we shall make use of has a clear and operational notion of the measure for pure states, it is harder to make statements about the non-Gaussianity of states that are mixed, such as the coherent thermal state th r . See section 7 for a discussion of the properties of the relative entropy measure.

Open system dynamics
All realistic systems experience decoherence. In optomecahnical systems, this manifests as photons leaking from the cavity or as damping of the mechanical motion. Given a sufficiently weakly coupled environment, we can model the open dynamics of the system with the help of the Lindblad equation [57]. We note, however, that there is increasing evidence that the standard master equation treatment breaks down, especially in the strong coupling regime [58]. Here, we shall only consider weak coupling in the presence of noise, and so the Lindblad equation is given by where L ĉ and L m are the Lindblad operators for the optics and mechanics, respectively. To model photon and phonon decay, we assume that L a c c k is the optical decoherence rate and m k is the phonon decoherence rate.
While analytic solutions for this particular choice of L m were obtained in [22], photon decay can currently only be modelled numerically. We therefore make use of the Python library QuTiP to simulate the noisy state evolution and its effect on the non-Gaussianity of the resulting state. We shall examine the non-Gaussianity for open systems in sections 5 and 6, but first, we define the measure of non-Gaussianity.

Measures of deviation from Gaussianity
Given a Hamiltonian Ĥ , and an initial Gaussian state 0 r ( ), we ask the following question: can we quantify how much the state r t ( ) deviates from a Gaussian state at time τ? This question stems from the following observation. The dynamics of our system is nonlinear. 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 [59].
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 [50]. A similar measure has been used to compute features of mixed systems [60].

Measures of deviation from Gaussianity: definition
Let us detail here the construction of the non-Gaussianity quantifier δ(τ) for our nonlinear dynamics. First, our initial state 0 r ( ) evolves into the state r t ( ) 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 r t ( ). Then, we construct a state G r t ( ), which is the Gaussian state defined by the first and second moments that coincide with those of r t ( ). Now, we recall that a Gaussian state is fully defined by its first and second moments. Therefore, if two Gaussian states 1 r and 2 r have equal first and second moments they are the same state [56,59]. In general, our state r t ( ) 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 r t ( ) deviates from G r t ( ): where S r (ˆ) is the von Neumann entropy of a state r , defined by S Tr ln r r r -(ˆ) ≔ (ˆˆ). This measure has been shown to capture the intrinsic non-Gaussianity of the system, and it vanishes if and only if r t ( ) is a Gaussian state [50]. 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 r t r = (ˆ( )) (ˆ( )). If we start from a pure state, then S 0 0 r = (ˆ( )) and S G d t r t = ( ) (ˆ( )). We discuss the case where the initial state is mixed in section 7.
3.2. Measures of deviation from Gaussianity: computation using the covariance matrix formalism Since G r t ( ) is a Gaussian state, its entropy can be exactly computed using the covariance matrix formalism [56,59]. 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 r involves finding the first of second moments of r anyway. While we could compute the entropy for G r 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 s t ( ) of the state G r t ( ). This matrix contains the second moments of the state G r t ( ) which in our specific choice of basis In summary, the state r t ( ) 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 δ(τ). The following analysis only holds when the system is pure, which in our case means that we assume that both the optics and mechanics start off as coherent states. As mentioned above, the symplectic eigenvalues ν ± satisfy ν ± 1 [59]. Therefore, we can conveniently write ν ± =1+δν ± , where 0  dn  captures any deviation from purity. If the state is pure at time τ=0 then it follows that ν ± (0)=1 [59]. 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 nonlinear. In this case, we would define Then, we would have that ν ± (0)=ν 0± and, again, linear evolution would imply that ν ± (τ)=ν 0± . The preceding statements imply that δν ± are functions of the nonlinear contributions alone. Thus, when the nonlinearity tends to vanish, then 0 dn   . Among the possible asymptotic regimes we have that dn  +¥  or that it becomes constant. These observations are important. To understand their implications we use the expression 1 n dn = +   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 nonlinear dynamics, we will have 1 dn   and therefore On the contrary, in the highly nonlinear (non-Gaussian) regime we have 1 dn   and therefore 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 D.

General results
We now proceed to evolve the initial state 0 r ( ) with the evolution operator U t ( ) in (5). To compute the amount of non-Gaussianity of this state δ(τ), we must first find the elements σ nm of the covariance matrix s, which we construct form the first and second moments of r t ( ). We do so in full generality, meaning that the light-matter coupling g t( ) can take the form of any time-dependent function.
We have computed the second moments and the covariance matrix s in appendix C. The second moments for the mechanical being in an initial coherent state and in an initial coherent thermal state can be found in (B.5) and (B.9) respectively. These can then be used to compute the covariance matrix elements σ nm where n, m take values 0, 1, 2, 3. We have explicitly computed the elements of s for a mechanical coherent state in (B.6).
Our challenge now is to compute the symplectic eigenvalues ν ± , given the expressions (B.6). The process of computing the eigenvalues can be simplified by using the expression 2 4 d e t 2 2 s n = D  D - ( ), which is based on the existence of symplectic invariants [59]. The definition of Δ is given in appendix B.2. The full analytic expression for the symplectic eigenvalues ν ± of s in (B.6) 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.3; small and large coherent state parameters respectively.

Asymptotic behaviour for a small optical coherent state parameter
We begin by looking at the case where 1 c 2 m  | | and where the mechanics is initially in a coherent state. Here, one can take the covariance matrix elements (B.6) and, after some algebra, show that the perturbative expansion of the symplectic eigenvalues gives This implies that the behaviour of δ(τ) for small c m | | goes as in perfect agreement with (16). This approximation suggests that δ(τ) scales with F ln These expressions do not hold if the mechanical element is initially mixed. However, we will find in the next section that initial phonon occupation only marginally affects the non-Gaussianity.

Asymptotic behaviour for a large optical coherent state parameter
We now investigate the case where 1 c m  | | . Our goal is to derive an analytic expression for the non-Gaussianity that can be used to analyse the overall features of δ(τ). Before making any quantitative evaluation, we recall that the measure will have the form (17), where now x c 2 m º | | . Let us proceed to demonstrate this result analytically for this specific case.
For large μ c and for the mechanics in the ground-state μ m =0, it is clear that whenever n 2 a q p ¹ for integer n, the matrix elements σ 31 , σ 21 and σ 41 in (B.6) vanish, due to the exponentials containing the factor c 2 m | | . 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. We have kept the full expression for σ 11 because it reproduces some key elements of d t ( ), which we shall discuss later. Therefore, we do not expect the thermal occupation of the mechanics to significantly affect the non-Gaussianity that can be accessed in this system. Note also that we need to keep the next leading order in each element of (20), 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 [61]. If the thermal element is in a coherent thermal state, however, the expectation values of a á ñ changes slightly and 11 s in (20) will look different. However, if we approximate 11 s as 1 2 11 c 2 s m » + | | , which follows from that a 0 á ñ~for very large μ c , the non-zero covariance matrix elements of the coherent thermal mechanical state are the same as in (20). We therefore conclude that an initially thermal coherent mechanical state will also exhibit most of the non-Gaussianity we will examine for coherent mechanical states.
With this simplified matrix, we are able to find a simple and analytic expression for the symplectic eigenvalues, which reads We note that both eigenvalues grow with c m | |, as expected from our analysis in section 3.3. The amount of non-Gaussianity for large μ c is now given by the following expression (17). Note that (22) is also valid for a time-dependent light-matter coupling g t( ). In all cases, the nonlinearity grows as ln c m | | to leading order. In sections 5.2 and 6.3 we will compare the asymptotic measure d t( ) with the full measure δ for different cases.

Applications: 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 nonlinear light-matter interaction is constant: To a large extent, this is the case for most experimental systems. The coefficients which determine the time-evolution are those found in (9), and we note that the function F, defined in (8), which appears in the covariance matrix elements σ nm is now given by F g 1 e ). We now proceed to compute the exact measure of non-Gaussianity δ(τ) for constant coupling g 0 and with the system initially in two coherent states. The exact expression is again too long and cumbersome to be reprinted here, but we plot the results in figures 2 and 3. In figure 2(a) we plot the measure of non-Gaussianity δ(τ) as a function of time τ for different values of the coherent state parameter μ c , over the period 0<τ<2 π. The other parameters are set to g 1 0 = and μ m =0. It is known that the full nonlinear dynamics is periodic (or recurrent) with period 2 π, see [25], whenever g 0 2 is an integer and this is clearly reflected in our plot. At τ=2π, the optics and mechanics are no longer entangled, and while the mechanics returns to its initial coherent or coherent thermal state (see supplemental note 1 in [25] for an explicit proof), the final optical state will depends on the value of g 0 . For example, when g 0.5 0 = , the cavity state becomes a superposition of coherent states at τ=2π, also known as a cat state [22]. However, if g 0 2 is integer, we obtain a phase factor of e 1 g 2 i 0 2 = p˜i n the optical state, and the optics returns to an initial state as well. This is the case in figure 2, where δ(2π)=0. We will make use of the asymptotic measure defined in section 4.2 to analyse this behaviour, see section 5.2. Furthermore, while it might seem that the non-Gaussianity peaks at τ=π, the measure d t ( ) exhibits a local minimum which grows increasingly narrow with larger μ c . This is apparent from figure 2(b) where we have shown a close-up of d t ( ) around τ=π for increasing values of μ c , and for g 1 0 = and μ m =0. The dip occurs because at τ=π, we find that g 2 a 0 2 q p = -˜and F g 2 0 = -˜. Thus, for integer g 0 2 , we have sin 2 0 a 2 q = and σ 11 becomes ) . The non-zero exponent causes the non-Gaussianity to temporarily decrease, and the same behaviour occurs in the other covariance matrix elements, resulting in the dip.
As already noted, increasing μ c yields a logarithmic increase in δ(τ), which is evident from the approximation in equation (22). Figure 2(a) also implies that for closed dynamics, the nonlinear 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. The appearance of these plateaus 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 2(a) closely resembles the linear entropy of the traced-out subsystems as found in [22,25].
To better understand the behaviour of δ(τ) for small times τ, we plot the behaviour of log 10 d t ( ) for τ = 1 for different values of μ c in figure 2(c). We note that δ(τ) increases quickly at first, but soon tends to a nearconstant value. This means that δ(τ) grows linearly for an interval of small times, which can be seen as the increasing and decreasing parts in figure 2(a).
Finally, we proceed to examine the scaling behaviour of δ(τ) at fixed time τ=π. Figure 3(a) shows a log 10 -log 10 plot of the measure δ(τ) as a function of the nonlinear coupling g 0 for different values of μ c . As g 0 increases, the amount of non-Gaussianity first grows linearly in the logarithm, then plateaus as g 0 increases further. The same behaviour occurs for larger μ c , only more rapidly. This suggests that if we wish to increase the non-Gaussianity substantially, it will become increasingly difficult to do so by increasing g 0 . As such, focusing on increasing the coupling g 0 will only give marginal returns. Similarly, 3(b) shows log 10 d t ( ) as a function of increasing μ c for various values of g 0 . Since δ(τ) in this regime is proportional to F 2 | | , it follows that δ(τ) grows quadratically with the light-matter coupling in this regime. and μ m =0. The measure increases exponentially at first before it plateaus towards a constant value, which is the overall behaviour we observe in(a). . The behaviour of the measure of non-Gaussianity δ(τ) at τ=π for systems with constant nonlinear coupling g 0 starting in coherent states. (a) A log-log plot of δ(τ) versus the rescaled coupling g 0 . As g 0 increases, the state becomes more and more non-Gaussian, polynomially at first but then it quickly tends towards a constant value. (b) A log-plot of δ(τ) versus the coherent state parameter μ c for different values of g 0 . δ(τ) first increases quickly, then plateaus towards a single value.

Large coherent state parameters
We derived an asymptotic form of δ(τ) in (22) for the case 1 c m  | | , which we called d t( ). As argued before, the behaviour of the measure δ(τ) in this regime depends crucially on the distance of θ a from the value 2π. In our present case we have that a 3 q t for τ=1 and g a 0 2 q t -˜for τ ?1. The functions that we decided to ignore To see how well the asymptotic form d t( ) in (22) 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 1 c m| | , the asymptotic measure d t( ) well approximates the exact value of δ(τ). In fact, it becomes even more accurate as the optical coherent state parameter μ c increases, which is to be expected given the nature of the approximation. The asymptotic form also becomes more accurate once we also increase g 0 , as evident in figure 4(b). For g 10 0 2 = , the approximation is almost entirely accurate. This occurs because the function θ a increases with g 0 , which further suppresses the off-diagonal covariance matrix elements at the beginning and end of each cycle.
Let us discuss the fact that the measure recurs with 2 t p = for integer g 0 2 which we can now address analytically by examining the asymptotic covariance matrix elements in (20). We find that F=0 for all τ=2π n with integer n. This means that σ 42 =0 and that σ 22 =1. We also find that integer, we find that sin 2 0 a 2 q = and the final covariance matrix element is σ 11 =1. This results in diag 1, 1, 1, 1 s = ( ) which corresponds to a coherent state, which is fully Gaussian. As a result, the non-Gaussianity vanishes. When g 0 2 is not an integer, some non-Gaussianity will be retained, but the fact that F=0 will still result in a reduction at τ=2π.

Non-Gaussianity in open systems with constant coupling
Any realistic system will suffer from decoherence. In figure 5 we have plotted the non-Guassianity δ as a function of time for an optomechanical system with open dynamics. Here, the cavity state and the mechanics are both in initial coherent states (11). and μ m =0. We have chosen a low value of μ c to ensure high numerical accuracy of the simulation, as larger values quickly lead to numerical instabilities. We note that the non-Gaussianity d t ( ) tends towards a steady value, which is clear from the fact that the higher values of decoherence start to coincide around τ=5π. We also note that around τ=2π n, for integer n the inclusion of noise appears to temporarily increase the non-Gaussianity. This could, however, be due to the fact that the relative entropy measure cannot distinguish between non-Gaussianity induced as a result of genuinely nonlinear dynamics or as a result of classical mixing of the states [62]. We discuss this further in section 7. Similarly, in figure 5

Applications: time-dependent coupling
In all physical systems, such as optomechanical cavities, the confining trap is not ideal. This means that, in general, the coupling g t( ) is time-dependent as a consequence of, for example, trap instabilities. Timedependent variations such as phase fluctuations in the laser beam used to trap a levitated bead will modulate the coupling.
In this work, we want to exploit the possibility of controlling the coupling g t( ) 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 [63], and have been used to shuttle ions across large distances, typically for the purpose of quantum information processing [64,65]. See also section 7 for additional methods by which the coupling can be modulated.

Modelling the trap modulation
We shall model a time-dependent light-matter coupling g t( ) by assuming that the coupling has the simple form 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 (D.1).
We can now compute δ(τ) for this time-dependent coupling for initial coherent states, and we display the results in figure 6. In 6(a) we plot δ(τ) versus τ 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 W  we recover the time-independent solution, as expected. Interestingly, when 0 0 W ¹ we see that we can achieve higher values for the nonlinear measure δ(τ). This is especially pronounced as 1 0 W  , 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.

Trap modulation on resonance
The functions (D.1) contain denominators of the form Ω 0 −1. Therefore, among all possible values of Ω 0 , we can ask what happens on resonance, i.e. when Ω 0 =1. Figure 6(a) already provides evidence that the system should behave markedly differently.
e have plotted in figure 6 the exact measure of non-Gaussianity δ(τ) in the resonant case for initially coherent states and for different values of μ c . As anticipated, here we no longer have recurrent behaviour. Instead, the nonlinearity increases as ln t. 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 6(c). We note that δ(τ) oscillates increasingly rapidly with larger ò but with decreasing amplitude for increasing τ, as F g 2 0 2 2 2  t | |˜becomes the dominant term for τ?1. As already noted, it is evident from figure 6 that the non-Gaussianity increases continuously. The nonlinear coupling in the Hamiltonian is derived by considering the effect of photon pressure on the mechanical element. Given that the overall photon number a a á ñ † is conserved, the coupling acts as a photon number displacement. If this coupling is time-dependent, this means that the photon pressure displaces with a time-dependence. When this occurs at the resonance, this linear displacement grows linearly in time. See also [48] for further insight once the rotating wave approximation has been applied.

Large coherent state parameters at resonance
Using the explicit form of the coefficients (24) for large τ and therefore we expect, as it happened in section 4.2, that most covariance matrix elements will vanish and will not contribute to the asymptotic form of δ(τ). This observation allows us to compute the symplectic eigenvalues, which read  (21). We again stress that we have retained the exact expression for σ 11 to capture some crucial features of the non-Gaussianity, such as δ(0)=0.
In figure 7, we compare the exact measure δ(τ) at resonance with the asymptotic form derived in (22). The solid lines represent the exact measure δ(τ) and the dashed lines represent the asymptotic expression. In figure 7(a) 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 (B.6) that suppress some elements for large μ c , unless θ a =n π. Similarly, in figure 7(b) 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.

Open system dynamics at resonance
If it is possible to continuously increase the non-Gaussianity, the system might have a certain tolerance to noise. That is, there is a level of noise at which the non-Gaussianity essentially reaches a steady-state. In figure 8 we have plotted the non-Gaussianity δ as a function of time for different values of photon and phonon decoherence. Figure 8(a) shows the system at resonance with photons leaking from the cavity with a rate c c m k k w = for parameters μ c =0.1, g 1 0 = , ò=0.5 and μ m =0. We note that 0.3 c k = yields what is essentially a steady-state of the non-Gaussianity. In figure 8(b) we note the same behaviour but for phonon decoherence with rate m m m k k w = .

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 coupling g 0 , the non-Gaussianity d t ( ) scales differently in two contrasting regimes: (i)for a weak optical input coherent state c m | |, the nonlinear character of the state grows as g ln Crucially, we also find that the amount of non-Gaussianity can be continuously increased by driving the light-matter coupling at mechanical resonance. This becomes especially useful in the presence of noise. 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. First, however, we will discuss the measure of non-Gaussianity that we have used in this work.

Choice of measure
In this work, we chose to work with a relative entropy measure of non-Gaussianity (see section 3) which was first defined in [50]. This measure has previously been extensively used to compute the non-Gaussianity of various states [66], as well as in an experimental setting where single photons were gradually added to a coherent state to increase its non-Gaussian character [62]. Several additional measures for the quantification of non-Gaussianity have been proposed in the literature, linking it to the Hilbert-Schmidt distance [67] or to quantum correlations [68]. Specifically, the relative entropy measure was shown to be more general than the Hilbert-Schmidt measure [67]. Furthermore, a connection has been put forward between the non-Gaussianity of a state and its Wigner function [69], and similarly there appears to be an intrinsic link between the quantum Fisher information and the lowest amount of non-Gaussianity of a state [70].
Most crucially, this measure is not upper-bounded. This means that, as opposed to, for example, an entanglement measure where the notion of a maximally entangled state is well-defined, there is no such thing as a maximally non-Gaussian state. This is reflected by our results, where taking μ c to infinity yields lim c d = ¥ m ¥ . As such, it is only possible to state that one state is more non-Gaussian than another. However, for pure states, there is the relation of the measure to the Hilbert-Schmidt measure. As such, the non-Gaussianity d t ( ) of pure states has strong operational implications [8]. For mixed states, the operational meaning is not clear because the measure cannot detect the difference between classical mixtures of Gaussian states, which can be easily prepared by classical mixtures of Gaussian states, and inherent non-Gaussianity due to some nonlinear evolution of pure states [62]. This means that the measure often needs to be used together with a measure of non-classicality, such as the negativity of the Wigner function. We know from previous work [22,25] that for a constant coupling, the system is maximally entangled at τ=π, which satisfies the occurrence of non-classicality in conjunction with the non-Gaussianity. The state is however fully disentangled at τ=2π, and in the case of open system dynamics, this feature of the measure becomes apparent. We note that the non-Gaussianity plotted in figure 5 spikes at times τ=2π n for integer n, which is when we usually have no entanglement. This implies that the addition of non-Gaussianity most likely comes from a classical mixture of coherent states that have slightly decohered.

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 singlephoton optomechnical coupling regime, the light-matter coupling g is small compared to the resonant frequency ω m and the optical decoherence rate κ c . Such experiments usually involve a strong laser drive, which tends to wash out the nonlinearity. In the strong single-photon coupling regime, nonlinear effects are in practice small but more significant. Under these conditions, a single photon displaces the mechanical oscillator by more than its zero-point uncertainty and weak optical fields tend to be used [71]. 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.
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.
Let us also briefly discuss our results with regard to linearised dynamics. This linearisation of dynamics is fundamentally different to the scenarios considered in this work. When linearising the dynamics, the system is opened and the field operators â are treated as flucutations around a strong optical field as such: a a a a  = + ¢ˆˆ, where a¢ are the fluctuations. In this work, we have retained the nonlinear dynamics, even when considering open system dynamics. Thus, while we observe that a large coherent state parameter μ c increases the non-Gaussianity, we cannot generalise this result to the linearised dynamics.

Methods of 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 coupling g t( ) is modulated. An explanation of this phenomena was provided in [48]. Consider the force F  exerted by the photons on the mechanics. For a number of n photons, this force is proportional to F n 1 2 µ +  ( ) , where 1/2 comes from the zero-point energy. When the light-matter coupling is constant, this force is constant, and thus we see the periodic evolution. However, when we modulate g t( ), the photon-pressure force F  acts periodically on the mechanics, and is amplified when pushing in tandem with the mechanical resonance.
While engineering the modulation is challenging, we shall explore several methods that can achieve it. The question is whether the modulation can be performed at mechanical resonance. As a basis for this discussion, we present a derivation of a time-dependent light-matter coupling for levitated nanobeads in appendix F, which is based on the work in [72]. There are several practical ways in which one may envisage to increase the nonlinearity 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 F, this phase is denoted by j(τ) 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 1 sin 2 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 [64,65] using Paul traps, which are customarily used for ions but which have also recently been used for trapping nanoparticles [63,73,74]. 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. To date, the micromotion is generally smaller than the mechanical frequency, ω d ω m , but current experimental efforts appear promising.
There are potentially many more ways in which the light-matter coupling could be modulated, including with optomechanically induced transparecny [75,76] and by using the Kerr effect to change the refractive index of the oscillator.
We conclude that the enhancement of the nonlinearity 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.

Detecting and measuring non-Gaussianity in optomechanical systems
In practise, how would one proceed to measure the amount of non-Gaussianity in the laboratory? As shown in [62], the measure of non-Gaussianity used in this work has been measured for the addition of single photons to a coherent state. This requires full state tomography and is thus an expensive process. There are however others ways to proceed. In [77] a witness of non-Gaussianity was proposed based on bounding the average photon number in the system from above. While they apply to a single system, they can probably be extended to bipartite systems as well.
Finally, we here suggest a simple method by which non-Gaussianity can be detected for pure states. We note that the von Neumann entropy S AB r (ˆ) of a bipartite state AB r is bounded by S S S (ˆ)|, through the Araki-Lieb inequality [78] where A r and B r are the reduced states of the optical and mechanical subsystems, respectively. Therefore, the measure of non-Gaussianity δ(τ) that we defined in (15) is lower-bounded by In this sense, this reduced measure acts as a sufficient (but not necessary) condition for non-Gaussianity. That is, finding that the measure is non-zero does tell us that the state is non-Gaussian, however it does not tell us the full magnitude of the non-Gaussianity. Furthermore, to compute this measure, one would still have to measure the second moments of the optical and mechanical subsystems. This does, however, require fewer measurements than full state tomography on the joint optical and mechanical system.

Conclusions
We have quantified the non-Gaussianity of initially Gaussian coherent states evolving under the standard, timedependent optomechanical Hamiltonian. We used a measure of non-Gaussianity based on the relative entropy of a state to characterise the deviation from Gaussianity of the full system. Our techniques allowed us to derive asymptotic expressions for small and large optical coherent-state amplitudes, see equation (19) and equation (22) respectively. We found that for coherent states with amplitude 1 c  m | | , the amount of non-Gaussianity grows logarithmically with the input average number of excitations c m | | and with the light-matter coupling. At resonance, we find that the non-Gaussianity is further enhanced by a logarithmic scaling with the time of interaction.
An important and promising aspect of our study consists in showing that the amount of non-Gaussianity in the system can be continuously increased by driving a time-modulated optomechanical coupling at mechanical resonance. This allows us to circumvent the usual periodic increase and decrease of non-Gaussianity, and we find that this behaviour effectively yields a non-Gaussian steady-state in the presence of noise. As such, this points to a practically accessible, mechanism to enhance the nonlinear character of optomechanical dynamics at a given light-matter coupling strength. We point out that certain systems, such as hybrid-trap systems, are particularly well-suited for this purpose, as their light-matter interaction is naturally modulated due to the trap characteristics. Finally, we also conclude that the mechanical system does not have to be cooled to the ground state in order to access significant amounts of the non-Gaussianity.
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.

Data availability statement
The figures in this work were generated using Mathematica. The numerical results in sections 5.3 and 6.4 were obtained using the Python library QuTiP. The code used to generate the data can be found at https://github. com/sqvarfort/QM-Nonlinearities.

Appendix A. Derivation of the dynamics and general tools
In this appendix, we will derive the coefficients in (7) that determine the time-evolution of the system. This follows the derivation in [49]. We will also show the explicit time-dependence of the second moments and discuss some methods related to computing the symplectic eigenvalues of the covariance matrix.

B.1. Properties of of the nonlinear Hamiltonian
Firstly, we remind the reader that the laboratory time t is rescaled by ω m . Finding a simple expression for U t ( ) is straight-forward when the light-matter coupling g g m w = is not time-dependent. If g g m m t w t w =( ) ( ) 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 [49]. We note that compared with [49], we have here swapped the definition of â and b, and we have a minus-sign in front of g t( ).
For the time-dependent Hamiltonian Ĥ in (1), the time-evolution operator is given by is the time-ordering operator. The basis for decoupling the operator is finding a Lie algebra of generators G î that induce the time-evolution. This Lie algebra must be closed under commutation, that is, either G G G , where c is a scalar. This will allow for the terms in U 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 as where F j are coefficients corresponding to each of the generators G ĵ . Our task is now to find the coefficients F j . We begin by defining the operators G ĵ in the algebra: 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: . This is the equation that determines the coefficients. We can now commute all the operators through, where we find and τ=ω m t. Depending on the form of g t( ), we can now use these equations to find a simplified form of Û t ( ).

B.2. Computing determinants of symplectic matrices
When computing the amount of non-Gaussianity in (15), it is useful to consider the symplectic eigenvalues of a Gaussian state [59]. In short, for an arbitrary covariance matrix s, they are the eigenvalues of the matrix i s W , ) is the symplectic form. There are other ways to define the symplectic eigenvalues though. In the following, we have to switch the basis of the operators to a more convenient one, but this does not affect the final result. The correct definition can be found in [56]. Let us write an arbitrary covariance matrix s in the particular basis q p q p , , , his expression can also be rewritten in more compact notation as We now compute the expectation values for initial optical coherent states and coherent and thermal coherent states of the mechanics.

C.1. Mechanical coherent states
As a result, the expectation values for the initial mechanical states as a coherent thermal state are given by The resulting covariance matrix elements can be computed from here.

Appendix C. Derivation of the asymptotic form of the symplectic eigenvalues
The symplectic eigenvalues ν ± have the expression ν ± =1+δν ± . We would like to see what is the form of the function f x ln ln when we compute f (ν ± ) and δν ± =1 or δν ± ?1. In the first case, δν ± =1, and we have

Appendix D. Coefficients for time-dependent light-matter coupling
In this appendix, we compute the coefficients used in section 6. Starting from (A.9), we assume that the coupling has the functional form g g 1 sin 0 0 . The algebra is straightforward, although cumbersome, and it leads us to the expressions: t 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 1 0 W  , which occurs when ω 0 =ω m .
It is straight-forward to see how the terms F N B a +ˆa nd F N B a -ˆs implify as . The long expression for F N a 2 is more challenging. We note that the terms independent of ò remain unchanged with Ω 0 . Thus, at resonance, these coefficients read: The classical contribution, x ( ) will yield a trapping frequency, while the operator terms E x ( ) will yield the light-matter interaction term for the levitated sphere. The cross-term, E x x 2  ( ) ( ) will generate elastic scattering processes inside the cavity which converts cavity photons and tweezer photons into free modes [72]. 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.

F.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 [79].
If we denote the radius of the cylinder by r, we can write down the trapping field as . 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 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

F.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 [72] 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 p ( ) , 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 j is a generic phase which determines the minimum of the potential seen by the bead. For laser-trapped nanobeads, this phase can be made timedependent, whereas for a Paul trap, it is static. We will leave out the time-dependence 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 x y z r , , = ( )leads to