Enhanced continuous generation of non-Gaussianity through optomechanical modulation

TL;DR

This paper investigates how to continuously enhance non-Gaussianity in quantum optomechanical systems by modulating the light-matter interaction, providing analytical insights and practical mechanisms for non-Gaussian state generation.

Contribution

It introduces a method to increase non-Gaussianity through resonant modulation of the optomechanical coupling, with analytical expressions and robustness against thermal noise.

Findings

Non-Gaussianity can be increased by modulating the coupling at mechanical resonance.

Initial thermal phonon occupation has minimal impact on non-Gaussianity.

Continuous non-Gaussianity enhancement is feasible even with noise.

Abstract

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-Gaussianity 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.

Equations

$$\hat{H}$$

$$\hat{H}/(\hslash {\omega }_{\mathrm{m}})$$

$$\hat{U}(\tau):=\overset{\leftarrow}{\mathcal{T}}\,\exp\biggl{[}-i\,\int_{0}^{\tau}d\tau^{\prime}\,\hat{H}(\tau^{\prime})\biggr{]},$$

$$\hat{U}(\tau )=\prod _j{\hat{U}}_j(\tau )=\prod _j{e}^{-i{F}_j(\tau ){\hat{G}}_j},$$

$$\hat{U}(t)=e^{-i{\hat{N}}_b\tau }{e}^{-i{F}_{{\hat{N}}_a^2}{\hat{N}}_a^2}{e}^{-i{F}_{{\hat{N}}_a{\hat{B}}_{+}}{\hat{N}}_a{\hat{B}}_{+}}{e}^{-i{F}_{{\hat{N}}_a{\hat{B}}_{-}}{\hat{N}}_a{\hat{B}}_{-}},$$

$${\hat{N}}_a$$

$${\hat{N}}_a{\hat{B}}_{+}$$

$$F_{{\hat{N}}_a^2}$$

$$F_{{\hat{N}}_a{\hat{B}}_{+}}$$

$$F_{{\hat{N}}_a{\hat{B}}_{-}}$$

$${\theta }_a$$

$$F$$

$$F_{{\hat{N}}_a^2}$$

$$F_{{\hat{N}}_a{\hat{B}}_{+}}$$

$$F_{{\hat{N}}_a{\hat{B}}_{-}}$$

$$\hat{U}(\tau )=e^{-i{\hat{N}}_b\tau }{e}^{i{\tilde{g}}_0^2[1-\text{sinc}(2\tau )]\tau {\hat{N}}_a^2}{e}^{i{\tilde{g}}_0\sin (\tau ){\hat{N}}_a{\hat{B}}_{+}}{e}^{i{\tilde{g}}_0(1-\cos (\tau )){\hat{N}}_a{\hat{B}}_{-}}.$$

$$|\Psi (0)⟩=|{\mu }_{\mathrm{c}}⟩\otimes |{\mu }_{\mathrm{m}}⟩.$$

$${\hat{\rho }}_{\mathrm{th}}=\frac{1}{\bar{n}\pi }\int {\mathrm{d}}^2\beta e^{-|\beta {|}^2/\bar{n}}|\beta ⟩⟨\beta |,$$

$$\hat{\rho }(0)=|{\mu }_{\mathrm{c}}⟩⟨{\mu }_{\mathrm{c}}|\otimes {\hat{\rho }}_{\mathrm{th}}$$

$$\dot{\hat{\rho }}=-i\left[\hat{H},\hat{\rho }\right]+{\hat{L}}_{\mathrm{c}}\hat{\rho }{\hat{L}}_{\mathrm{c}}^{†}+{\hat{L}}_{\mathrm{m}}\hat{\rho }{\hat{L}}_{\mathrm{m}}^{†}-\frac{1}{2}\{{\hat{L}}_{\mathrm{c}}^{†}{\hat{L}}_{\mathrm{c}},\hat{\rho }\}-\frac{1}{2}\{{\hat{L}}_{\mathrm{m}}^{†}{\hat{L}}_{\mathrm{m}},\hat{\rho }\},$$

$$\delta (\tau ):=S({\hat{\rho }}_{\text{G}}(\tau ))-S(\hat{\rho }(\tau )),$$

$$\delta (\tau )\approx$$

$$\delta (\tau )\sim (N_{+}+N_{-})\ln x.$$

$${\nu }_{+}\sim$$

$${\nu }_{-}\sim$$

$$\delta (\tau )\sim -\left(1+\left(1-2e^{-|F{|}^2}\right)|F{|}^2\right)|{\mu }_{\mathrm{c}}{|}^2\ln |{\mu }_{\mathrm{c}}|,$$

$${\sigma }_{11}$$

$${\sigma }_{22}$$

$${\sigma }_{42}$$

$${\nu }_{+}$$

$${\nu }_{-}$$

$$\delta (\tau )\sim s_V\left(1+2|{\mu }_{\mathrm{c}}{|}^2\left(1-e^{-4|{\mu }_{\mathrm{c}}{|}^2{\sin }^2{\theta }_a/2}{e}^{-|F{|}^2}\right)\right)+s_V(\sqrt{4|{\mu }_{\mathrm{c}}{|}^2|F{|}^2+1}),$$

$$\tilde{g}(\tau )={\tilde{g}}_0\left(1+ϵ\sin ({\Omega }_0\tau )\right).$$

$$F_{{\hat{N}}_a^2}$$

$$\displaystyle\quad\quad\quad\quad\quad\quad\quad+\epsilon^{2}\,\bigl{(}6\,\tau-4\sin(2\,\tau)+\sin(2\,\tau)\,\cos(2\,\tau)\bigr{)}\bigr{]}$$

$$F_{{\hat{N}}_a{\hat{B}}_{+}}$$

Full text

1 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 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 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 general-dyne detection, are classically simulatable, as their Wigner functions may be mimicked by classical probability distributions [19].111We 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.

In optomechanical systems [20], 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 [21, 22, 23]. The nonlinear features of optomechanical systems have been frequently explored in the context of metrology, such as force sensing [24] and gravimetry [25, 26]. In particular, the nonlinear coupling enables the creation of optical cat states in the form of superpositions of coherent states [21, 22]. These cat states can also be transferred to the mechanics [27], which opens up the possibility of using massive superpositions for testing fundamental phenomena such as collapse theories [28] and, potentially, signatures of gravitational effects on quantum systems at low energies [29, 30]. In addition to cat states, other non-Gaussian states such as compass states [31, 32] and hypercube states [33], and also all been found to have excellent sensing capabilities. This combination of sensing with nonlinear state and fundamental applications makes it imperative to explore the nonlinear properties of the mechanical systems.

A number of different optomechanical systems have been experimentally implemented, including Fabry–Pérot cavities with a moving-end mirror [34], as shown in Figure 1, levitated nano-diamonds [35, 36], membrane-in-the-middle configurations [37] and optomechanical crystals [38, 39]. While several experiments have demonstrated genuine nonlinear behaviour (see for example [40, 41, 42, 43]), most experimental settings can however be fully modelled with linear dynamics [44, 43], and it is generally difficult to access the fully nonlinear regime. As a result, 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, by using a large-amplitude, strongly detuned mechanical parametric drive [45], or by modulating the spring constant [46]. 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 [47].

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 further increased? One such proposal was put forward in [48] where it was suggested that the nonlinearity in electromechanical systems could be enhanced by several orders of magnitude by modulating the light–matter coupling. This is achieved by driving the system close to mechanical resonance and takes a simple form in the rotating-wave approximation. Here, we seek to fully quantify the non-Gaussianity of the exact, nonlinear optomechanical state for both ideal and open systems. 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 time-dependent systems [49], and employ a specific measure of non-Gaussianity based on the relative entropy of the state [50]. 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 find that the non-Gaussianity scales logarithmically with the coherent state parameter of the optical system, and we illustrate how this behaviour differs for small and large coherent states. Most importantly, however, we find that the non-Gaussianity of the state can be continuously increased by modulating the light–matter coupling strength at the mechanical resonant frequency. Such a continuous increase might prove especially useful given open system dynamics, where it would allow for the creation of an effectively non-Gaussian steady state. Moreover, we find 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  A, B, C, and D

2 Dynamics

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

2.1 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 $\hat{a}$ and $\hat{b}$ of the cavity and mechanical modes respectively, obey the canonical commutation relations $[\hat{a},\hat{a}^{\dagger}]=[\hat{b},\hat{b}^{\dagger}]=1$, 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:

[TABLE]

where $\omega_{\mathrm{c}}$ and $\omega_{\mathrm{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 $\omega_{\mathrm{m}}$, therefore introducing the dimensionless time $\tau:=\omega_{\mathrm{m}}\,t$, the dimensionless frequency $\Omega:=\omega_{\mathrm{c}}/\omega_{\mathrm{m}}$, and the dimensionless coupling $\tilde{g}(\tau):=g(t\omega_{\mathrm{m}})/\omega_{\mathrm{m}}$. This choice will prove convenient throughout the rest of this work, and dimensions can be restored when necessary by multiplying by $\omega_{\mathrm{m}}$. This renormalisation effectively is equivalent to the use of time $\tau$ and the Hamiltonian

[TABLE]

To determine the action of this Hamiltonian on initial states, we now proceed to solve the dynamics induced by (2).

2.2 Time evolution of the system

We now need an expression for the time evolution operator $\hat{U}(\tau)$ for a system evolving with (2). The unitary time-evolution operator reads

[TABLE]

where $\overset{\leftarrow}{\mathcal{T}}$ is the time-ordering operator [52]. This expression simplifies dramatically when the Hamiltonian $\hat{H}$ is independent of time, in which case one simply has \hat{U}(\tau)=\exp\bigl{(}-i\,\hat{H}\,\tau\bigr{)}. As we will here consider time-dependent light–matter couplings $\tilde{g}(\tau)$, 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 $\hat{U}_{j}$. This is possible if there exists a finite set of generators $\hat{G}_{j}$ that form a closed Lie algebra under commutation [52]. We thus write the evolution operator (3) as

[TABLE]

where $F_{j}(\tau)$ are generally time-dependent coefficients determining the influence of the generator $\hat{G}_{j}$ 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 $\hat{a}$ and $\hat{b}$ have been swapped and that in our case the coupling $\tilde{g}(\tau)$ is preceded by a minus-sign. For clarity, we therefore present the full derivation for the case considered here in A.

With these techniques, we find that the time-evolution operator $\hat{U}(\tau)$ can be cast into the convenient form

[TABLE]

where the operators, given by,

[TABLE]

form a closed Lie algebra under commutation, and where the coefficients that determine the evolution in (5) are given by

[TABLE]

Note that in (5), we have transformed into a frame rotating with $\Omega_{\mathrm{c}}\,\hat{N}_{a}$ in order to neglect the phase-term $e^{-i\Omega\tau}$. Given an explicit form of $\tilde{g}(\tau)$, it is then possible to write down a full solution for $\hat{U}(\tau)$. 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:

[TABLE]

These quantities will be useful when discussing features of the non-Gaussianity.

2.3 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 $\tilde{g}(\tau)=\tilde{g}_{0}$. The functions (2.2) simplify to

[TABLE]

These coefficients allow us to write the time evolution operator $\hat{U}(\tau)$ as

[TABLE]

This 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)).

2.4 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.

• i)

Coherent states. We start by considering the case when both the optical and the mechanical modes are in a coherent state, which we denote $\ket{\mu_{\mathrm{c}}}$ and $\ket{\mu_{\mathrm{m}}}$ respectively. These states satisfy the relations $\hat{a}\ket{\mu_{\mathrm{c}}}=\mu_{\mathrm{c}}\ket{\mu_{\mathrm{c}}}$ and $\hat{b}\ket{\mu_{\mathrm{m}}}=\mu_{\mathrm{m}}\ket{\mu_{\mathrm{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 $\mu_{\mathrm{m}}=0$. The initial state $\ket{\Psi(0)}$ of the compound system will therefore be

[TABLE]

• i)

Thermal coherent states. 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 state is obtained simply by integrating over the coherent state parameter with an appropriate kernel [55]. We define the thermal state $\hat{\rho}_{\mathrm{th}}$ as

[TABLE]

where $\bar{n}$ is the average thermal phonon occupation of the state, $\ket{\beta}$ is a coherent state, and the integration occurs over the full complex space. Assuming that the optical mode in the coherent state $\ket{\mu_{\mathrm{c}}}$, the full initial state is therefore given by

[TABLE]

By starting in an initial Gaussian state, we ensure that any non-Gaussianity revealed by our work is due to the nonlinear 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 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 $\hat{\rho}_{\mathrm{th}}$. See Section 7 for a discussion of the properties of the relative entropy measure.

2.5 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

[TABLE]

where $\hat{L}_{\mathrm{c}}$ and $\hat{L}_{\mathrm{m}}$ are the Lindblad operators for the optics and mechanics, respectively. To model photon and phonon decay, we assume that $\hat{L}_{\mathrm{c}}=\sqrt{\kappa_{\mathrm{c}}}\,\hat{a}$ and $\hat{L}_{\mathrm{m}}=\sqrt{\kappa_{\mathrm{m}}}\,\hat{b}$, where $\kappa_{\rm{c}}$ is the optical decoherence rate and $\kappa_{\rm{m}}$ is the phonon decoherence rate.

While analytic solutions for this particular choice of $\hat{L}_{\mathrm{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.

3 Measures of deviation from Gaussianity

Given a Hamiltonian $\hat{H}$, and an initial Gaussian state $\hat{\rho}(0)$, we ask the following question: can we quantify how much the state $\hat{\rho}(\tau)$ deviates from a Gaussian state at time $\tau$? 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 [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 $\delta$, 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].

3.1 Measures of deviation from Gaussianity: definition

Let us detail here the construction of the non-Gaussianity quantifier $\delta(\tau)$ for our nonlinear dynamics. First, our initial state $\hat{\rho}(0)$ evolves into the state $\hat{\rho}(\tau)$ at time $\tau$. With our full solutions for the dynamics in Section 2.2, we can find analytic expressions for the first and second moments of $\hat{\rho}(\tau)$. Then, we construct a state $\hat{\rho}_{\textrm{G}}(\tau)$, which is the Gaussian state defined by the first and second moments that coincide with those of $\hat{\rho}(\tau)$. Now, we recall that a Gaussian state is fully defined by its first and second moments. Therefore, if two Gaussian states $\hat{\rho}_{1}$ and $\hat{\rho}_{2}$ have equal first and second moments they are the same state [56, 59]. In general, our state $\hat{\rho}(\tau)$ at time $\tau$ will not be Gaussian and therefore cannot be specified fully by its first and second moments. This implies that we can introduce a measure $\delta(\tau)$ that quantifies how $\hat{\rho}(\tau)$ deviates from $\hat{\rho}_{\textrm{G}}(\tau)$:

[TABLE]

where $S(\hat{\rho})$ is the von Neumann entropy of a state $\hat{\rho}$, defined by $S(\hat{\rho}):=-\Tr(\hat{\rho}\,\ln\hat{\rho})$. This measure has been shown to capture the intrinsic non-Gaussianity of the system, and it vanishes if and only if $\hat{\rho}(\tau)$ is a Gaussian state [50]. In other words, if at all times $\delta(\tau)=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(\hat{\rho}(\tau))=S(\hat{\rho}(0))$. If we start from a pure state, then $S(\hat{\rho}(0))=0$ and $\delta(\tau)=S(\hat{\rho}_{\textrm{G}}(\tau))$. 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 $\hat{\rho}_{\textrm{G}}(\tau)$ 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 $\hat{\rho}_{\mathrm{G}}$ involves finding the first of second moments of $\hat{\rho}$ anyway. While we could compute the entropy for $\hat{\rho}_{\mathrm{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\times 4$ covariance matrix $\boldsymbol{\sigma}(\tau)$ of the state $\hat{\rho}_{\textrm{G}}(\tau)$. This matrix contains the second moments of the state $\hat{\rho}(\tau)$ which in our specific choice of basis $\hat{\mathbb{X}}=(\hat{a},\hat{b},\hat{a}^{\dagger},\hat{b}^{\dagger})^{\mathrm{T}}$ is defined through its elements $\sigma_{nm}(\tau):=\langle\{\hat{X}_{n},\hat{X}_{m}^{\dagger}\}\rangle-2\,\langle\hat{X}_{n}\rangle\langle\hat{X}_{m}^{\dagger}\rangle$, where $\{\bullet,\bullet\}$ is the anti-commutator, we have defined the expectation value of an operator $\langle\bullet\rangle:=\Tr\{\bullet\,\rho_{\textrm{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(\hat{\rho}_{\textrm{G}}(\tau))$ we require the symplectic eigenvalues $\{\pm\nu_{+}(\tau),\pm\nu_{-}(\tau)\}$ of $\boldsymbol{\sigma}(\tau)$, where the property $\nu_{\pm}(\tau)\geq 1$ holds for all states. The symplectic eigenvalues can be computed by finding the eigenvalues of the object $i\,\boldsymbol{\Omega}\,\boldsymbol{\sigma}(\tau)$, where $\boldsymbol{\Omega}$ is the symplectic form given by $\boldsymbol{\Omega}=\textrm{diag}(-i,-i,i,i)$ in this basis. The von Neumann entropy $S(\boldsymbol{\sigma})$ is then given in this formalism by $S(\boldsymbol{\sigma})=s_{V}(\nu_{+})+s_{V}(\nu_{-})$, where the binary entropy $s_{V}(x)$ is defined by $s_{V}(x):=\frac{x+1}{2}\ln\left(\frac{x+1}{2}\right)-\frac{x-1}{2}\ln\left(\frac{x-1}{2}\right)$. In summary, the state $\hat{\rho}(\tau)$ is non-Gaussian at time $\tau$ if and only if $\delta(\tau)>0$.

3.3 Measures of deviation from Gaussianity: General behaviour

We shall now infer some general characteristics of the measure of non-Gaussianity $\delta(\tau)$. 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 $\nu_{\pm}$ satisfy $\nu_{\pm}\geq 1$ [59]. Therefore, we can conveniently write $\nu_{\pm}=1+\delta\nu_{\pm}$, where $\delta\nu_{\pm}\geq 0$ captures any deviation from purity. If the state is pure at time $\tau=0$ then it follows that $\nu_{\pm}(0)=1$ [59]. If the evolution is linear, then it is also the case that $\nu_{\pm}(\tau)=1$ for all $\tau$. For closed dynamics, the symplectic eigenvalues may only change if the evolution is non-linear. In this case, we would define $\nu_{\pm}=\nu_{0\pm}+\delta\nu_{\pm}$ with $\nu_{0\pm}>1$. Then, we would have that $\nu_{\pm}(0)=\nu_{0\pm}$ and, again, linear evolution would imply that $\nu_{\pm}(\tau)=\nu_{0\pm}$. The preceding statements imply that $\delta\nu_{\pm}$ are functions of the nonlinear contributions alone. Thus, when the non-linearity tends to vanish, then $\delta\nu_{\pm}\rightarrow 0$. Among the possible asymptotic regimes we have that $\delta\nu_{\pm}\rightarrow+\infty$ or that it becomes constant.

These observations are important. To understand their implications we use the expression $\nu_{\pm}=1+\delta\nu_{\pm}$ to compute the general deviation from Gaussianity as $\delta(\tau)=s_{V}(1+\delta\nu_{+})+s_{V}(1+\delta\nu_{-})$. Using this form, we see that in the nearly linear (Gaussian) regime with only small contributions from the nonlinear dynamics, we will have $\delta\nu_{\pm}\ll 1$ and therefore

[TABLE]

On the contrary, in the highly nonlinear (non-Gaussian) regime we have $\delta\nu_{\pm}\gg 1$ and therefore $\delta(\tau)\approx\ln\frac{\delta\nu_{+}}{2}+\ln\frac{\delta\nu_{-}}{2}$. If the symplectic eigenvalues depend on a large parameter $x\gg 1$, then one will in general find that they have the asymptotic form $\nu_{\pm}\sim x^{N_{\pm}}\sum^{N_{\pm}}_{n=0}\nu_{\pm}^{(n)}\,x^{-n}$ for some appropriate real coefficients $\nu_{\pm}^{(n)}$, where $N_{\pm}$ constitutes the upper limits of the sum [61].

A careful asymptotic expansion of the measure of nonlinearity in this regime gives

[TABLE]

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 C.

4 General results

We now proceed to evolve the initial state $\hat{\rho}(0)$ with the evolution operator $\hat{U}(\tau)$ in (5). To compute the amount of non-Gaussianity of this state $\delta(\tau)$, we must first find the elements $\sigma_{nm}$ of the covariance matrix $\boldsymbol{\sigma}$, which we construct form the first and second moments of $\hat{\rho}(\tau)$. We do so in full generality, meaning that the light–matter coupling $\tilde{g}(\tau)$ can take the form of any time-dependent function.

We have computed the second moments and the covariance matrix $\boldsymbol{\sigma}$ in B. The second moments for the mechanical being in an initial coherent state and in an initial coherent thermal state can be found in (B.1) and (B.2) respectively. These can then be used to compute the covariance matrix elements $\sigma_{nm}$ where $n,m$ take values $0,1,2,3$. We have explicitly computed the elements of $\mathbf{\sigma}$ for a mechanical coherent state in (B.1).

Our challenge now is to compute the symplectic eigenvalues $\nu_{\pm}$, given the expressions (B.1). The process of computing the eigenvalues can be simplified by using the expression $2\,\nu_{\pm}^{2}=\Delta\pm\sqrt{\Delta^{2}-4\,\textrm{det}(\boldsymbol{\sigma})}$, which is based on the existence of symplectic invariants [59]. The definition of $\Delta$ is given in A.2. The full analytic expression for the symplectic eigenvalues $\nu_{\pm}$ of $\boldsymbol{\sigma}$ in (B.1) is too long and cumbersome to be printed here. It also does not yield any immediate insight into the behaviour of $\delta(\tau)$. 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.

4.1 Asymptotic behaviour for a small optical coherent state parameter

We begin by looking at the case where $|\mu_{\mathrm{c}}|^{2}\ll 1$ and where the mechanics is initially in a coherent state. Here, one can take the covariance matrix elements (B.1) and, after some algebra, show that the perturbative expansion of the symplectic eigenvalues gives

[TABLE]

This implies that the behaviour of $\delta(\tau)$ for small $|\mu_{\mathrm{c}}|$ goes as

[TABLE]

in perfect agreement with (16). This approximation suggests that $\delta(\tau)$ scales with $\sim|F|^{2}|\mu_{\mathrm{c}}|^{2}\ln|\mu_{\mathrm{c}}|$ to leading order.

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.

4.2 Asymptotic behaviour for a large optical coherent state parameter

We now investigate the case where $|\mu_{\mathrm{c}}|\gg 1$. Our goal is to derive an analytic expression for the non-Gaussianity that can be used to analyse the overall features of $\delta(\tau)$. Before making any quantitative evaluation, we recall that the measure will have the form (17), where now $x\equiv|\mu_{\mathrm{c}}|^{2}$. Let us proceed to demonstrate this result analytically for this specific case.

For large $\mu_{\mathrm{c}}$ and for the mechanics in the ground-state $\mu_{\mathrm{m}}=0$, it is clear that whenever $\theta_{a}\neq 2\pi n$ for integer $n$, the matrix elements $\sigma_{31}$, $\sigma_{21}$ and $\sigma_{41}$ in (B.1) vanish, due to the exponentials containing the factor $|\mu_{\mathrm{c}}|^{2}$. Therefore, far (enough) from the times where $\theta_{a}=2\pi n$ we are left with the following covariance matrix elements

[TABLE]

and all other elements are zero. We have kept the full expression for $\sigma_{11}$ because it reproduces some key elements of $\delta$ 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 (4.2), which is a constant in the case of $\sigma_{11}$ and $\sigma_{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 $\braket{a}$ changes slightly and $\sigma_{11}$ in (4.2) will look different. However, if we approximate $\sigma_{11}$ as $\sigma_{11}\approx 1+2\,|\mu_{\mathrm{c}}|^{2}$, which follows from that $\braket{a}\sim 0$ for very large $\mu_{\mathrm{c}}$, the non-zero covariance matrix elements of the coherent thermal mechanical state are the same as in (4.2). 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

[TABLE]

We note that both eigenvalues grow with $|\mu_{\mathrm{c}}|$, as expected from our analysis in Section 3.3. The amount of non-Gaussianity for large $\mu_{\mathrm{c}}$ is now given by the following expression

[TABLE]

which scales asymptotically as $\delta(\tau)\sim\tilde{\delta}(\tau):=4\ln\,|\mu_{\mathrm{c}}|$, in perfect agreement with (17).

Note that (22) is also valid for a time-dependent light–matter coupling $\tilde{g}(\tau)$. In all cases, the nonlinearity grows as $\ln{|\mu_{\mathrm{c}}|}$ to leading order. In Sections 5.2 and 6.3 we will compare the asymptotic measure $\tilde{\delta}(\tau)$ with the full measure $\delta$ for different cases.

5 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: $\tilde{g}(\tau)=\tilde{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 (2.3), and we note that the function $F$, defined in (2.2), which appears in the covariance matrix elements $\sigma_{nm}$ is now given by $F=\tilde{g}_{0}\,(1-e^{-i\tau})$.

We now proceed to compute the exact measure of non-Gaussianity $\delta(\tau)$ for constant coupling $\tilde{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 Figure 2 and Figure 3. In Figure 2a we plot the measure of non-Gaussianity $\delta(\tau)$ as a function of time $\tau$ for different values of the coherent state parameter $\mu_{\mathrm{c}}$, over the period $0<\tau<2\,\pi$. The other parameters are set to $\tilde{g}_{0}=1$ and $\mu_{\mathrm{m}}=0$. It is known that the full nonlinear dynamics is periodic (or recurrent) with period $2\,\pi$, see [25], whenever $\tilde{g}^{2}_{0}$ is an integer and this is clearly reflected in our plot.

At $\tau=2\pi$, 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 $\tilde{g}_{0}$. For example, when $\tilde{g}_{0}=0.5$, the cavity state becomes a superposition of coherent states at $\tau=2\pi$, also known as a cat state [22]. However, if $\tilde{g}_{0}^{2}$ is integer, we obtain a phase factor of $e^{2\pi i\tilde{g}_{0}^{2}}=1$ in the optical state, and the optics returns to an initial state as well. This is the case in Figure 2, where $\delta(2\pi)=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 $\tau=\pi$, the measure $\delta$ exhibits a local minimum which grows increasingly narrow with larger $\mu_{\mathrm{c}}$. This is apparent from Figure 2b where we have shown a close-up of $\delta$ around $\tau=\pi$ for increasing values of $\mu_{\mathrm{c}}$, and for $\tilde{g}_{0}=1$ and $\mu_{\mathrm{m}}=0$. The dip occurs because at $\tau=\pi$, we find that $\theta_{a}=-2\pi\tilde{g}^{2}_{0}$ and $F=-2\tilde{g}_{0}$. Thus, for integer $\tilde{g}_{0}^{2}$, we have $\sin^{2}{\theta_{a}/2}=0$ and $\sigma_{11}$ becomes $\sigma_{11}=1+2|\mu_{\mathrm{c}}|^{2}\left(1-e^{-4\tilde{g}_{0}^{2}}\right)$. 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 $\mu_{\mathrm{c}}$ yields a logarithmic increase in $\delta(\tau)$, which is evident from the approximation in Eq. (22). Figure 2a 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 $\tau=2\pi$, 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 $\delta(\tau)$ in Figure 2a closely resembles the linear entropy of the traced-out subsystems as found in  [22, 25].

To better understand the behaviour of $\delta(\tau)$ for small times $\tau$, we plot the behaviour of $\log_{10}{\delta(\tau)}$ for $\tau\ll 1$ for different values of $\mu_{\mathrm{c}}$ in Figure 2c. We note that $\delta(\tau)$ increases quickly at first, but soon tends to a near-constant value. This means that $\delta(\tau)$ grows linearly for an interval of small times, which can be seen as the increasing and decreasing parts in Figure 2a.

Finally, we proceed to examine the scaling behaviour of $\delta(\tau)$ at fixed time $\tau=\pi$. Figure 3a shows a $\log_{10}$–$\log_{10}$ plot of the measure $\delta(\tau)$ as a function of the nonlinear coupling $\tilde{g}_{0}$ for different values of $\mu_{\mathrm{c}}$. As $\tilde{g}_{0}$ increases, the amount of non-Gaussianity first grows linearly in the logarithm, then plateaus as $\tilde{g}_{0}$ increases further. The same behaviour occurs for larger $\mu_{\mathrm{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 $\tilde{g}_{0}$. As such, focusing on increasing the coupling $\tilde{g}_{0}$ will only give marginal returns. Similarly, 3b shows $\log_{10}{\delta(\tau)}$ as a function of increasing $\mu_{\mathrm{c}}$ for various values of $\tilde{g}_{0}$.

5.1 Small coherent state parameters

For a small amplitude coherent state of the optics, with $|\mu_{\mathrm{c}}|^{2}\ll 1$, and with the mechanics in a coherent state, we found in (19) that $\delta(\tau)$ scales with $\sim|F|^{2}|\mu_{\mathrm{c}}|^{2}\ln|\mu_{\mathrm{c}}|$. Given the explicit form of $F$, we see that it scales with $F\propto\tilde{g}_{0}$. Since $\delta(\tau)$ in this regime is proportional to $|F|^{2}$, it follows that $\delta(\tau)$ grows quadratically with the light–matter coupling in this regime.

5.2 Large coherent state parameters

We derived an asymptotic form of $\delta(\tau)$ in (22) for the case $|\mu_{\mathrm{c}}|\gg 1$, which we called $\tilde{\delta}(\tau)$. As argued before, the behaviour of the measure $\delta(\tau)$ in this regime depends crucially on the distance of $\theta_{a}$ from the value $2\pi$. In our present case we have that $\theta_{a}\sim\tau^{3}$ for $\tau\ll 1$ and $\theta_{a}\sim-\tilde{g}_{0}^{2}\tau$ for $\tau\gg 1$. The functions that we decided to ignore (except for $\sigma_{11}$) in the derivation of $\tilde{\delta}(\tau)$ are of the form $f_{|\mu_{\mathrm{c}}|}(\theta)=(1-\exp[-\beta\,|\mu_{\mathrm{c}}|^{2}\sin(\theta_{a}/2)])$ or $f_{|\mu_{\mathrm{c}}|}(\theta)=(1-\exp[-\beta^{\prime}\,|\mu_{\mathrm{c}}|^{2}\sin^{2}(\theta_{a}/2)])$ where $\beta$ and $\beta^{\prime}$ are irrelevant numerical constant of order $1$. We focus on $f_{|\mu_{\mathrm{c}}|}(\theta)$ and note that a similar argument applies for the other function as well. Finally, we ignore the transient regime of $\tau\ll 1$ and focus on times $\theta_{a}\sim-\tilde{g}_{0}^{2}\theta$.

To see how well the asymptotic form $\tilde{\delta}(\tau)$ in  (22) approximates the exact measure, we have plotted both the exact form of $\delta(\tau)$ (solid lines) with the asymptotic form (dashed lines) in Figure 4. We note that, even for $|\mu_{\mathrm{c}}|\sim 1$, the asymptotic measure $\tilde{\delta}(\tau)$ well approximates the exact value of $\delta(\tau)$. In fact, it becomes even more accurate as the optical coherent state parameter $\mu_{\mathrm{c}}$ increases, which is to be expected given the nature of the approximation. The asymptotic form also becomes more accurate once we also increase $\tilde{g}_{0}$, as evident in Figure 4b. For $\tilde{g}_{0}=10^{2}$, the approximation is almost entirely accurate. This occurs because the function $\theta_{a}$ increases with $\tilde{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 $\tau=2\pi$ for integer $\tilde{g}_{0}^{2}$ which we can now address analytically by examining the asymptotic covariance matrix elements in (4.2). We find that $F=0$ for all $\tau=2\pi n$ with integer $n$. This means that $\sigma_{42}=0$ and that $\sigma_{22}=1$. We also find that $\theta_{n}(2\pi)=-4\pi\tilde{g}_{0}^{2}$. Thus, if $\tilde{g}_{0}^{2}$ is integer, we find that $\sin^{2}{\theta_{a}/2}=0$ and the final covariance matrix element is $\sigma_{11}=1$. This results in $\mathbf{\sigma}=\mathrm{diag}(1,1,1,1)$ which corresponds to a coherent state, which is fully Gaussian. As a result, the non-Gaussianity vanishes. When $\tilde{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 $\tau=2\pi$.

5.3 Non-Gaussianity in open systems with constant coupling

Any realistic system will suffer from decoherence. In Figure 5 we have plotted the non-Guassianity $\delta(\tau)$ 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). Figure 5a shows the non-Gaussianity for increasing values of the photon decoherence rate $\bar{\kappa}_{\mathrm{c}}=\kappa_{\mathrm{c}}/\omega_{\mathrm{m}}$ with Lindblad operator $\hat{L}_{\mathrm{c}}=\sqrt{\bar{\kappa}_{\mathrm{c}}}\,\hat{a}$ and values $\mu_{\mathrm{c}}=0.1$, $\tilde{g}_{0}=1$ and $\mu_{\mathrm{m}}=0$. We have chosen a low value of $\mu_{\mathrm{c}}$ to ensure high numerical accuracy of the simulation, as larger values quickly lead to numerical instabilities. We note that the non-Gaussianity $\delta(\tau)$ tends towards a steady value, which is clear from the fact that the higher values of decoherence start to coincide around $\tau=5\pi$. We also note that around $\tau=2\pi 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 5b we have plotted the non-Gaussianity $\delta(\tau)$ for increasing values of phonon decoherence rate $\bar{\kappa}_{\mathrm{m}}$ with Lindblad operator $\hat{L}_{\mathrm{m}}=\sqrt{\bar{\kappa}_{\mathrm{m}}}\,\hat{b}$ and the same values as before.

6 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 $\tilde{g}(\tau)$ is time-dependent as a consequence of, for example, trap instabilities. Time-dependent 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 $\tilde{g}(\tau)$ 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.

6.1 Modelling the trap modulation

We shall model a time-dependent light–matter coupling $\tilde{g}(\tau)$ by assuming that the coupling has the simple form

[TABLE]

Here, $\tilde{g}_{0}$ is the expected value of the coupling, $\epsilon$ is the amplitude of oscillation and $\Omega_{0}:=\omega_{0}/\omega_{\textrm{m}}$ is the dimensionless frequency that determines how the coupling oscillates in time. We can insert this ansatz in the general expressions (2.2) and obtain an explicit form for this case. The full expressions for the coefficients in (2.2) are again very long and cumbersome, and we do not print them here. They are listed in (D).

We can now compute $\delta(\tau)$ for this time-dependent coupling for initial coherent states, and we display the results in Figure 6. In 6a we plot $\delta(\tau)$ vs. $\tau$ for different values of the oscillation frequency $\Omega_{0}$. Note that we here include a larger range of $\tau$ to capture potentially recurring behaviour. In the limit $\Omega_{0}\rightarrow 0$ we recover the time-independent solution, as expected. Interestingly, when $\Omega_{0}\neq 0$ we see that we can achieve higher values for the nonlinear measure $\delta(\tau)$. This is especially pronounced as $\Omega_{0}\rightarrow 1$, where the trap oscillation frequency is equal to the mechanical frequency $\omega_{\mathrm{m}}$, for which $\delta(\tau)$ ceases to oscillate periodically, but instead steadily increases. We discuss this case in detail in the following section.

6.2 Trap modulation on resonance

The functions (D) contain denominators of the form $\Omega_{0}-1$. Therefore, among all possible values of $\Omega_{0}$, we can ask what happens on resonance, i.e., when $\Omega_{0}=1$. Figure 6a already provides evidence that the system should behave markedly differently.

At resonance, where $\Omega_{0}=1$, the functions (D) take the relatively simple form

[TABLE]

We have plotted in Figure 6 the exact measure of non-Gaussianity $\delta(\tau)$ in the resonant case for initially coherent states and for different values of $\mu_{\mathrm{c}}$. As anticipated, here we no longer have recurrent behaviour. Instead, the non-linearity increases as $\ln{\tau}$. Formally, this growth can continue for arbitrarily large times $\tau$, however, the maximum time $\tau$ that can be achieved in practice is limited by the coherence time of the experiment. Similarly, we plotted $\delta(\tau)$ for various values of $\epsilon$ in Figure 6c. We note that $\delta(\tau)$ oscillates increasingly rapidly with larger $\epsilon$ but with decreasing amplitude for increasing $\tau$, as $|F|^{2}\sim\tilde{g}_{0}^{2}\,\epsilon^{2}\,\tau^{2}$ becomes the dominant term for $\tau\gg 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 $\langle\hat{a}^{\dagger}\hat{a}\rangle$ 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 resonance, this linear displacement grows linearly in time. See also [48] for further insight once the rotating wave approximation has been applied.

6.3 Large coherent state parameters at resonance

Using the explicit form of the coefficients (6.2), we note that $|F|^{2}=F_{\hat{N}_{a}\,\hat{B}_{-}}^{2}+F_{\hat{N}_{a}\,\hat{B}_{+}}^{2}$ has the asymptotic behaviour $|F|^{2}\sim\frac{1}{4}\tilde{g}_{0}^{2}\,\epsilon^{2}\,\tau^{2}$ for $\tau\gg 1$. This implies that $\exp[-|F|^{2}]\ll 1$ for large $\tau$ 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 $\delta(\tau)$. This observation allows us to compute the symplectic eigenvalues, which read $\nu_{+}=1+2|\mu_{\mathrm{c}}|^{2}\left(1-e^{-4\,|\mu_{\mathrm{c}}|^{2}\sin^{2}{\theta_{a}/2}}\,e^{-|F|^{2}}\right)$ and $\nu_{-}=\sqrt{1+4\,|\mu_{\mathrm{c}}|^{2}\,|F|^{2}}$, and they match the expressions (4.2). We again stress that we have retained the exact expression for $\sigma_{11}$ to capture some crucial features of the non-Gaussianity, such as $\delta(0)=0$.

In Figure 7, we compare the exact measure $\delta(\tau)$ at resonance with the asymptotic form derived in (22). The solid lines represent the exact measure $\delta(\tau)$ and the dashed lines represent the asymptotic expression. In Figure 7a we compare them for different values of $\mu_{\mathrm{c}}$. We note that, except for at very small $\tau$, the asymptotic form is entirely accurate and gets even more precise for increasing values of $\mu_{\mathrm{c}}$. This is a consequence, as we noted before, of the exponentials in (B.1) that suppress some elements for large $\mu_{\mathrm{c}}$, unless $\theta_{a}=n\,\pi$. Similarly, in Figure 7b we have plotted $\delta(\tau)$ and its asymptotic form for different values of the oscillation amplitude $\epsilon$. Again, the suppression of the exponentials with increasing $\epsilon$ means that larger values of $\epsilon$ yield a more accurate expression.

6.4 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 $\delta$ as a function of time for different values of photon and phonon decoherence. Figure 8a shows the system at resonance with photons leaking from the cavity with a rate $\bar{\kappa}_{\mathrm{c}}=\kappa_{\mathrm{c}}/\omega_{\mathrm{m}}$ for parameters $\mu_{\mathrm{c}}=0.1$, $\tilde{g}_{0}=1$, $\epsilon=0.5$ and $\mu_{\mathrm{m}}=0$. We note that $\bar{\kappa}_{\mathrm{c}}=0.3$ yields what is essentially a steady-state of the non-Gaussianity. In Figure 8b we note the same behaviour but for phonon decoherence with rate $\bar{\kappa}_{\mathrm{m}}=\kappa_{\mathrm{m}}/\omega_{\mathrm{m}}$.

7 Discussion and practical implementations

We have employed a measure of non-Gaussianity $\delta(\tau)$ 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 $\tilde{g}_{0}$, the non-Gaussianity $\delta(\tau)$ scales differently in two contrasting regimes: (i) For a weak optical input coherent state $\lvert\mu_{\mathrm{c}}\rvert$, the nonlinear character of the state grows as $\tilde{g}_{0}^{2}\lvert\mu_{\mathrm{c}}\lvert^{2}\,\ln{\lvert\mu_{\mathrm{c}}\lvert}$ if the mechanics is also in a coherent state, (ii) conversely, for large $\lvert\mu_{\mathrm{c}}\lvert$, the nonlinear character of the state grows logarithmically with the quantity $\tilde{g}_{0}|\mu_{\mathrm{c}}|$, which also holds when the mechanical element is not fully cooled. The same general scaling with $|\mu_{\mathrm{c}}|$ occurs when $\tilde{g}(\tau)$ is time-dependent.

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.

7.1 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 $\mu_{\mathrm{c}}$ to infinity yields $\lim_{\mu_{\mathrm{c}}\rightarrow\infty}\delta=\infty$. 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 $\delta(\tau)$ 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 $\tau=\pi$, which satisfies the occurrence of non-classicality in conjunction with the non-Gaussianity. The state is however fully disentangled at $\tau=2\pi$, 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 $\tau=2\pi 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.

7.2 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 $\omega_{\mathrm{m}}$ and the optical decoherence rate $\kappa_{\mathrm{c}}$. Such experiments usually involve a strong laser drive, which tends to wash out the non-linearity. 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 $\hat{a}$ are treated as flucutations around a strong optical field as such: $\hat{a}\rightarrow\hat{a}=\alpha+\hat{a}^{\prime}$, where $\hat{a}^{\prime}$ 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 $\mu_{\mathrm{c}}$ increases the non-Gaussianity, we cannot generalise this result to the linearised dynamics.

7.3 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 $\tilde{g}(\tau)$ is modulated. An explanation of this phenomena was provided in [48]. Consider the force $\vec{F}$ exerted by the photons on the mechanics. For a number of $n$ photons, this force is proportional to $\vec{F}\propto(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 $\tilde{g}(\tau)$, the photon-pressure force $\vec{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  E, which is based on the work in [72]. 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  E, this phase is denoted by $\varphi(\tau)$ 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 $\varphi(\tau)=\frac{\pi}{2}\left(1+\epsilon\sin{\Omega_{0}\tau}\right)$, with $\Omega_{0}=\omega_{0}/\omega_{\mathrm{m}}$, and where $\omega_{0}$ is the phase modulation frequency, we obtain the expression used in Section 6. If, then, the phase frequency is resonant with $\Omega_{0}=1$, it should be possible to increase the non-Gaussianity even further.

• ii)

Paul traps. The shuttling of ions has been demonstrated [65, 64] 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 $\omega_{\mathrm{m}}$. Since the micromotion moves the bead around the potential minimum with a frequency $\omega_{\mathrm{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 $\omega_{\mathrm{d}}$ equal to $\omega_{\mathrm{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, $\omega_{\mathrm{d}}\leq\omega_{\mathrm{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 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 $\tau$, 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.

7.4 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(\hat{\rho}_{AB})$ of a bipartite state $\hat{\rho}_{AB}$ is bounded by $S(\hat{\rho}_{AB})\geq|S(\hat{\rho}_{A})-S(\hat{\rho}_{B})|$, through the Araki–Lieb inequality [78] where $\hat{\rho}_{A}$ and $\hat{\rho}_{B}$ are the reduced states of the optical and mechanical subsystems, respectively. Therefore, the measure of non-Gaussianity $\delta(\tau)$ that we defined in (15) is lower-bounded by

[TABLE]

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.

8 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 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 $\lvert\mu_{\mathrm{c}}\rvert\geq 1$, the amount of non-Gaussianity grows logarithmically with the input average number of excitations $|\mu_{\mathrm{c}}|$ 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.

Acknowledgments

We thank Fabienne Schneiter, Daniel Braun, Nathanaël Bullier, Antonio Pontin, Peter F. Barker, Ryuji Takagi, Francesco Albarelli, Marco G. Genoni, James Bateman and the reviewers for helpful comments and discussions.

SQ acknowledges support from the EPSRC Centre for Doctoral Training in Delivering Quantum Technologies and thanks the University of Vienna for its hospitality. DR would like to thank the Humboldt Foundation for supporting his work with their Feodor Lynen Research Fellowship. This work was supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 732894 (FET Proactive HOT).

Data availability statement

The figures in this work were generated using Mathematica. The numerical results in Section 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.