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

TL;DR

This paper presents an analytical method for solving the time evolution of nonlinear optomechanical systems with arbitrary parameters, and analyzes how mechanical squeezing influences the system's non-Gaussianity over time.

Contribution

The authors develop a Lie algebra-based approach to solve the dynamics of nonlinear optomechanical systems with time-dependent parameters, enabling analysis of non-Gaussianity evolution.

Findings

Adding mechanical squeezing generally reduces non-Gaussianity.

Sinusoidal squeezing leads to Mathieu equations for system dynamics.

Non-Gaussianity increases with time and squeezing amplitude in certain regimes.

Abstract

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

Equations

$$\hat{H}$$

$$\frac{\hat{H}}{\hslash {\omega }_{\mathrm{m}}}=$$

$$\hat{U}(t)=\overleftarrow{\mathcal{T}}\exp \left[-\frac{i}{\hslash }{\int }_0^{t}d{t}^{\prime }\hat{H}(t^{\prime })\right],$$

$$\hat{U}(t)=\prod _n{\hat{U}}_n(t)=\prod _n\exp [-i{F}_n(t){\hat{G}}_n],$$

$$\mathbf{S}=\prod _{n=1}^{N(2N+1)}{\mathbf{S}}_n,$$

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

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

$${\hat{B}}_{+}^{(2)}$$

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

$$\hat{U}(\tau ):=$$

$$\times e^{-i{F}_{{\hat{N}}_a{\hat{B}}_{-}}{\hat{N}}_a{\hat{B}}_{-}},$$

$${\hat{U}}_{\mathrm{sq}}$$

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

$$-2{\int }_0^{\tau }\mathrm{d}{\tau }^{\prime }\tilde{\mathcal{G}}({\tau }^{\prime })\Im \xi {\int }_0^{{\tau }^{\prime }}\mathrm{d}{\tau }^{\prime \prime }{\tilde{\mathcal{D}}}_1({\tau }^{\prime \prime })\Re \xi ({\tau }^{\prime \prime }),$$

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

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

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

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

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

$$\xi (\tau )=P_{11}(\tau )-i{\int }_0^{\tau }\mathrm{d}{\tau }^{\prime }{P}_{22}(\tau ),$$

$$\alpha (\tau )=$$

$$\beta (\tau )=$$

$${\ddot{P}}_{11}+(1+4{\tilde{\mathcal{D}}}_2(\tau ))P_{11}$$

$${\ddot{P}}_{22}-4\frac{{\dot{\tilde{\mathcal{D}}}}_2(\tau )}{1+4{\tilde{\mathcal{D}}}_2(\tau )}{\dot{P}}_{22}+(1+4{\tilde{\mathcal{D}}}_2(\tau ))P_{22}$$

$${\ddot{I}}_{P_{22}}+(1+4{\tilde{\mathcal{D}}}_2(\tau ))I_{P_{22}}=0,$$

$$I_{P_{22}}={\int }_0^{\tau }\mathrm{d}{\tau }^{\prime }{P}_{22}({\tau }^{\prime }).$$

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

$$|\Psi (\tau )⟩=$$

$$\displaystyle e^{-i\,\left(\Omega_{\textrm{c}}\,\tau+F_{\hat{N}_{a}}+F_{\hat{N}_{a}\,\hat{B}_{+}}\,F_{\hat{B}_{-}}+F_{\hat{N}_{a}\,\hat{B}_{-}}\,F_{\hat{B}_{+}}+\Im\left(K_{\hat{N}_{a}}\,\mu_{\rm{m}}\right)\right)\,n}\,\ket{n}\otimes\ket{\phi_{n},\tilde{\mathcal{D}}_{2}(\tau)}\biggr{]}\,,$$

$${\hat{\rho }}_{\text{Mech.}}(\tau )$$

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

$$\delta (\tau )=S({\hat{\rho }}_{\mathrm{G}}(\tau )),$$

$$|S({\hat{\rho }}_A)-S({\hat{\rho }}_B)|\le S({\hat{\rho }}_{AB})\le S({\hat{\rho }}_A)+S({\hat{\rho }}_B),$$

$${\nu }_{\mathrm{Op}}^2=$$

$$\displaystyle-e^{-4|\mu_{\rm{c}}|^{2}\sin^{2}\theta}\,e^{-4|K_{\hat{N}_{a}}|^{2}}2\,e^{-3|K_{\hat{N}_{a}}|^{2}}\,\Re\left\{e^{i\theta}\,e^{|\mu_{\mathrm{c}}|^{2}(e^{2i\theta}-1)}\,e^{2|\mu_{\rm{c}}|^{2}(e^{-i\theta}-1)}\right\}\biggr{)}\,,$$

$${\nu }_{\mathrm{Me}}^2=$$

$${\nu }_{\mathrm{Op}}<\sqrt{1+4|{\mu }_{\mathrm{c}}{|}^2+4|{\mu }_{\mathrm{c}}{|}^4},$$

$${\nu }_{\mathrm{Op},\max }\sim$$

Full text

1 Introduction

The mathematical understanding of optomechanical systems operating in the nonlinear quantum regime is a major topic of current interest. While most experiments effectively undergo linear dynamics, governed by quadratic Hamiltonians that emerge following a ‘linearisation’ procedure [1, 2, 3], many experiments now operate in the fully nonlinear regime [4, 5, 6] where this procedure fails. It is therefore highly desirable to provide a complete and analytic characterisation of the fully nonlinear system dynamics. Analytic solutions have previously been found for a constant light–matter coupling [7, 8] and, more recently, the time-dependent case was solved [9].

The inherently nonlinear interaction between the optical field and the mechanical element in an optomechanical system allows for the generation of non-Gaussian states. Starting from a broad class of initial states, including coherent states, the vacuum, and thermal states, this is only possible in the nonlinear regime; in contrast, quadratic Hamiltonians take input Gaussian states to output Gaussian states. As such, investigating the non-Gaussianity of optomechanical states can only be performed once the time-evolution in the nonlinear regime has been solved, which is the primary aim of this work. Interestingly, a number of non-classical and non-Gaussian states have been found to constitute an important resource for sensing. Schrödinger cat states [8, 7], compass states [10, 11] and hypercube states [12] – which are all non-Gaussian and highly non-classical states – have all been found to have applications for sensing. More generally, the detection and generation of non-Gaussianity in optomechanical systems has been extensively studied in theoretical proposals [13, 14, 15] as well as in experiments [4, 16, 5]. Beyond optomechanics, the presence of a nonlinear element is also key to a number of quantum information tasks, such as obtaining a universal gate set for quantum computing [17, 18], teleportation [19], distillation of entanglement [20, 21], error correction [22], and non-Gaussianity has been explored as the basis of an operational resource theory [23, 24, 25].

Optomechanical systems offer a natural nonlinear coupling which, if strong enough, may lead to substantial non-Gaussianity in the evolved state. It is therefore essential to better understand the dynamics of such systems, with special emphasis on the interplay between nonlinearities and other Hamiltonian terms in this dynamics. An important question to be answered is thus how do the different aspects of an optomechanical system affect the non-Gaussianity of the state at a given time? A preliminary study of non-Gaussianity in standard optomechanical systems provided the first tools to approach this question [9], however, optomechanical systems can exhibit additional, potentially more interesting, effects. An important non-classical effect that can be included into optomechanical systems is squeezing of the optical or mechanical modes. The addition of squeezing has been found to be beneficial for sensing since it increases the sensitivity in a specific field quadrature. For example, it has been shown that squeezed light injected into LIGO can be used to enhance the detection of gravitational waves [26]. Similarly, mechanical squeezing can aid the amplification and measurement of weak mechanical signals [27].

In this work, we study the non-Gaussianity of a quantum system of two bosonic modes characterised by an optomechanical Hamiltonian with the standard cubic light–matter interaction term, and with the addition of a mechanical displacement term and a mechanical squeezing term. We extend a recently developed solution of the time evolution operator induced by a plain optomechanical Hamiltonian [9] to include the additional terms of interest here. Interestingly, for time-dependent squeezing modulated at resonance, we find that the dynamics are governed by the well-studied Mathieu equation. We subsequently derive perturbaive solutions and show that these coincide with the physically intuitive rotating-wave approximation for large times. The decoupling methods used in this work have a long tradition in quantum theory [28, 29, 30] and were recently applied to problems such as the one at hand [31, 9, 32]. We use the resulting analytic solutions to compute the amount of non-Gaussianity of the state using a measure of relative entropy [33, 34] for both a constant and a time-dependent mechanical squeezing parameter.

Our results indicate that the non-Gaussian character of an initially coherent state decreases in general with an increasing squeezing parameter. However, when the squeezing is applied periodically at mechanical resonance, the non-Gaussianity increases approximately linearly with time and the amplitude of the squeezing. The competition between the amount of squeezing and the strength of the nonlinear term is difficult to compute explicitly; instead, we provide asymptotic expressions in terms of upper and lower bounds to the non-Gaussianity in different regimes. A conclusive answer requires further investigation, potentially providing a concise expression where such competition can be easily understood.

The paper is structured as follows. In Section 2, we introduce the nonlinear Hamiltonian with mechanical squeezing. This is followed by Section 3 where we provide a short introduction to the methods used to solve the dynamics. The full derivation can be found in B. Following this, we review the measure of non-Gaussianity and derive expressions for an asymptotic expression and a reduced measure in Section 4. In Section 5, we then specialise to two specific cases and compute the amount of non-Guassianity for constant squeezing (Section 5.2), and for modulated squeezing (Section 5.3). Finally, we conclude with a discussion in Section 6 and some final remarks in Section 7.

2 Dynamics

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

2.1 Hamiltonian

In this work we consider the two-mode Hamiltonian

[TABLE]

where $\hat{H}_{0}:=\hbar\,\omega_{\mathrm{c}}\hat{a}^{\dagger}\hat{a}+\hbar\,\omega_{\mathrm{m}}\,\hat{b}{}^{\dagger}\hat{b}$ is the free Hamiltonian, while $\omega_{\mathrm{c}}$ and $\omega_{\mathrm{m}}$ are the frequencies of the cavity mode and the mechanical resonator respectively.

The Hamiltonian (1) describes the dynamics of a number of different systems. For example, $\mathcal{G}(t)$ appears in optomechanical systems as a standard coupling term due to radiation pressure obtained for Fabry-Pérot cavities, where one end of the cavity is a mirror that can move freely [35]. Such coupling appears also within systems with a central translucent membrane in a rigid optical cavity [36], levitated nanodiamonds [37] or optomechanical crystals [38, 39]. A depiction of a levitated nanosphere interacting with cavity modes can be found in Figure 1.

The Hamiltonian (1) reduces to the standard optomechanical Hamiltonian when $\mathcal{D}_{1}=\mathcal{D}_{2}=0$.222We note here that non-Gaussianity for the case $\mathcal{D}_{1}=\mathcal{D}_{2}=0$ has been already studied [9]. The term weighted by the coupling $\mathcal{D}_{1}$ corresponds to an externally imposed displacement of the mechanical part, which can be induced by a piezoelectric element connected to its support or by an external acceleration, such as that caused by the gravitational force acting on the mechanical element [40, 41]. The term governed by $\mathcal{D}_{2}$ can be thought of as a modulation of the trap frequency and leads to squeezing of the mechanics., which can be externally imposed employing another strong optical field or an electrostatic force [42].

2.2 Dimensionless dynamics

To understand which dimensionless parameters are relevant to the dynamics of the system, we start by introducing dimensionless quantities and rescaling the Hamiltonian. Such a rescaling also serves to simplify the notation and any graphical representation of the system dynamics. We achieve this by dividing the functions in the Hamiltonian by the mechanical frequency $\omega_{\mathrm{m}}$. The action corresponds to switching from the laboratory time $t$ to $\tau=\omega_{\mathrm{m}}\,t$, where $\tau$ is the new, dimensionless time. The optical frequency becomes $\Omega_{\mathrm{c}}:=\omega_{\mathrm{c}}/\omega_{\mathrm{m}}$. In addition, the couplings in the Hamiltonian become $\tilde{\mathcal{G}}(\tau)=\mathcal{G}(\omega_{\mathrm{m}}\,t)/\omega_{\mathrm{m}}$, $\tilde{\mathcal{D}}_{1}(\tau)=\mathcal{D}_{1}(\omega_{\mathrm{m}}\,t)/\omega_{\mathrm{m}}$ and $\tilde{\mathcal{D}}_{2}(\tau)=\mathcal{D}_{2}(\omega_{\mathrm{m}}\,t)/\omega_{\mathrm{m}}$. We also rescale the Hamiltonian by $\hbar$, meaning that $\hat{H}$ becomes

[TABLE]

which is the Hamiltonian that we will be working with.

3 Solving the dynamics

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

3.1 Continuous variables and covariance-matrix formalism

When solving the dynamics, we employ methods from the continuous variable formalism [43, 3]. Specifically, the methods are used to solve the time-evolution of the quadratic part of the system, and to describe its action on the nonlinear light–matter interaction term. We briefly review the continuous variable formalism here.

In recent years, thanks to progress in the mathematical framework provided by the covariance matrix formalism [43, 3], it has become clear that Gaussian states constitute an extremely valuable toolkit to investigate quantum information processing in quantum setups, and in relativistic ones as well [44]. The main advantage is that the covariance matrix formalism provides a powerful set of mathematical tools to treat Gaussian states of bosonic fields that undergo linear transformations of the creation and annihilation operators fully analytically [43]. Ultimately, Gaussian states are the paramount resource for continuous variables quantum information processing and computation [17] and have become a standard feature in most quantum optics laboratories. However, it should also be pointed out that these methods can be used to describe the evolution of operators in the Heisenberg picture, even when the states considered are not Gaussian.

In quantum mechanics, the initial state $\hat{\rho}_{i}$ of a system of $N$ bosonic modes with operators $\{\hat{a}_{n},\hat{a}^{{\dagger}}_{n}\}$ evolves to a final state $\hat{\rho}_{f}$ through the standard Schrödinger equation $\hat{\rho}_{f}=\hat{U}\,\hat{\rho}_{i}\,\hat{U}^{{\dagger}}$, where $\hat{U}$ implements the transformation of interest, such as time evolution. If the state $\hat{\rho}$ is Gaussian and the Hamiltonian $\hat{H}$ is quadratic in the operators, it is convenient to introduce the vector $\hat{\mathbb{X}}=(\hat{a}_{1},\ldots,\hat{a}_{N},\hat{a}^{{\dagger}}_{1},\ldots,\hat{a}^{{\dagger}}_{N})^{\rm{T}}$, where $\rm{T}$ denotes the transpose of the vector. Similarly, the vector of first moments $d:=\langle\hat{\mathbb{X}}\rangle$ and the covariance matrix $\boldsymbol{\sigma}$ are defined by $\sigma_{nm}:=\langle\{\hat{X}_{n},\hat{X}^{{\dagger}}_{m}\}\rangle-2\langle\hat{X}_{n}\rangle\langle\hat{X}_{m}^{{\dagger}}\rangle$, where $\{\cdot,\cdot\}$ stands for anticommutator and all expectation values of an operator $\hat{\mathcal{A}}$ are defined by $\langle\hat{\mathcal{A}}\rangle:=\text{Tr}(\hat{\mathcal{A}}\,\hat{\rho})$.

In this language, the canonical commutation relations read $[\hat{X}_{n},\hat{X}_{m}^{{\dagger}}]=i\,\Omega_{nm}$, where the $2N\times 2N$ matrix $\boldsymbol{\Omega}$ is known as the symplectic form [43]. We then notice that, while arbitrary states of bosonic modes are, in general, characterised by infinite real parameters, a Gaussian state is uniquely determined by its first and second moments, $d_{n}$ and $\sigma_{nm}$ respectively [43]. Furthermore, unitary transformations quadratic in the annihilation and creation operators, such as Bogoliubov transformations [45], preserve the Gaussian character of a Gaussian state and can always be represented by a $2N\times 2N$ symplectic matrix $\boldsymbol{S}$ that preserves the symplectic form, i.e., $\boldsymbol{S}^{{\dagger}}\,\boldsymbol{\Omega}\,\boldsymbol{S}=\boldsymbol{S}\,\boldsymbol{\Omega}\,\boldsymbol{S}^{{\dagger}}=\boldsymbol{\Omega}$.333Note that $\boldsymbol{\sigma}$ is complex in our choice of basis, which implies taking the Hermitian conjugate of $\boldsymbol{S}$. In a similar manner, the symplectic matrix $\boldsymbol{S}$ that encodes the evolution of a state is generated by the Hamiltonian matrix $\boldsymbol{H}$, which is defined by $\hat{H}=\hat{\mathbb{X}}^{\dagger}\boldsymbol{H}\hat{\mathbb{X}}/2+\hat{\mathbb{X}}^{\dagger}d$. The symplectic matrix becomes $\boldsymbol{S}=\mathrm{exp}\left[\boldsymbol{\Omega}\boldsymbol{H}\right]$.

The Schrödinger equation can be translated in this language to the simple equation $\boldsymbol{\sigma}_{f}=\boldsymbol{S}\,\boldsymbol{\sigma}_{i}\,\boldsymbol{S}^{{\dagger}}$ for the second moments, and $\boldsymbol{r}_{f}=\boldsymbol{S}\,\boldsymbol{r}_{i}$ for the first moments, which shifts the problem of usually untreatable operator algebra to simple $2N\times 2N$ matrix multiplication. Here, the indices $i$ and $f$ denote the initial and final state, respectively. Finally, Williamson’s theorem guarantees that any $2N\times 2N$ Hermitian matrix, such as the covariance matrix $\boldsymbol{\sigma}$, can be decomposed as $\boldsymbol{\sigma}=\boldsymbol{S}^{{\dagger}}\,\boldsymbol{\nu}_{\oplus}\,\boldsymbol{S}$, where $\boldsymbol{S}$ is an appropriate symplectic matrix. The diagonal matrix $\boldsymbol{\nu}_{\oplus}=\textrm{diag}(\nu_{1},\dots,\nu_{N},\nu_{1},\dots,\nu_{N})$ is known as the Williamson form of the state and $\nu_{n}:=\coth(\frac{\hbar\,\omega_{n}}{2\,k_{B}\,T})\geq 1$ (where we have introduced normal frequencies $\omega_{n}$ and a nominal temperature $T$) are the symplectic eigenvalues of the state [46].

Williamson’s form $\boldsymbol{\nu}_{\oplus}$ contains information about the local and global mixedness of the state of the system [43]. The state is pure if $\nu_{n}=1$ for all $n$ and is mixed otherwise. As an example, the thermal state $\boldsymbol{\sigma}_{th}$ of a $N$-mode bosonic system is simply given by its Williamson form, i.e., $\boldsymbol{\sigma}_{th}=\boldsymbol{\nu}_{\oplus}$.

3.2 Decoupling of a time evolution operator

The time evolution of a system with time-dependent Hamiltonian $\hat{H}(t)$ is

[TABLE]

where $\overset{\leftarrow}{\mathcal{T}}$ is the time ordering operator. This expression simplifies dramatically when the Hamiltonian $\hat{H}$ is time independent, in which case one obtains $\hat{U}(t)=\exp[-\frac{i}{\hbar}\,\hat{H}\,t]$ as a solution to the time-dependent Schrödinger equation. We are, however, interested in Hamiltonians with time-dependent parameters. Any Hamiltonian can be cast in the form $\hat{H}=\sum_{n}\hbar\,g_{n}(t)\,\hat{G}_{n}$, where the $\hat{G}_{n}$ are time independent, Hermitian operators and the $g_{n}(t)$ are time-dependent real functions. The choice of $\hat{G}_{n}$ need not be unique, and if this is the case, a specific choice is motivated by the specific aims of the problem.

We say that the time evolution operator (3) has been decoupled if it can be written as [28, 29]

[TABLE]

where the real functions $F_{n}(t)$ are in general time-dependent. It has been shown that these functions can be found as solutions to a set of differential equations and are determined solely by the parameters $g_{n}(t)$ of the Hamiltonian [28]. The order of the operators in (4) is not unique; a different order changes the form of the functions $F_{n}(t)$, but the not the expectation value of physical quantities. A more detailed outline of these decoupling techniques may be found in A.

It is possible to obtain an even more explicit decoupling (4) in the context of Gaussian states and linear (i.e., quadratic in the operators) interactions. Given a set of $N$ bosonic modes, there are $N\,(2\,N+1)$ independent quadratic Hermitian operators, which we can denote $\hat{G}_{n}$, that can be formed by arbitrary quadratic combinations of the creation and annihilation operators [47].444For example, $\hat{G}_{1}=\hat{a}_{1}^{{\dagger}}\hat{a}_{1}+\hat{a}_{1}\hat{a}_{1}^{{\dagger}}$ or $\hat{G}_{8}=\hat{a}_{2}^{{\dagger}}\hat{a}_{5}^{{\dagger}}+\hat{a}_{5}\hat{a}_{2}$, where the numbering and ordering of the generators $\hat{G}_{n}$ is a matter of convenience. Work in this direction has also been done in [48] We also recall that any unitary transformation induced by a quadratic operator, including the quadratic time evolution operator (3), can be represented by a $2N\times 2N$ symplectic matrix $\boldsymbol{S}$. Combining all of this together, it can be shown [47] that the symplectic matrix $\boldsymbol{S}$ that represents the time evolution operator (3) takes the form

[TABLE]

where the symplectic matrices $\boldsymbol{S}_{n}$ are given by $\boldsymbol{S}_{n}:=\exp[F_{n}(t)\,\boldsymbol{\Omega}\,\boldsymbol{G}_{n}]$ and the matrices $\boldsymbol{G}_{n}$ can be obtained through $\hat{G}_{n}=\frac{1}{2}\,\mathbb{X}^{{\dagger}}\,\boldsymbol{G}_{n}\,\mathbb{X}$, with the restriction that the generator matrix $\boldsymbol{G}_{n}$ must be Hermitian. The techniques to obtain the real, time-dependent functions $F_{n}(t)$ are the same as in the more general case described above. More details can be found in A.

3.3 Decoupling algebra of the nonlinear Hamiltonian

Decoupling of the Hamiltonian (1) can be done using different choices of the Hermitian operators $\hat{G}_{n}$. Here, we find it convenient to consider the closed finite $9$-dimensional Lie algebra generated by the following set of Hermitian basis operators

[TABLE]

which form the smallest set of operators in the Lie algebra that generate the Hamiltonian (2).555In particular, the Hamiltonian (2) is generated by a linear combination of the Hermitian operators $\hat{N}_{a}$, $\hat{N}_{b}$, $\hat{B}_{+}$, $\hat{B}^{(2)}_{+}$ and $\hat{N}_{a}\,\hat{B}_{+}$.

A generic time evolution operator $\hat{U}(\tau)$ induced by an arbitrary Hamiltonian cannot in general be written in the form (4) for finite number of operators $\hat{G}_{n}$. A finite decoupling (4) is however possible when the operators forms a finite Lie algebra that is closed under commutation. This is the case for the Hamiltonian in (1), since the commutator of any two elements in the algebra (3.3) yields a linear combination of the elements of the algebra. This allows us to make an informed ansatz for the evolution operator as we will see below.

3.4 Decoupling of a nonlinear time-dependent optomechanical Hamiltonian

In order to achieve the main aim of this work, we need an analytical expression of the decoupling (4) given our Hamiltonian (1). While we will show that we can always obtain a formal expression for the evolution, the coefficients that make up the evolution cannot always be computed analytically, as will be clear for certain choices of the mechanical squeezing function $\tilde{D}_{2}(\tau)$.

We find it convenient to proceed by collecting all quadratic terms – including the squeezing term with $\tilde{\mathcal{D}}_{2}$ in (2) – as a separate operator which we call $\hat{U}_{\mathrm{sq}}(\tau)$. This choice allows us to study the action of the quadratic and nonlinear parts separately, which can be solved through different means. Since we are interested in computing the first and second moments of the system for the purpose of computing the non-Gaussianity, it is straight-forward to apply $\hat{U}_{\mathrm{sq}}(\tau)$ to the operator $\hat{b}$ as a symplectic transformation.

We now make an ansatz for the time-evolution operator $\hat{U}(\tau)$ as a finite product of the operators in the algebra:

[TABLE]

where we have defined an evolution operator $\hat{U}_{\rm{sq}}$ as a quadratic evolution operator of the mechanical degree of freedom:

[TABLE]

Here, we have effectively divided the Hamiltonian into a quadratic contribution $\hat{U}_{\rm{sq}}(\tau)$ and a remaining nonlinear contribution with the addition of linear term proportional to $\hat{B}_{\pm}$.

The coefficients in the decoupling above can now be obtained in terms of definite integrals. The full calculations can be found in B. We obtain

[TABLE]

where we have introduced the function

[TABLE]

and where $P_{11}(\tau)$ and $P_{22}(\tau)$ are defined below.

The only problem that we encounter is computing a decoupled form of $\hat{U}_{\rm{sq}}$ in (8). In fact, it has been shown that decoupling of the evolution operator does not yield analytical results except in very specific cases [49]. For our purposes, this is not problematic, because we can calculate the action of $\hat{U}_{\rm{sq}}$ on the first and second moments analytically using the covariance matrix formalism.

3.5 Action of the single-mode squeezing component

Although it is not possible to obtain an analytical decoupling of (8), it is possible to obtain an expression for its action on the operators $\hat{b}$ and $\hat{b}^{\dagger}$. First of all, we note that a Bogoliubov transformation of a single mode operator always has the general expression $\hat{U}_{\mathrm{sq}}^{\dagger}\,\hat{b}\,\hat{U}_{\mathrm{sq}}=\alpha(\tau)\,\hat{b}+\beta(\tau)\,\hat{b}^{\dagger}$, see [49]. The challenge is to find an explicit expression for the Bogoliubov coefficients $\alpha(\tau)$ and $\beta(\tau)$, which satisfy the only nontrivial Bogoliubov identity $|\alpha(\tau)|^{2}-|\beta(\tau)|^{2}=1$. In B we show that the Bogoliubov coefficients $\alpha(\tau)$ and $\beta(\tau)$ can be obtained through

[TABLE]

whose explicit form can be obtained once an explicit expression of the functions $P_{11}(\tau)$ and $P_{22}(\tau)$ is found. Given the previously defined function $\xi(\tau)$ in (10), we also find $\alpha(\tau)=\frac{1}{2}(\xi(\tau)+i\dot{\xi}(\tau))$ and $\beta(\tau)=\frac{1}{2}(\xi^{*}(\tau)+i\dot{\xi}^{*}(\tau))$, where dotted functions imply differentiation with respect to $\tau$.

The two functions $P_{11}$ and $P_{22}$ are determined by the two following uncoupled differential equations:

[TABLE]

where the dot stands for a derivative with respect to $\tau$ and the initial conditions are $P_{11}(0)=P_{22}(0)=1$ and $\dot{P}_{11}(0)=\dot{P}_{22}(0)=0$. Furthermore, the second equation in (3.5) can be written as

[TABLE]

which now has boundary conditions $I_{P_{22}}(0)=0$ and $\dot{I}_{P_{22}}=1$, and where

[TABLE]

The solutions to $P_{11}$ and $P_{22}$ (or $I_{P_{22}}$) can then be used in the expressions (3.4), (10) and (3.5) to find the full dynamics of the state. While the solutions must in general be obtained numerically, we anticipate that there are scenarios, such as constant $\tilde{\mathcal{D}}_{2}$, where the equations above admit analytical solutions.

3.6 Initial state

In this work, we assume that both the optical and mechanical modes are initially in a coherent states, namely $\ket{\mu_{\mathrm{c}}}$ and $\ket{\mu_{\mathrm{m}}}$ respectively, defined as the eigenstates of the annihilation operators, i.e., $\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 optical fields, this is generally a good assumption. On the other hand, within optomechanical systems the mechanical element is typically found initially in a thermal state. Our choice of initial coherent state can be generalised to that of a thermal state in a straight-forward manner, that is, by integrating over the coherent state parameter with an appropriate kernel (as any thermal state may be written as Gaussian average of coherent states, as per its P-representation). Restricting ourselves hence to a single coherent state also for the mechanical oscillator, the initial state $|\Psi(0)\rangle$ reads

[TABLE]

We now proceed to apply (3.4) to this state.

3.7 Full state evolution for general dynamics

For completeness, we present here the full state derived under the evolution with $\hat{U}(\tau)$ for two initially coherent states (15):

[TABLE]

where we have defined $K:=F_{\hat{B}_{-}}+iF_{\hat{B}_{+}}$ and $K_{\hat{N}_{a}}:=F_{\hat{N}_{a}\,\hat{B}_{-}}+iF_{\hat{N}_{a}\,\hat{B}_{+}}$, where $|\phi_{n}(\tau),\tilde{\mathcal{D}}_{2}(\tau)\rangle$ is a mechanical coherent squeezed state where $|\phi_{n}(\tau),\tilde{\mathcal{D}}_{2}(\tau)\rangle=\hat{U}_{\mathrm{sq}}(\tau)|\phi_{n}(\tau)\rangle$, and where $|\phi_{n}(\tau)\rangle$ is a coherent state with $\phi_{n}(\tau):=K^{*}+n\,K_{\hat{N}_{a}}^{*}+\mu_{\mathrm{m}}$. Note that, in the above, we have kept the dependence on $\tau$ implicit but, in general, all exponentials will oscillate in time. We also note that the state (3.7) contains all terms that have been considered in the literature before, including the contributions from a constant nonlinear light–matter term [7], a time-dependent light–matter term [9], and a linear, mechanical displacement term [40]. The main addition here is $\hat{U}_{\mathrm{sq}}$, which takes on the trivial form $\hat{U}_{\rm{sq}}=e^{-i\tau\hat{b}^{\dagger}\hat{b}}$ only if $\tilde{\mathcal{D}}_{2}\neq 0$.

We note here that the expression of (3.7) allows us to compute the reduced state of the mechanics $\hat{\rho}_{\text{Mech.}}(\tau)$ at any time $\tau$, which reads

[TABLE]

We are now ready to consider the non-Gaussianity of the evolved state.

4 Measures of deviation from Gaussianity

The time evolution (3.4) is not linear. Therefore, an initial Gaussian state will evolve, in general, to a non-Gaussian state. Here we ask: is it possible to quantify the deviation from Gaussianity of the state evolving from an initial Gaussian state?

To answer this question we need to find one or more suitable measures of deviation from Gaussianity. In this work we choose to employ a measure based on the comparison between the entropy of the final state and that of a reference Gaussian state [33]. This measure can be understood simply as follows: Let us assume that our initial state $\hat{\rho}(0)$ evolves into the state $\hat{\rho}(\tau)$ at time $\tau$. We can analytically compute the first and second moments of $\hat{\rho}(\tau)$. We then consider a Gaussian state, which we call $\hat{\rho}_{\mathrm{G}}(\tau)$, with the same first and second moments as $\hat{\rho}(\tau)$. It has been shown that this reference state $\hat{\rho}_{\rm{G}}(\tau)$ is indeed the Gaussian state that is closest to $\hat{\rho}(\tau)$ [34]. In general, since $\hat{\rho}(\tau)$ is not determined uniquely by its first and second moments, as is the case for Gaussian states, the two states do not coincide, i.e., $\hat{\rho}(\tau)\neq\hat{\rho}_{\mathrm{G}}(\tau)$.

One way to quantify the difference between two states $\hat{\rho}$ and $\hat{\rho}_{\rm{G}}$ is via the relative entropy $S(\hat{\rho},\hat{\rho}_{\rm{G}})$. It has been shown that the relative entropy $S(\hat{\rho}(\tau),\hat{\rho}_{\mathrm{G}}(\tau))$ is equivalent to the difference between the local von-Neumann entropies of the states [33]. The measure of non-Gaussianity $\delta(\tau)$ can therefore be defined as

[TABLE]

where $S(\hat{\rho})$ is the usual von Neumann entropy of a state $\hat{\rho}$, defined by $S(\hat{\rho}):=-\textrm{Tr}(\hat{\rho}\,\ln\hat{\rho})$.

Since the reference state $\hat{\rho}_{\rm{G}}(\tau)$ is Gaussian, it is fully characterised by its first and second moments. We therefore turn to the continuous variable formalism and consider the covariance matrix $\boldsymbol{\sigma}$ of $\hat{\rho}_{\rm{G}}(\tau)$. Furthermore, we can define the von Neumann entropy $S(\boldsymbol{\sigma})$ of the state as given by $S(\boldsymbol{\sigma})=\sum_{j}s_{V}(\nu_{j})$, where $j$ runs over all the modes, $\nu_{j}$ are the symplectic eigenvalues of $\boldsymbol{\sigma}$ and $s_{V}(\nu_{j})$ is the binary entropy of the state 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)$. The symplectic eigenvalues are defined as $\nu_{j}=|\lambda_{j}|$, where $\lambda_{j}$ are the eigenvalues of the matrix $i\,\boldsymbol{\Omega}\boldsymbol{\sigma}$, where $\boldsymbol{\Omega}$ is the $4\times 4$ symplectic form. Note that, for all physical states, the eigenvalues satisfy $1\leq\nu_{j}$. It follows from the above that a state is non-Gaussian at time $\tau$ if and only if $\delta(\tau)\neq 0$.

An alternative interpretation of this measure is as a quantification of the impurity of $\hat{\rho}_{\mathrm{G}}(\tau)$. While the initial state $\hat{\rho}(\tau)$ remains pure throughout the evolution, such that $S(\hat{\rho}(0))=S(\hat{\rho}(\tau))=0$, the constructed Gaussian reference state $\hat{\rho}_{\mathrm{G}}(\tau)$ does not remain pure. This is not due to external noise, but occurs because we are, loosely speaking, ‘approximating’ the actual state with the Gaussian subset of states.

In this work, we consider unitary dynamics only. If the initial state $\hat{\rho}(0)$ is pure at $\tau=0$, it stays pure throughout its evolution, and the measure thus reduces to

[TABLE]

where $\hat{\rho}_{\mathrm{G}}(\tau)$ is the Gaussian reference state constructed form the first and second moments of $\hat{\rho}(\tau)$. Our challenge is therefore to compute the symplectic eigenvalues $\nu_{j}$ in order to be able to find the expression of $\hat{\rho}_{\mathrm{G}}$. Using the expression for the decoupled time evolution operator (8), we can obtain all of the elements of $\boldsymbol{\sigma}$. These expressions are cumbersome and can be found in C. The expression for the symplectic eigenvalues are too involved and we choose not to print them.

Before we proceed, we also consider the effect of mechanical squeezing on the symplectic eigenvalues. In the continuous variable formalism, a squeezing operation can be represented as a symplectic transformation $\boldsymbol{S}$ acting on the covariance matrix $\boldsymbol{\sigma}$ through congruence: $\boldsymbol{S}\,\boldsymbol{\sigma}\,\boldsymbol{S}^{\dagger}$. All symplectic transformations leave the symplectic eigenvalues $\nu_{j}$ invariant when acted upon in this way. In this work, however, we consider the inclusion of mechanical squeezing as a term in the Hamiltonian, which acts on the fully non-Gaussian state $\hat{\rho}(\tau)$. The presence of the nonlinearity means that the squeezing term acts non-trivially on the full state and can actually affect the symplectic eigenvalues of the Gaussian reference state. The mechanical squeezing parameter $\tilde{D}_{2}(\tau)$ affects all $F$-coefficients, meaning that not only the mechanical subsystem but also the optical subsystem will be affected.

5 Application: Non-Gaussianity for optomechanical systems

In this section, we demonstrate the applicability of our techniques by computing the non-Gaussianity of an optomechanical system. The solutions allow us to consider both constant and time-dependent light–matter couplings, however, in order to obtain explicit results we choose to set $\tilde{\mathcal{G}}(\tau)\equiv\tilde{g}_{0}$ constant throughout this work and refer the reader to [9] for a thorough analysis of the non-Gaussianity of the optomechanical state given a time-dependent light–matter coupling. Furthermore, we set $\tilde{\mathcal{D}}_{1}=0$ throughout the remainder of this work. Since the second moments are not affected by a displacement term, the non-Gaussian character of a state remains unchanged[3].

We consider two cases in this section: one where we assume that the mechanical squeezing parameter is constant, and one where the mechanical squeezing is periodic. Our goal is to derive some general bounds on the non-Gaussianity of the state for each case. First, however, we provide some bounds on the total amount of non-Gaussianity.

5.1 Bounding the full measure

The exact expression for $\delta(\tau)$ is long and cumbersome due to the complex expressions of the covariance matrix elements (C). We therefore provide bounds to the measure that can be expressed as simple analytic functions. Since the full measure $\delta(\tau)$ is an entropy, it can be bounded from above and below by the means of the Araki–Lieb inequality [50], which reads

[TABLE]

where $\hat{\rho}_{AB}$ is the full bipartite state and $\hat{\rho}_{A}$ and $\hat{\rho}_{B}$ are the traced-out subsystems. This inequality allows us to bound the behaviour of the full measure $\delta(\tau)$ in terms of the subsystem entropies. We therefore proceed to define the lower and upper bounds as $\delta_{\rm{min}}(\tau):=|S(\hat{\rho}_{A})-S(\hat{\rho}_{B})|$ and $\delta_{\rm{max}}(\tau):=S(\hat{\rho}_{A})+S(\hat{\rho}_{B})$.

In our case, the subsystems are the traced out optical state $\hat{\rho}_{\rm{Op}}$ and the traced out mechanical state $\hat{\rho}_{\rm{Me}}$. To quantify the entropy of the subsystems, we must find the symplectic eigenvalues of the optical and mechanical subsystems, which we call $\nu_{\rm{Op}}$ and $\nu_{\rm{Me}}$ respectively. Lengthy algebra (see C), the use of the Bogoliubov identities $|\alpha|^{2}=1+|\beta|^{2}$ and $\alpha\,\beta^{*}=\alpha^{*}\,\beta$, and observing that $|E_{\hat{B}_{+}\hat{B}_{-}}|^{2}=e^{-|K_{\hat{N}_{a}}|^{2}}$ (see C for a definition of $E_{\hat{B}_{+}\hat{B}_{-}}$ and its appearance in the first and second moments) allow us to find

[TABLE]

where we recall that $K_{\hat{N}_{a}}:=F_{\hat{N}_{a}\,\hat{B}_{-}}+i\,F_{\hat{N}_{a}\,\hat{B}_{+}}$, and where we have defined $\theta(\tau)=2\left(F_{\hat{N}_{a}^{2}}+F_{\hat{N}_{a}\hat{B}_{+}}F_{\hat{N}_{a}\hat{B}_{-}}\right)$.

The optical symplectic eigenvalue (5.1) is bounded by

[TABLE]

which can be inferred by noting that $K_{\hat{N}_{a}}$ is generally given by an oscillating function multiplied by the strength of the optomechanical coupling $\tilde{g}_{0}$. For specific $\tau$ which ensures that $|K_{\hat{N}_{a}}|^{2}\neq 0$, and then considering $\tilde{g}_{0}\gg 1$, the exponentials in $\nu_{\rm{Op}}$ in (5.1) are suppressed, which means we are left with $\nu_{\rm{Op}}\sim\sqrt{1+4\,|\mu_{\mathrm{c}}|^{2}+4|\mu_{\mathrm{c}}|^{4}}$.

When $S(\hat{\rho}_{\rm{Op}})\gg S(\hat{\rho}_{\rm{Me}})$ or $S(\hat{\rho}_{\rm{Op}})\ll S(\hat{\rho}_{\rm{Me}})$, the bipartite entropy of the Gaussian reference state $S(\hat{\rho}_{\rm{G}})$ is approximately equal to one of the subsystem entropies. To determine when this is the case, we consider the maximum values of $\nu_{\rm{Op}}$ and $\nu_{\rm{Me}}$. In general, when $|\mu_{\rm{c}}|^{2}\gg 1$, and when $|K_{\hat{N}_{a}}|^{2}\gg 1$, which requires $\tilde{g}_{0}\gg 1$ and specific values of $\tau$, the eigenvalues $\nu_{\mathrm{Op}}$ and $\nu_{\mathrm{Me}}$ tend to their maximum values $\nu_{\rm{Op},\rm{max}}$ and $\nu_{\rm{Me},\rm{max}}$, which are

[TABLE]

We note that there are three distinct scenarios which arise from the comparison of the coherent state parameter $|\mu_{\rm{c}}|^{2}$ and the function $K_{\hat{N}_{a}}$:

• i)

First, we assume that $1\ll|K_{\hat{N}_{a}}|\ll 2|\mu_{\mathrm{c}}|$, which implies $\delta(\tau)\sim S(\hat{\rho}_{\rm{Op}})=s_{V}(\nu_{\rm{Op}})$. Here, the non-Gaussianity is well-approximated by

[TABLE]

• ii)

Secondly, we assume that $1\ll 2|\mu_{\rm{c}}|\ll|K_{\hat{N}_{a}}|$, which implies that $\delta(\tau)\sim S(\nu_{\rm{Me}})=s_{V}(\nu_{\rm{Me}})$. Thus we find that

[TABLE]

• iii)

Finally, when $|K_{\hat{N}_{a}}|\sim 2|\mu_{\mathrm{c}}|$ and $|\mu_{\mathrm{c}}|\gg 1$, we have $S(\hat{\rho}_{A})\sim S(\hat{\rho}_{B})$. In this case, the Araki–Lieb bound is not very informative since the left-hand-side is zero and must evaluate the non-Gaussianity exactly.

Note that the first two cases might occur only for short periods of time $\tau$ since $K_{\hat{N}_{a}}$ is oscillating. Furthermore, we note that the squeezing parameter $\tilde{\mathcal{D}}_{2}(\tau)$ affects the peak value of the non-Gaussianity because it enters into $|K_{\hat{N}_{a}}|$ through the $F$-coefficients (3.4). The dependence is non-trivial, but we will consider the analytic case for constant squeezing below. However, in general, when $|\mu_{\rm{c}}|\gg|K_{\hat{N}_{a}}|$, we see from (24) that the non-Gaussianity is independent of $\tilde{\mathcal{D}}_{2}(\tau)$ and can be accurately modelled by the standard optomechanical Hamiltonian without mechanical squeezing.

Let us now consider two specific cases where the squeezing term is either constant or modulated.

5.2 Applications: Constant squeezing parameter

Here we assume that the rescaled squeezing parameter is constant, with $\tilde{\mathcal{D}}_{2}(\tau)=\tilde{d}_{2}$. This case is equivalent to the case where the mechanical oscillation frequency $\omega_{\mathrm{m}}$ is shifted by a constant amount and where the initial state is a squeezed coherent state, see D. We begin by deriving analytic expressions for the coefficients in (3.4) given this choice of parameters.

5.2.1 Decoupled dynamics

We use the methods discussed in Section 3 to start by solving the differential equations (3.5). We find the solutions $P_{11}=P_{22}=\cos{\zeta\tau}$, where we define $\zeta:=\sqrt{1+4\,\tilde{d}_{2}}$. This, in turn, yields the following Bogoliubov coefficients (defined in (3.5)):

[TABLE]

Furthermore, we find $\xi(\tau)=\cos{\zeta\tau}-\frac{i}{\zeta}\sin{\zeta\tau}$, which in turn can be integrated to obtain the coefficients (3.4), which now read

[TABLE]

where $\rm{sinc}(x)=\sin(x)/x$. Since $\tilde{\mathcal{D}}_{1}=0$, all other coefficients are zero. The functions (5.2.1) now fully determine the time evolution through (3.4).

5.2.2 Quadratures

To gain intuition about the evolution of the system, we include plots of the optical quadratures. These can be found in Figure 2. The quadratures are the expectation values of $\hat{x}_{1}=(\hat{a}^{\dagger}+\hat{a})/\sqrt{2}$ and $\hat{p}_{1}=i(\hat{a}^{\dagger}-\hat{a})/\sqrt{2}$ and would correspond to classical trajectories in phase space. The full expression for the expectation values $\braket{\hat{x}_{1}}$ and $\braket{\hat{p}_{1}}$ can be found in (C) in C.

In Figures 2a, 2b, 2c and 2d, we have plotted the quadratures for $\mu_{\mathrm{c}}=1$, $\mu_{\mathrm{m}}=0$, $\tilde{g}_{0}=1$ and increasing values of $\tilde{d}_{2}$. While it is generally difficult to engineer a coupling of this magnitude, these values are chosen as example values to demonstrate the scaling behaviour of $\delta(\tau)$. Similarly in Figures 2e, 2f, 2g and 2h, we have plotted the quadratures for $\mu_{\mathrm{c}}=1$, $\mu_{\mathrm{m}}=1$, $\tilde{g}_{0}=1$ and again increasing values of $\tilde{d}_{2}$. To show the directionality of the evolution, the colour of the curve starts as light blue for $\tau=0$ and becomes increasingly darker as $\tau$ increases. We observe that the addition of mechanical squeezing causes the system to trace out highly complex trajectories, compared with the case when $\tilde{d}_{2}=0$.

5.2.3 Measure of non-Gaussianity

We now proceed to compute the non-Gaussianity $\delta(\tau)$, defined in (18), of the state evolving at constant squeezing parameter. A fully analytic expression for $\delta(\tau)$ exists but is too cumbersome to include here. Instead, we plot the measure of non-Gaussianity in Figure 3. In the first row of Figure 3, we present a comparison between the full measure $\delta$ (Figure 3d) and the lower and upper bounds $\delta_{\rm{min}}$ and $\delta_{\rm{max}}$ provided by the Araki–Lieb inequality in Figures 3e and 3f.

We note that the non-Gaussianity increases for large light–matter coupling $\tilde{g}_{0}$ and large coherent state parameter $\mu_{\mathrm{c}}$. This feature was also observed for standard optomechanical systems in [9]. However, the most striking feature here is that the larger $\tilde{d}_{2}$ is, the less non-Gaussian the system becomes. To understand why this is the case, we examine the dependence on $\tilde{d}_{2}$ in the function $|K_{\hat{N}_{a}}|$, since this determines the behaviour of the non-Gaussianity in certain regime, as discussed in Section 5.1. Using the expression (5.2.1) we find

[TABLE]

For large $\tilde{d}_{2}$, and therefore large $\zeta$, the first term inside the brackets dominates and for $\zeta\tau\neq n\pi$ with integer $n$, we are left with $|K_{\hat{N}_{a}}|\sim\tilde{g}_{0}\,\sin^{2}(\zeta\,\tau)/\zeta^{2}$. In general, we find $\lim_{\tilde{d}_{2}\rightarrow\infty}|K_{\hat{N}_{a}}|^{2}=0$. The consequences for the non-Gaussianity are difficult to predict given the complexity of the expressions, but we note that the mechanical symplectic eigenvalue $\nu_{\rm{Me}}$ decreases, while the optical symplectic eigenvalue $\nu_{\rm{Op}}$ increases.

Furthermore, the quantity $\theta(\tau)=2\,\left(F_{\hat{N}_{a}^{2}}+F_{\hat{N}_{a}\,\hat{B}_{+}}\,F_{\hat{N}_{a}\,\hat{B}_{-}}\right)$ is given by

[TABLE]

We find that $\lim_{\tilde{d}_{2}\rightarrow\infty}\theta(\tau)=0$. We then look at the symplectic eigenvalues (5.1) in this limit. We find that $\nu_{\rm{Me}}\rightarrow 1$, and $\nu_{\rm{Op}}\rightarrow 1$, which means that both the upper and the lower bounds of the non-Gaussianity tend to zero, and hence $\delta(\tau)\rightarrow 0$ as $\tilde{d}_{2}$ increases. We conclude that increasing the amount of constant squeezing in the system reduces the overall non-Gaussianity of the state.

5.3 Applications: Modulated squeezing parameter

In this section, we consider a modulated squeezing term. The dimensionless squeezing $\tilde{\mathcal{D}}_{2}(\tau)=\mathcal{D}_{2}(t)/\omega_{\rm{m}}$ is time-dependent and of the form

[TABLE]

where $\tilde{d}_{2}=d_{2}/\omega_{\mathrm{m}}$ is the amplitude of the squeezing and $\Omega_{0}$ denotes the time-scale of squeezing.666We remind the reader that our rescaled quantities require us to use $\tilde{d}_{2}=d_{2}/\omega_{\mathrm{m}}$ and we define $\Omega_{0}=\omega_{0}/\omega_{\mathrm{m}}$.

The differential equations in (3.5) are not generally analytically solvable for arbitrary choices of $\tilde{\mathcal{D}}_{2}(\tau)$. However, for the choice of $\tilde{\mathcal{D}}_{2}(\tau)$ in (30), both equations have a known form. Consider the differential equation for $P_{11}$, which we reprint here for convenience,

[TABLE]

Equation (31) is that of a parametric oscillator, which is used elsewhere in physics to describe, for example, a driven pendulum. As shown in B, the equation for the integral of $P_{22}$ (61) takes the same form.

The equation (31) is known as the Mathieu equation. In its most general form, and using conventional notation, it reads:

[TABLE]

where $a,q,$ and $x$ are real parameters. The general solutions to this equation are linear combinations of functions known as the Mathieu cosine $C(a,q,x)$ and Mathieu sine $S(a,q,x)$, the exact form of which will be determined by the boundary conditions for $y$.

To find which values the $a$, $q$ and $x$ parameters correspond to, we note that the cosine-term in $\tilde{\mathcal{D}}_{2}(\tau)$ has the argument $\Omega_{0}\,\tau$, which means that we must rescale time $\tau$ as $\tau^{\prime}=\Omega_{0}\tau/2$. Inserting our expression for $\tilde{\mathcal{D}}_{2}(\tau)$ and using the chain-rule to change variables from $\tau$ to $\tau^{\prime}$, we rewrite the equation for $P_{11}$ as

[TABLE]

where we identify the variables $a=4/\Omega_{0}^{2}$, and $q=-8\,\tilde{d}_{2}/\Omega_{0}^{2}$. The boundary conditions $P_{11}(0)=1$ and $\dot{P}_{11}(0)=0$ will yield the Mathieu cosine $C(a,q,x)$, and for $I_{P_{22}}$ as the solution, and the boundary conditions $I_{P_{22}}(0)=0$ and $\dot{I}_{P_{22}}(0)=1$ will yield the Mathieu sine $S(a,q,x)$ as the solution. For our specific choice of $\mathcal{D}_{2}(\tau)$ in (30), the system is resonant at $\Omega_{0}=2$ (see E), which means that $a=1$ and $q=-2\tilde{d}_{2}$.

5.4 Approximate solutions at resonance

The Mathieu equations are notoriously difficult to evaluate numerically. Instead, we use a two-scale method to derive perturbative solutions to $P_{11}$ and $I_{P_{22}}$. The perturbative solutions are valid for $\tilde{d}_{2}\ll 1$ and make use of specific resonance conditions to ensure that the solutions do not diverge. See E for the full derivation, where we also show that these approximate solutions correspond exactly to the more physically intuitive rotating-wave approximation (RWA) when $\tau\gg 1$. For smaller values of $\tau$, the approximate solutions are still valid, but they cannot be interpreted as equivalent to the RWA.

The squeezing term is resonant when $\Omega_{0}=2$. We find that the approximate solutions for $P_{11}$ and $I_{P_{22}}$ (the integral of $P_{22}$) are given by, respectively,

[TABLE]

We then compute $\xi(\tau)$ in (131). We assume that $\tilde{d}_{2}\tau\ll 1$ to find

[TABLE]

where we have expanded the hyperbolic functions to second order. By using the relations between $\xi(\tau)$ and the Bogoliubov coefficients (72), we find that the Bogoliubov condition is approximately satisfied as:

[TABLE]

With this expression, we can now compute the non-zero $F$-coefficients (3.4):

[TABLE]

where we have discarded terms with $\tilde{d}_{2}^{3}$. With these expressions, we are ready to compute the non-Gaussianity of the system when the squeezing is applied at mechanical resonance.

5.4.1 Measure of non-Gaussianity at resonance

We first compute the full measure of non-Gaussianity $\delta(\tau)$ and plot the results in Figure 4. Figure 4a shows the full measure $\delta(\tau)$, the lower bound $\delta_{\rm{min}}(\tau)$ and the upper bound $\delta_{\rm{max}}(\tau)$ as a function of $\tau$ for the parameter $\tilde{g}_{0}=1$, $\mu_{\rm{c}}=1$, $\tilde{d}_{2}=0.1$, and $\mu_{\rm{m}}=0$ as a function of time $\tau$ and the squeezing $\tilde{d}_{2}$. The second plot in Figure 4b also shows the full measure $\delta(\tau)$, the lower bound $\delta_{\rm{min}}(\tau)$ and the upper bound $\delta_{\rm{max}}(\tau)$ as a function of $\tilde{g}_{0}$ at $\tau=\pi$, $\tilde{d}_{2}=0.1$, $\mu_{\rm{c}}=1$, and $\mu_{\rm{m}}=0$. We find that the non-Gaussianity increases with $\tilde{g}_{0}$, as expected.

In Figure 4, we considered $\tilde{d}_{2}=0.1$; a value consistent with the validity of the approximate solutions to the Mathieu equation. For this value, the non-Gaussianity is found to increase very slightly with $\tilde{d}_{2}$. To demonstrate this, we consider the regime where $1\ll 2|\mu_{\rm{c}}|\ll|K_{\hat{N}_{a}}|$, which occurs when $2|\mu_{\rm{c}}|\ll\tilde{g}_{0}$ for specific values of $\tau$. In this regime, the non-Gaussianity was approximately given by $s_{V}\left(2\,|K_{\hat{N}_{a}}||\mu_{\rm{c}}|\right)$ (25). Given the functions (5.4), we find that

[TABLE]

where we have again removed terms proportional to $\tilde{d}_{2}^{3}$ and $\tilde{d}_{2}^{4}$. The behaviour of $|K_{\hat{N}_{a}}|^{2}$ is markedly different compared with the constant case. Firstly, while $|K_{\hat{N}_{a}}|^{2}$ still oscillates, it also increases with $\tau$ and $\tilde{d}_{2}$. If we consider the leading term with $\tau^{2}$, we find that the non-Gaussianity scales with $\delta\sim\ln(\tau\,\tilde{d}_{2}\,\tilde{g}_{0})$, which confirms that in this specific regime, the non-Gaussianity increases logarithmically with $\tau$, $\tilde{d}_{2}$, and $\tilde{g}_{0}$. We conclude that squeezing is not necessarily detrimental to the non-Gaussianity if the squeezing is modulated at resonance, although more work needs to be done to ascertain the full interplay between the two effects.

6 Discussion

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

6.1 Advantages over direct numerical simulations

With our techniques, we have shown that it is possible to analytically solve the dynamics of a nonlinear optomechanical system even when the mechanical squeezing is time-dependent. To emphasise this point, we wish to compare our approach, which relies on numerically solving the differential equations in (3.5), with a general numerical method using a standard higher-order Runge–Kutta solver to evolve a state in a truncated Hilbert space, e.g. using the Python library QuTiP [51].

When the dynamics is solved with a Runge–Kutta method, the continuous variable (pure) states are represented as finite-dimensional vectors in a truncated Hilbert space. When the system is nonlinear, information about the state is quickly distributed across large sectors of the Hilbert space. If the computational Hilbert space is too small, numerical inaccuracies quickly enter into the evolution, as a result of truncating the vectors. It follows that the dimension of the Hilbert space must be large enough to prevent this, which requires significant amounts of computer memory. It is also very difficult to consider parameters of the magnitude $\tilde{g}_{0}=10$ and $\tilde{d}_{2}=10$, as done in this work, since these cause the system to evolve very rapidly and, consequently, require the evolution of the system to be calculated using smaller and smaller time intervals.

The methods developed here excel at treating systems numerically for large parameters $\tilde{g}_{0},\tilde{d}_{2},\mu_{\mathrm{c}}$ and $\mu_{\mathrm{m}}$. However, we note that it becomes increasingly difficult to numerically evaluate the dynamics at longer times $\tau$ when the system of differential equations (3.5) is numerically solved for arbitrary functions $\tilde{\mathcal{D}}_{2}(\tau)$. The difficulty is primarily caused by the double integral that determines the coefficient $F_{\hat{N}_{a}^{2}}$ in (3.4), which must be evaluated numerically. For each value of $\tau$, the integral will be evaluated from [math] to the final $\tau^{\prime}$, and then from [math] to $\tau$. As a result, the integrals take an increasingly long time to evaluate for large $\tau$. We therefore conclude that the key strength in our method lies in evaluating the state of the system at early times $\tau\in(0,2\pi)$ for large parameters $\mu_{\mathrm{c}},\tilde{g}_{0}$, and $\tilde{d}_{2}$. We also emphasise that, the computation using our methods is numerically exact, which a naive computation using QuTiP or a similar library is not.

To conclude, our methods allow for the evaluation of the state of the system with large parameters, e.g. $\tilde{g}_{0}=100$ and $\tilde{d}_{2}=10$, which would be nearly impossible to perform with QuTiP or a similar library unless one had access to significantly more computational resources.

6.2 Competing behaviours of nonlinearity and squeezing

We concluded from Figure 3 that the addition of a constant squeezing term has a detrimental effect on the non-Gaussianity of the system. We also noted that inclusion of a constant squeezing term is equivalent to changing the mechanical trapping frequency $\omega_{\rm{m}}$ to a specific value and starting the computation with a squeezed coherent state (see D). With this interpretation, our results also show that an initially squeezed state evolving under the optomechanical Hamiltonian can be expected to exhibit less non-Gaussianity compared with coherent states. The reason for this overall behaviour can be found by simple inspection of the total Hamiltonian. If a strong squeezing term is included in the Hamiltonian (1), it dominates over the interaction term, leading to a decrease in the non-Gaussianity. However, such a process is not fully monotonic, since an increase of the squeezing parameter does not always decrease the non-Gaussianity. This is, however, reasonable, as it cannot be expected that only the relative weight of the two parts of the Hamiltonian matter; the precise dynamics is much more complex, and the non-Gaussianity depends on the entire state, which is driven by the full Hamiltonian.

The finding that the non-Gaussianity increases with both time $\tau$ and $\tilde{d}_{2}$ when modulated at mechanical resonance is interesting and warrants further investigation. We leave this to future work.

7 Conclusion

In this work, we solved the time-evolution of a nonlinear optomechanical system with a time-dependent mechanical displacement term and a time-dependent mechanical single-mode squeezing term. We found analytic expressions for all first and second moments of the quadratures of the nonlinear system and used them to compute the amount of non-Gaussianity of the state. We considered both constant and modulated squeezing parameter, and found that a squeezing parameter modulated at resonance results in the Mathieu equations, for which we provide approximate solutions equivalent to the rotating-wave approximation.

In general, we find that the relationship between the squeezing and non-Gaussianity is highly nontrivial. The inclusion of a mechanical squeezing term in the Hamiltonian, which is equivalent to starting with a coherent squeezed state evolving with the standard optomechanical Hamiltonian with a shifted mechanical frequency, decreases the overall non-Gaussianity of the state. If the squeezing term is modulated at mechanical resonance, however, we found that the non-Gaussianity increases with both time and the squeezing parameter in specific regimes. These results hold interesting implications for quantum control of nonlinear optomechanical systems.

Our results also suggest that the combination of non-Gaussian resources and mechanical squeezing may not necessarily be beneficial if the application relies specifically on the non-Gaussian character of the state. However, more work is needed to conclude if this has a significant effect on, for example, sensing applications. More work is also necessary to properly study the instabilities of the full solutions to the Mathieu equations and how they affect the dynamics. The effect of squeezing the optical rather than mechanical mode is another question we defer to future work.

The decoupling methods demonstrated here constitute an important step towards fully characterising nonlinear systems with mechanical squeezing and can be used both to model experimental systems and to test numerical methods. The results presented here apply to any system with the same characteristic cubic Hamiltonian interaction term and single mode squeezing term. More broadly, the Lie algebra method can be used to solve any system where a finite set of Lie algebra operators has been identified. Our work can further be extended to more complicated quadratic Hamiltonians of bosonic modes, such as Dicke-like models [52], which would allow for applications in other areas of physics to be developed.

Acknowledgements

We would like to thank Antonio Pontin, Peter F. Barker, Robert Delaney, Doug Plato, and Ivette Fuentes for useful 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). DEB acknowledges support from the CEITEC Nano RI.