Combining Floquet and Lyapunov techniques for time-dependent problems in optomechanics and electromechanics

In the past years, cavity optomechanics and electromechanics [1] became a mature field. Starting with first experiments demonstrating ground-state cooling of mechanical motion [2, 3], electromagnetic fields are now commonly used to manipulate and measure quantum states of mechanical resonators. Examples of nonclassical mechanical states include squeezed states, prepared using reservoir engineering [4–7] or parametric effects [8], as well as Fock states and their superpositions [9–13]. Preparation of these states presents a crucial first step towards creating more complex mechanical states such as truly macroscopic mechanical superpositions [14–16] or towards generation of entanglement in optomechanical and electromechanical systems [17–20].

In many applications, the system of interest is described by a Hamiltonian that is quadratic in the canonical operators of radiation and mechanical motion. As the unitary dynamics are accompanied by single-excitation loss and gain and—in case of backaction-evading measurements—homodyne detection (nonunitary processes which are linear in the canonical operators), the phase-space description of the quantum state retains a Gaussian character [21]. We can therefore characterise the state by the first and second statistical moments of the quasiprobability density—the mean vector and covariance matrix of the canonical operators—instead of the density operator. This presents a great advantage since we can formulate equations of motion for the mean and covariance matrix [22, 23] which can be solved more efficiently than the master equation that describes the evolution of the density operator. Particularly the covariance matrix is important for investigating Gaussian quantum features of the resulting quantum state such as quadrature squeezing or entanglement [21]; in the steady state, its form can be obtained by solving an algebraic equation of the Lyapunov (for unconditional dynamics) or Riccati (for conditional evolution) type [22, 23].

A crucial obstacle for a more widespread application of these techniques is the explicit time dependence of the driving electromagnetic fields. Dissipative preparation of mechanical states [4, 6, 7, 14, 16–18, 20, 24–28] and tomographic backaction-evading measurements of mechanical motion [5, 29–36] rely on driving the system with multiple fields at different frequencies. Similarly, parametric squeezing requires modulation of the optical spring [8, 37–39]. Finally, optomechanical systems with multiple mechanical modes [40–47] often involve driving at the sideband of each mechanical mode; this effectively produces driving at the sum and difference of the mechanical frequencies.

In all of these approaches, the resulting optomechanical Hamiltonian is time-dependent. The steady-state Lyapunov equation can then only be applied under the rotating wave approximation (RWA) which neglects fast oscillating terms in the interaction and only keeps those that are resonant. This approximation can work well as long as both the cavity damping rate and optomechanical coupling strength are much smaller than the mechanical frequency. When one of these quantities becomes comparable to the mechanical frequency—or, in the multimode case, the difference of two mechanical frequencies—the RWA breaks down and the fast oscillating terms can no longer be neglected. In this case, Floquet techniques can be used to find the steady state by expanding the system operators into their Fourier components [48–52]. Their application is, however, often cumbersome and lacks a simple, general approach.

In this article, we show how the Lyapunov and Floquet techniques can be combined to obtain a straightforward approach to finding the steady state of optomechanical and electromechanical systems with periodic Hamiltonians. We derive a general method which can be directly applied to the Lyapunov equation of a time-dependent quantum system to obtain a time-independent Lyapunov equation in the Floquet space. In this process, the Floquet formalism is used to introduce frequency components associated with the periodicity of the time-dependent Hamiltonian. Importantly, the procedure does not require the derivation of the Floquet master equation which is the common approach; instead, the frequency components of the relevant quantities can be directly read off from the time-dependent equations of motion for the canonical operators. This feature makes our strategy particularly easy to implement for a wide variety of critical problems in optomechanics and electromechanics.

After presenting the general method, we apply it to the problem of dissipative mechanical squeezing. We consider both the usual squeezing obtained by cooling a mechanical Bogoliubov mode [4] as well as the recently proposed combination of dissipative and parametric squeezing possible for levitated particles [38] and analyse the resulting squeezing in both cases beyond the usual RWA. This approximation is invalidated not only by the nonzero sideband ratio (the ratio between the cavity decay rate and the mechanical frequency) but also by the ultrastrong coupling regime which is being reached in an evergrowing number of optomechanical and electromechanical systems [53–62]. A detailed systematic analysis of the effects that the counterrotating terms will thus have important consequences not only for immediate applications in state-of-the-art optomechanical and electromechanical systems but also for advanced state preparation techniques, allowing a more detailed analysis of protocols that combine Gaussian and non-Gaussian elements [14–16].

2.1. Lyapunov equation

The dynamics of an open quantum system can be described by the master equation (ℏ = 1)

$$\begin{equation}\dot {\hat{\rho }}=-\text{i}\left[\hat{H},\hat{\rho }\right]+\sum _{n}\mathcal{D}\left[{\hat{{\jmath}}}_{n}\right]\hat{\rho },\end{equation} \tag{ 1 }$$

where $\mathcal{D}\left[\hat{{\jmath}}\right]\hat{\rho }=\hat{{\jmath}}\hat{\rho }{\hat{{\jmath}}}^{{\dagger}}-\frac{1}{2}{\hat{{\jmath}}}^{{\dagger}}\hat{{\jmath}}\hat{\rho }-\frac{1}{2}\hat{\rho }{\hat{{\jmath}}}^{{\dagger}}\hat{{\jmath}}$ is the Lindblad superoperator with collapse operator $\hat{{\jmath}}$ [63]. The full dynamics can, in principle, be solved from this equation but this is not always feasible in practice. Optomechanical systems are one such example; the large thermal occupation of the mechanical bath requires prohibitively large Hilbert-space dimensions in the Fock-state basis to be considered for numerical simulations.

Full description of the dynamics without solving the master equation (1) is possible for Gaussian systems, which are characterised by Gaussian quasiprobability distributions in phase space and as such can be fully described by the first and second statistical moments of these distributions. A system is Gaussian when its Hamiltonian Ĥ is quadratic and the collapse operators ${\hat{{\jmath}}}_{n}$ linear in canonical operators; then, we can solve for the covariance matrix describing the state of the system instead of the density matrix. (Note that the first moments remain zero when the Hamiltonian contains no linear terms and the collapse operators no constant terms.)

The covariance matrix Γ obeys the Lyapunov equation

$$\begin{equation}\dot {\mathbf{\Gamma }}=\mathbf{A}\mathbf{\Gamma }+\mathbf{\Gamma }{\mathbf{A}}^{\top }+\mathbf{N},\end{equation} \tag{ 2 }$$

where

$$\begin{equation}{{\Gamma}}_{ij}=\langle {\hat{r}}_{i}{\hat{r}}_{j}+{\hat{r}}_{j}{\hat{r}}_{i}\rangle -2\langle {\hat{r}}_{i}\rangle \langle {\hat{r}}_{j}\rangle ,\end{equation} \tag{ 3a }$$

$$\begin{equation}\hat{\mathbf{r}}={\left({\hat{q}}_{1},{\hat{p}}_{1},{\hat{q}}_{2},{\hat{p}}_{2},\dots ,\;{\hat{q}}_{N},{\hat{p}}_{N}\right)}^{\top },\end{equation} \tag{ 3b }$$

with the canonical commutation relation

$$\begin{equation}\left[\hat{\mathbf{r}},{\hat{\mathbf{r}}}^{\top }\right]=\text{i}{\sigma }_{N}=\text{i}\left(\underset{N}{\underbrace{\sigma \oplus \sigma \oplus \cdots \oplus \sigma }}\right),\quad \sigma =\left(\begin{pmatrix}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{pmatrix}\right),\end{equation} \tag{ 4 }$$

where σ is the one-mode symplectic matrix. The Lyapunov equation—including precise forms for the drift (A) and diffusion (N) matrices can be found from the master equation [23]. Starting from Langevin equations

$$\begin{equation}\dot {\hat{\mathbf{r}}}=\mathbf{D}\hat{\mathbf{r}}+\hat{\xi },\end{equation} \tag{ 5 }$$

one finds that the drift matrix is identical to the dynamical matrix, A = D; the diffusion matrix can be expressed in terms of the vector of noise operators $\hat{\xi }={\left({\hat{q}}_{1,\mathrm{i}\mathrm{n}},{\hat{p}}_{1,\mathrm{i}\mathrm{n}},\dots ,\;{\hat{q}}_{N,\mathrm{i}\mathrm{n}},{\hat{p}}_{N,\mathrm{i}\mathrm{n}}\right)}^{\top }$ as

$$\begin{equation}\mathbf{N}=\left\langle \hat{\xi }\left(t\right)\hat{\xi }{\left(t\right)}^{\top }\right\rangle .\end{equation} \tag{ 6 }$$

The Lyapunov equation (2) is a powerful tool in optomechanics where it can be used to evaluate the dynamics of the system or to find its steady state from the algebraic form of the Lyapunov equation

$$\begin{equation}\mathbf{A}\mathbf{\Gamma }+\mathbf{\Gamma }{\mathbf{A}}^{\top }+\mathbf{N}=0.\end{equation} \tag{ 7 }$$

The process of solving equation (2) or (7) is straightforward as long as the system operators are not explicitly time dependent. When the Hamiltonian or coupling to the environment becomes time dependent, the drift matrix (and possibly also the diffusion matrix) becomes time dependent as well. While the differential Lyapunov equation (2) can still be solved in this case, a straightforward solution of the algebraic Lyapunov equation (7) is not possible.

2.2. The Floquet–Lyapunov method

A common approach to problems with periodic time-dependence is based on the Floquet formalism, which allows us to transform a set of periodic linear differential equations into a larger set of linear differential equations without explicit time dependence [64]. Here, we will show how to use this method to obtain a time-independent Lyapunov equation for explicitly time-dependent dynamics. For simplicity, we will assume that only the Hamiltonian is τ-periodic, Ĥ(t + τ) = Ĥ(t), leading to a periodic drift matrix A(t + τ) = A(t), while the coupling to the environment (and the diffusion matrix N) remains explicitly time independent. We will only present a general recipe for obtaining a time-independent Lyapunov equation in this section; the mathematical justification of these steps is presented in appendix A.

We start by expressing the operators and the drift matrix in terms of their Fourier components at frequency ω = 2π/τ,

$$\begin{equation}\hat{\mathbf{r}}={\hat{\mathbf{r}}}^{\left(0\right)}+\sqrt{2}\sum _{n=1}^{\infty }\left[{\hat{\mathbf{r}}}_{\text{c}}^{\left(n\right)}\;\mathrm{cos}\left(n\omega t\right)+{\hat{\mathbf{r}}}_{\text{s}}^{\left(n\right)}\;\mathrm{sin}\left(n\omega t\right)\right],\end{equation} \tag{ 8a }$$

$$\begin{equation}\mathbf{A}\left(t\right)={\mathbf{A}}^{\left(0\right)}+\sqrt{2}\sum _{n=1}^{\infty }\left[{\mathbf{A}}_{\text{c}}^{\left(n\right)}\;\mathrm{cos}\left(n\omega t\right)+{\mathbf{A}}_{\text{s}}^{\left(n\right)}\;\mathrm{sin}\left(n\omega t\right)\right],\end{equation} \tag{ 8b }$$

$$\begin{equation}\hat{\xi }={\hat{\xi }}^{\left(0\right)}+\sqrt{2}\sum _{n=1}^{\infty }\left[{\hat{\xi }}_{\text{c}}^{\left(n\right)}\;\mathrm{cos}\left(n\omega t\right)+{\hat{\xi }}_{\text{s}}^{\left(n\right)}\;\mathrm{sin}\left(n\omega t\right)\right],\end{equation} \tag{ 8c }$$

The frequency components ${\hat{\mathbf{r}}}^{\left(0\right)},{\hat{\mathbf{r}}}_{\text{c,}\;\text{s}}^{\left(n\right)}$ fulfil the canonical commutation relation, $\left[{\hat{\mathbf{r}}}^{\left(0\right)},{\hat{\mathbf{r}}}^{\left(0\right)\top }\right]=\left[{\hat{\mathbf{r}}}_{\text{k}}^{\left(n\right)},{\hat{\mathbf{r}}}_{\text{k}}^{\left(n\right)\top }\right]=\text{i}{\sigma }_{N}$ (here, k ∈ {c,s} and operators from different frequency components commute), and obey the Langevin equations

$$\begin{equation}{\dot {\hat{\mathbf{r}}}}^{\left(0\right)}={\mathbf{A}}^{\left(0\right)}{\hat{\mathbf{r}}}^{\left(0\right)}+\sum _{m=1}^{\infty }\left({\mathbf{A}}_{\text{c}}^{\left(m\right)}{\hat{\mathbf{r}}}_{\text{c}}^{\left(m\right)}+{\mathbf{A}}_{\text{s}}^{\left(m\right)}{\hat{\mathbf{r}}}_{\text{s}}^{\left(m\right)}\right)+{\hat{\xi }}^{\left(0\right)},\end{equation} \tag{ 9a }$$

$$\begin{align}\hfill {\dot {\hat{\mathbf{r}}}}_{\text{c}}^{\left(n\right)}& ={\mathbf{A}}^{\left(0\right)}{\hat{\mathbf{r}}}_{\text{c}}^{\left(n\right)}+{\mathbf{A}}_{\text{c}}^{\left(n\right)}{\hat{\mathbf{r}}}^{\left(0\right)}-n\omega {\mathbf{r}}_{\text{s}}^{\left(n\right)}+\frac{1}{\sqrt{2}}\sum _{m=1}^{\infty }{\mathbf{A}}_{\text{c}}^{\left(m\right)}\left({\hat{\mathbf{r}}}_{\text{c}}^{\left(m-n\right)}+{\hat{\mathbf{r}}}_{\text{c}}^{\left(m+n\right)}+{\hat{\mathbf{r}}}_{\text{c}}^{\left(n-m\right)}\right)\hfill \\ \hfill & \quad +\;\frac{1}{\sqrt{2}}\sum _{m=1}^{\infty }{\mathbf{A}}_{\text{s}}^{\left(m\right)}\left({\hat{\mathbf{r}}}_{\text{s}}^{\left(m-n\right)}+{\hat{\mathbf{r}}}_{\text{s}}^{\left(m+n\right)}-{\hat{\mathbf{r}}}_{\text{s}}^{\left(n-m\right)}\right)+{\hat{\xi }}_{\text{c}}^{\left(n\right)},\hfill \end{align} \tag{ 9b }$$

$$\begin{align}\hfill {\dot {\hat{\mathbf{r}}}}_{\text{s}}^{\left(n\right)}& ={\mathbf{A}}^{\left(0\right)}{\hat{\mathbf{r}}}_{\text{s}}^{\left(n\right)}+{\mathbf{A}}_{\text{s}}^{\left(n\right)}{\hat{\mathbf{r}}}^{\left(0\right)}+n\omega {\mathbf{r}}_{\text{c}}^{\left(n\right)}+\frac{1}{\sqrt{2}}\sum _{m=1}^{\infty }{\mathbf{A}}_{\text{c}}^{\left(m\right)}\left(-{\hat{\mathbf{r}}}_{\text{s}}^{\left(m-n\right)}+{\hat{\mathbf{r}}}_{\text{s}}^{\left(m+n\right)}+{\hat{\mathbf{r}}}_{\text{s}}^{\left(n-m\right)}\right)\hfill \\ \hfill & \quad +\;\frac{1}{\sqrt{2}}\sum _{m=1}^{\infty }{\mathbf{A}}_{\text{s}}^{\left(m\right)}\left({\hat{\mathbf{r}}}_{\text{c}}^{\left(m-n\right)}-{\hat{\mathbf{r}}}_{\text{c}}^{\left(m+n\right)}+{\hat{\mathbf{r}}}_{\text{c}}^{\left(n-m\right)}\right)+{\hat{\xi }}_{\text{s}}^{\left(n\right)}.\hfill \end{align} \tag{ 9c }$$

The individual frequency components of the noise operators have the same correlation properties as the initial noise, $\left\langle {\hat{\xi }}^{\left(0\right)}{\hat{\xi }}^{\left(0\right)\top }\right\rangle =\left\langle {\hat{\xi }}_{\text{k}}^{\left(n\right)}{\hat{\xi }}_{\text{k}}^{\left(n\right)\top }\right\rangle =\mathbf{N}$, and are mutually uncorrelated, $\left\langle {\hat{\xi }}^{\left(0\right)}{\hat{\xi }}_{\text{k}}^{\left(n\right)\top }\right\rangle =\left\langle {\hat{\xi }}_{\text{k}}^{\left(n\right)}{\hat{\xi }}_{\text{l}}^{\left(m\right)\top }\right\rangle =0$ for n ≠ m or k ≠ l. Note also that in equation (9), only the frequency components with ±m ± n > 0 should be included when evaluating the sums; effectively, we have ${\hat{\mathbf{r}}}_{\text{k}}^{\left(n\right)}=0$ for n ⩽ 0.

Collecting the frequency components in the vector ${\hat{\mathbf{r}}}_{\mathrm{F}}={\left({\hat{\mathbf{r}}}^{\left(0\right)},{\hat{\mathbf{r}}}_{\text{c}}^{\left(1\right)},{\hat{\mathbf{r}}}_{\text{s}}^{\left(1\right)},{\hat{\mathbf{r}}}_{\text{c}}^{\left(2\right)},{\hat{\mathbf{r}}}_{\text{s}}^{\left(2\right)},\dots \right)}^{\top }$ (the components ${\hat{\mathbf{r}}}^{\left(0\right)},{\hat{\mathbf{r}}}_{\text{k}}^{\left(n\right)}$, which are column vectors, are arranged into a longer column vector), we can formulate the Langevin equation

$$\begin{equation}{\dot {\hat{\mathbf{r}}}}_{\mathrm{F}}={\mathbf{A}}_{\mathrm{F}}{\hat{\mathbf{r}}}_{\mathrm{F}}+{\hat{\xi }}_{\mathrm{F}},\end{equation} \tag{ 10 }$$

where we introduced the drift matrix

$$\begin{equation}{\mathbf{A}}_{\mathrm{F}}=\left(\begin{pmatrix}{cccccc}\hfill {\mathbf{A}}^{\left(0\right)}\hfill & \hfill {\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill {\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \hfill \\ \hfill {\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill -\omega {\mathbf{I}}_{N}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}+{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill \hfill \\ \hfill {\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill \omega {\mathbf{I}}_{N}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}-\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left(-{\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}-{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill \dots \hfill \\ \hfill {\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}+{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left(-{\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill {\mathbf{A}}^{\left(0\right)}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(4\right)}\hfill & \hfill -2\omega {\mathbf{I}}_{N}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(4\right)}\hfill & \hfill \hfill \\ \hfill {\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}-{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill 2\omega {\mathbf{I}}_{N}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(4\right)}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}-\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(4\right)}\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill {\vdots}\hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{pmatrix}\right),\end{equation} \tag{ 11 }$$

with I_N denoting the N-mode (i.e., 2N-dimensional) identity matrix and ${\hat{\xi }}_{\mathrm{F}}={\left({\hat{\xi }}^{\left(0\right)},{\hat{\xi }}_{\text{c}}^{\left(1\right)},{\hat{\xi }}_{\text{s}}^{\left(1\right)},\dots \right)}^{\top }$ being the vector of noise operators in terms of their frequency components. In these expressions, we use the subscript F to signify that the quantities are defined in the Floquet space. We can now formulate the Lyapunov equation in the Floquet space

$$\begin{equation}{\dot {\mathbf{\Gamma }}}_{\mathrm{F}}={\mathbf{A}}_{\mathrm{F}}{\mathbf{\Gamma }}_{\mathrm{F}}+{\mathbf{\Gamma }}_{\mathrm{F}}{\mathbf{A}}_{\mathrm{F}}^{\top }+{\mathbf{N}}_{\mathrm{F}}\end{equation} \tag{ 12 }$$

with the time-independent drift matrix given by equation (11) and diffusion matrix N_F = diag(N, N, ...). The covariance matrix can be expressed in the block form

$$\begin{equation}{\mathbf{\Gamma }}_{\mathrm{F}}=\left(\begin{pmatrix}{ccc}\hfill {\mathbf{\Gamma }}^{\left(0\right)}\hfill & \hfill {\mathbf{\Gamma }}_{\text{c}}^{\left(01\right)}\hfill & \hfill {\mathbf{\Gamma }}_{\text{s}}^{\left(01\right)}\hfill \\ \hfill {\mathbf{\Gamma }}_{\text{c}}^{\left(10\right)}\hfill & \hfill {\mathbf{\Gamma }}_{\text{c}}^{\left(1\right)}\hfill & \hfill {\mathbf{\Gamma }}_{\text{cs}}^{\left(11\right)}\hfill & \hfill \dots \hfill \\ \hfill {\mathbf{\Gamma }}_{\text{s}}^{\left(10\right)}\hfill & \hfill {\mathbf{\Gamma }}_{\text{sc}}^{\left(11\right)}\hfill & \hfill {\mathbf{\Gamma }}_{\text{s}}^{\left(1\right)}\hfill \\ \hfill \hfill & \hfill {\vdots}\hfill & \hfill \hfill & \hfill \ddots \hfill \end{pmatrix}\right),\end{equation} \tag{ 13 }$$

where the blocks are defined as the covariances between the corresponding frequency components of the quadrature operators; we thus have

$$\begin{equation}\begin{aligned}\hfill {{\Gamma}}_{ij}^{\left(0\right)}=\left\langle {\hat{r}}_{i}^{\left(0\right)}{\hat{r}}_{j}^{\left(0\right)}+{\hat{r}}_{j}^{\left(0\right)}{\hat{r}}_{i}^{\left(0\right)}\right\rangle -2\left\langle {\hat{r}}_{i}^{\left(0\right)}\right\rangle \left\langle {\hat{r}}_{j}^{\left(0\right)}\right\rangle ,\quad \quad & {{\Gamma}}_{\text{k},ij}^{\left(0n\right)}=\left\langle {\hat{r}}_{i}^{\left(0\right)}{\hat{r}}_{\text{k},j}^{\left(n\right)}+{\hat{r}}_{\text{k},j}^{\left(n\right)}{\hat{r}}_{i}^{\left(0\right)}\right\rangle -2\left\langle {\hat{r}}_{i}^{\left(0\right)}\right\rangle \left\langle {\hat{r}}_{\text{k},j}^{\left(n\right)}\right\rangle ,\hfill \\ \hfill {{\Gamma}}_{\text{k},ij}^{\left(n\right)}=\left\langle {\hat{r}}_{\text{k},i}^{\left(n\right)}{\hat{r}}_{\text{k},j}^{\left(n\right)}+{\hat{r}}_{\text{k},j}^{\left(n\right)}{\hat{r}}_{\text{k},i}^{\left(n\right)}\right\rangle -2\left\langle {\hat{r}}_{\text{k},i}^{\left(n\right)}\right\rangle \left\langle {\hat{r}}_{\text{k},j}^{\left(n\right)}\right\rangle ,\quad \quad & {{\Gamma}}_{\text{kl},ij}^{\left(mn\right)}=\left\langle {\hat{r}}_{\text{k},i}^{\left(m\right)}{\hat{r}}_{\text{l},j}^{\left(n\right)}+{\hat{r}}_{\text{l},j}^{\left(n\right)}{\hat{r}}_{\text{k},i}^{\left(m\right)}\right\rangle -2\left\langle {\hat{r}}_{\text{k},i}^{\left(m\right)}\right\rangle \left\langle {\hat{r}}_{\text{l},j}^{\left(n\right)}\right\rangle .\hfill \end{aligned}\end{equation} \tag{ 14 }$$

The steady state can now be calculated from equation (12) by setting ${\dot {\mathbf{\Gamma }}}_{\mathrm{F}}=0$, which gives a linear, time-independent algebraic equation.

The system now becomes effectively infinite dimensional and must be suitably truncated for numerical simulations. This truncation has to ensure that the solution (contained in the zeroth Brillouin zone, i.e., the covariance block Γ⁽⁰⁾) converges. As we shall show in the following examples, only a small number of Brillouin zones is typically sufficient, making our approach more feasible than direct solution of the master equation or finding the long-time limit of the differential Lyapunov equation with time-dependent drift.

2.3. Illustration: optomechanical cooling

To illustrate our formalism, we now consider optomechanical sideband cooling and show how the Floquet–Lyapunov technique can be used to obtain the correct steady-state mechanical occupation in the rotating frame. For an optical cavity (with annihilation operator ĉ) coupled to a mechanical oscillator (annihilation operator $\hat{b}$), the linearised optomechanical Hamiltonian (in the lab frame of the mechanical motion) takes the form [1]

$$\begin{equation}\hat{H}={\Delta}{\hat{c}}^{{\dagger}}\hat{c}+{\omega }_{\text{m}}{\hat{b}}^{{\dagger}}\hat{b}+g\left({\hat{c}}^{{\dagger}}+\hat{c}\right)\left({\hat{b}}^{{\dagger}}+\hat{b}\right),\end{equation} \tag{ 15 }$$

where Δ = ω_c − ω_L is the detuning between the cavity frequency ω_c and driving frequency ω_L, ω_m is the mechanical frequency, and g is the coupling constant. This is a time-independent Hamiltonian and the dynamics can be solved directly by using the Lyapunov equation.

Introducing the canonical operators ${\hat{q}}_{1}=\left(\hat{c}+{\hat{c}}^{{\dagger}}\right)/\sqrt{2}$, ${\hat{p}}_{1}=-\text{i}\left(\hat{c}-{\hat{c}}^{{\dagger}}\right)/\sqrt{2}$ for the cavity field and ${\hat{q}}_{2}$, ${\hat{p}}_{2}$ for the mechanics (defined similarly), we can describe the dynamics using the Langevin equations

$$\begin{equation}{\dot {\hat{q}}}_{1}=-\kappa {\hat{q}}_{1}+{\Delta}{\hat{p}}_{1}+\sqrt{2\kappa }{\hat{q}}_{1,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 16a }$$

$$\begin{equation}{\dot {\hat{p}}}_{1}=-{\Delta}{\hat{q}}_{1}-\kappa {\hat{p}}_{1}-2g{\hat{q}}_{2}+\sqrt{2\kappa }{\hat{p}}_{1,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 16b }$$

$$\begin{equation}{\dot {\hat{q}}}_{2}=-\gamma {\hat{q}}_{2}+{\omega }_{\text{m}}{\hat{p}}_{2}+\sqrt{2\gamma }{\hat{q}}_{2,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 16c }$$

$$\begin{equation}{\dot {\hat{p}}}_{2}=-2g{\hat{q}}_{1}-{\omega }_{\text{m}}{\hat{q}}_{2}-\gamma {\hat{p}}_{2}+\sqrt{2\gamma }{\hat{p}}_{2,\mathrm{i}\mathrm{n}}.\end{equation} \tag{ 16d }$$

Here, κ and γ are the damping rates for the cavity field and the mechanical mode and the input noise operators of the optical and mechanical baths have the correlations

$$\begin{equation}\left\langle {\hat{q}}_{1,\mathrm{i}\mathrm{n}}\left(t\right){\hat{q}}_{1,\mathrm{i}\mathrm{n}}\left({t}^{\prime }\right)\right\rangle =\left\langle {\hat{p}}_{1,\mathrm{i}\mathrm{n}}\left(t\right){\hat{p}}_{1,\mathrm{i}\mathrm{n}}\left({t}^{\prime }\right)\right\rangle =\delta \left(t-{t}^{\prime }\right),\end{equation} \tag{ 17a }$$

$$\begin{equation}\left\langle {\hat{q}}_{2,\mathrm{i}\mathrm{n}}\left(t\right){\hat{q}}_{2,\mathrm{i}\mathrm{n}}\left({t}^{\prime }\right)\right\rangle =\left\langle {\hat{p}}_{2,\mathrm{i}\mathrm{n}}\left(t\right){\hat{p}}_{2,\mathrm{i}\mathrm{n}}\left({t}^{\prime }\right)\right\rangle =\left(2\overline{n}+1\right)\delta \left(t-{t}^{\prime }\right)\end{equation} \tag{ 17b }$$

with $\overline{n}$ denoting the average occupation of the thermal mechanical bath; all other correlations are zero. From these equations, we obtain the Lyapunov equation

$$\begin{equation}\dot {\mathbf{\Gamma }}=\mathbf{A}\mathbf{\Gamma }+\mathbf{\Gamma }{\mathbf{A}}^{\top }+\mathbf{N},\end{equation} \tag{ 18a }$$

$$\begin{equation}\mathbf{A}=\left(\begin{pmatrix}{cccc}\hfill -\kappa \hfill & \hfill {\Delta}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -{\Delta}\hfill & \hfill -\kappa \hfill & \hfill -2g\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\gamma \hfill & \hfill {\omega }_{\text{m}}\hfill \\ \hfill -2g\hfill & \hfill 0\hfill & \hfill -{\omega }_{\text{m}}\hfill & \hfill -\gamma \hfill \end{pmatrix}\right),\end{equation} \tag{ 18b }$$

$$\begin{equation}\mathbf{N}=\mathrm{diag}\left[2\kappa ,2\kappa ,2\gamma \left(2\overline{n}+1\right),2\gamma \left(2\overline{n}+1\right)\right].\end{equation} \tag{ 18c }$$

The optimal cooling performance is then achieved for driving on the lower mechanical sideband, Δ = ω_m, and with a sideband-resolved system, κ ≪ ω_m [65, 66]. This is due to the passive beam-splitter part of the optomechanical interaction being resonantly enhanced while the active two-mode-squeezing part (responsible for depositing energy into the mechanical mode) is far off resonant. The steady-state mechanical occupation, obtained by setting $\dot {\mathbf{\Gamma }}=0$ in equation (18), can be found as ${n}_{\text{f}}=\frac{1}{4}\left({{\Gamma}}_{33}+{{\Gamma}}_{44}-2\right)$ and is plotted figure 1 as the solid blue line.

Zoom In Zoom Out Reset image size

To demonstrate the Floquet formalism, we start from the Hamiltonian (15) driven on the red sideband, Δ = ω_m, and move to the interaction picture with respect to ${\hat{H}}_{0}={\omega }_{\text{m}}\left({\hat{c}}^{{\dagger}}\hat{c}+{\hat{b}}^{{\dagger}}\hat{b}\right)$. Thus, the Hamiltonian becomes

$$\begin{equation}\hat{H}=g\left({\hat{c}}^{{\dagger}}\hat{b}+\hat{c}{\hat{b}}^{{\dagger}}\right)+g\left({\hat{c}}^{{\dagger}}{\hat{b}}^{{\dagger}}\;{\text{e}}^{2\text{i}{\omega }_{\text{m}}t}+\hat{c}\hat{b}\;{\text{e}}^{-2\text{i}{\omega }_{\text{m}}t}\right).\end{equation} \tag{ 19 }$$

By applying the RWA (valid for high mechanical frequencies, ω_m ≫ κ, g), we are left with the beam-splitter Hamiltonian

$$\begin{equation}\hat{H}=g\left({\hat{c}}^{{\dagger}}\hat{b}+\hat{c}{\hat{b}}^{{\dagger}}\right).\end{equation} \tag{ 20 }$$

This results in the drift matrix

$$\begin{equation}\mathbf{A}=\left(\begin{pmatrix}{cccc}\hfill -\kappa \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill g\hfill \\ \hfill 0\hfill & \hfill -\kappa \hfill & \hfill -g\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill g\hfill & \hfill -\gamma \hfill & \hfill 0\hfill \\ \hfill -g\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\gamma \hfill \end{pmatrix}\right),\end{equation} \tag{ 21 }$$

while the diffusion matrix N remains the same as before. The RWA clearly demonstrates the mechanism of cooling—mechanical excitations are swapped to the cavity field from which they leak out through the cavity mirrors—but does not properly estimate the final mechanical occupation since it does not include the effective temperature of the optical bath induced by the two-mode squeezing part of the optomechanical interaction. This is visible in regimes where the cavity damping or the optomechanical coupling start approaching the mechanical frequency, invalidating the assumptions necessary for the RWA (see the black dotted line in figure 1, calculated in the same way as for the full model).

The full rotating-frame Hamiltonian (19) is, however, periodic with frequency 2ω_m so we can include the effect of the counterrotating terms in our simulations by applying the Floquet–Lyapunov approach. We begin by forming the equations of motion for the canonical operators,

$$\begin{equation}{\dot {\hat{q}}}_{1}=-\kappa {\hat{q}}_{1}+g{\hat{p}}_{2}-g{\hat{p}}_{2}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+g{\hat{q}}_{2}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\kappa }\;{\hat{q}}_{1,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 22a }$$

$$\begin{equation}{\dot {\hat{p}}}_{1}=-\kappa {\hat{p}}_{1}-g{\hat{q}}_{2}-g{\hat{q}}_{2}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)-g{\hat{p}}_{2}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\kappa }\;{\hat{p}}_{1,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 22b }$$

$$\begin{equation}{\dot {\hat{q}}}_{2}=-\gamma {\hat{q}}_{2}+g{\hat{p}}_{1}-g{\hat{p}}_{1}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+g{\hat{q}}_{1}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\gamma }\;{\hat{q}}_{2,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 22c }$$

$$\begin{equation}{\dot {\hat{p}}}_{2}=-\gamma {\hat{p}}_{2}-g{\hat{q}}_{1}-g{\hat{q}}_{1}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)-g{\hat{p}}_{1}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\gamma }\;{\hat{p}}_{2,\mathrm{i}\mathrm{n}}.\end{equation} \tag{ 22d }$$

From these equations, we can readily read off the frequency components of the Floquet-space drift matrix,

$$\begin{equation}{\mathbf{A}}^{\left(0\right)}={\mathbf{A}}_{0}+g{\mathbf{J}}_{-},\quad \quad {\mathbf{A}}_{\text{c}}^{\left(1\right)}=-\frac{g}{\sqrt{2}}{\mathbf{J}}_{+},\quad \quad {\mathbf{A}}_{\text{s}}^{\left(1\right)}=\frac{g}{\sqrt{2}}{\mathbf{G}}_{1},\end{equation} \tag{ 23 }$$

where A₀ = diag(−κ, −κ, −γ, −γ) is the part of the drift matrix associated with damping,

$$\begin{equation}{\mathbf{J}}_{{\pm}}=\left(\begin{pmatrix}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\pm}1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\pm}1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{pmatrix}\right),\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad {\mathbf{G}}_{1}=\left(\begin{pmatrix}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{pmatrix}\right).\end{equation} \tag{ 24 }$$

Keeping the first three Brillouin zones (i.e., defining the Floquet vector ${\hat{\mathbf{r}}}_{\mathrm{F}}={\left({\hat{\mathbf{r}}}^{\left(0\right)},{\hat{\mathbf{r}}}_{\text{c}}^{\left(1\right)},{\hat{\mathbf{r}}}_{\text{s}}^{\left(1\right)}\right)}^{\top }$), we have the Lyapunov equation

$$\begin{equation}{\dot {\mathbf{\Gamma }}}_{\mathrm{F}}={\mathbf{A}}_{\mathrm{F}}{\mathbf{\Gamma }}_{\mathrm{F}}+{\mathbf{\Gamma }}_{\mathrm{F}}{\mathbf{A}}_{\mathrm{F}}^{\top }+{\mathbf{N}}_{\mathrm{F}},\end{equation} \tag{ 25a }$$

$$\begin{equation}{\mathbf{A}}_{\mathrm{F}}=\left(\begin{pmatrix}{ccc}\hfill {\mathbf{A}}^{\left(0\right)}\hfill & \hfill {\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill \\ \hfill {\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}\hfill & \hfill -2{\omega }_{\text{m}}{\mathbf{I}}_{2}\hfill \\ \hfill {\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill 2{\omega }_{\text{m}}{\mathbf{I}}_{2}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}\hfill \end{pmatrix}\right),\end{equation} \tag{ 25b }$$

$$\begin{equation}{\mathbf{N}}_{\mathrm{F}}=\mathrm{diag}\left(\mathbf{N},\mathbf{N},\mathbf{N}\right).\end{equation} \tag{ 25c }$$

The steady-state mechanical occupation is calculated from the last two diagonal elements of the covariance block Γ⁽⁰⁾ which describe the mechanical variances. When plotted alongside the lab-frame solution (see the dashed yellow line in figure 1), it show perfect agreement with the full model.

In the previous example of optomechanical cooling, the use of the Floquet method is unnecessary as the system can be represented in a frame where the Hamiltonian is time independent. In the rest of this paper, we focus on more complicated driving schemes for which the optomechanical system remains time-dependent regardless of the reference frame. This is generally the case when the system is driven with multiple fields at different frequencies or, equivalently, using a single drive with a modulated amplitude. We focus on the simplest nontrivial example—one mechanical mode and one cavity field subject to a two-tone drive. We perform a detailed analysis of the mechanical squeezing that can be achieved in this system as originally proposed by Kronwald et al [4] in section 3.1. In section 3.2, we then analyse a recently proposed scheme for generating steady-state mechanical squeezing in levitated systems by a combination of parametric and dissipative squeezing [38].

3.1. Mechanical squeezing with a two-tone drive

Following reference [4], we start from the full optomechanical Hamiltonian under two-tone driving,

$$\begin{equation}\hat{H}={\omega }_{\text{c}}{\hat{c}}^{{\dagger}}\hat{c}+{\omega }_{\text{m}}{\hat{b}}^{{\dagger}}\hat{b}+{g}_{0}{\hat{c}}^{{\dagger}}\hat{c}\left(\hat{b}+{\hat{b}}^{{\dagger}}\right)+\left({\eta }_{+}{\text{e}}^{-\text{i}{\omega }_{+}t}+{\eta }_{-}{\text{e}}^{-\text{i}{\omega }_{-}t}\right){\hat{c}}^{{\dagger}}+\mathrm{H}.\mathrm{c}.,\end{equation} \tag{ 26 }$$

where g₀ is the single-photon coupling strength and η_± are the amplitudes of the drives at frequencies ω_± = ω_c ± ω_m. Under this driving, the cavity field acquires a large classical amplitude ${\alpha }_{+}{\text{e}}^{-\text{i}{\omega }_{+}t}+{\alpha }_{-}{\text{e}}^{-\text{i}{\omega }_{-}t}$; introducing the linearised coupling rates g_± = g₀ α_± (assuming, for simplicity, ${\alpha }_{{\pm}}\in \mathbb{R}$) and moving to the interaction picture with respect to ${\hat{H}}_{0}={\omega }_{-}{\hat{c}}^{{\dagger}}\hat{c}+{\omega }_{\text{m}}{\hat{b}}^{{\dagger}}\hat{b}$, we obtain the interaction Hamiltonian

$$\begin{equation}\hat{H}={g}_{-}\left({\hat{c}}^{{\dagger}}\hat{b}+{\hat{b}}^{{\dagger}}\hat{c}\right)+{g}_{+}\left(\hat{c}\hat{b}+{\hat{c}}^{{\dagger}}{\hat{b}}^{{\dagger}}\right)+{g}_{-}\left(\hat{c}\hat{b}\;{\text{e}}^{-2\text{i}{\omega }_{\text{m}}t}+{\hat{c}}^{{\dagger}}{\hat{b}}^{{\dagger}}{\text{e}}^{2\text{i}{\omega }_{\text{m}}t}\right)+{g}_{+}\left({\hat{c}}^{{\dagger}}\hat{b}{\text{e}}^{-2\text{i}{\omega }_{\text{m}}t}+{\hat{b}}^{{\dagger}}\hat{c}\;{\text{e}}^{2\text{i}{\omega }_{\text{m}}t}\right).\end{equation} \tag{ 27 }$$

To follow the Floquet approach, we first formulate the equations of motion

$$\begin{equation}{\dot {\hat{q}}}_{1}=-\kappa {\hat{q}}_{1}+\left({g}_{-}-{g}_{+}\right){\hat{p}}_{2}-\left({g}_{-}-{g}_{+}\right){\hat{p}}_{2}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+\left({g}_{-}-{g}_{+}\right){\hat{q}}_{2}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\kappa }{\hat{q}}_{1,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 28a }$$

$$\begin{equation}{\dot {\hat{p}}}_{1}=-\kappa {\hat{p}}_{1}-\left({g}_{-}+{g}_{+}\right){\hat{q}}_{2}-\left({g}_{-}+{g}_{+}\right){\hat{q}}_{2}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)-\left({g}_{-}+{g}_{+}\right){\hat{p}}_{2}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\kappa }{\hat{p}}_{1,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 28b }$$

$$\begin{equation}{\dot {\hat{q}}}_{2}=-\gamma {\hat{q}}_{2}+\left({g}_{-}-{g}_{+}\right){\hat{p}}_{1}-\left({g}_{-}-{g}_{+}\right){\hat{p}}_{1}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+\left({g}_{-}+{g}_{+}\right){\hat{q}}_{1}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\gamma }{\hat{q}}_{2,\mathrm{i}\mathrm{n}},\end{equation} \tag{ 28c }$$

$$\begin{equation}{\dot {\hat{p}}}_{2}=-\gamma {\hat{p}}_{2}-\left({g}_{-}+{g}_{+}\right){\hat{q}}_{1}-\left({g}_{-}+{g}_{+}\right){\hat{q}}_{1}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)-\left({g}_{-}-{g}_{+}\right){\hat{p}}_{1}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\sqrt{2\gamma }{\hat{p}}_{2,\mathrm{i}\mathrm{n}}.\end{equation} \tag{ 28d }$$

Moving to the Floquet basis, we obtain the drift matrix

$$\begin{equation}{\mathbf{A}}_{\mathrm{F}}=\left(\begin{matrix}{ccccc}{\mathbf{A}}^{\left(0\right)}& {\mathbf{A}}_{\text{c}}^{\left(1\right)}& {\mathbf{A}}_{\text{s}}^{\left(1\right)}& 0& 0\\ {\mathbf{A}}_{\text{c}}^{\left(1\right)}& {\mathbf{A}}^{\left(0\right)}& -2{\omega }_{\text{m}}{\mathbf{I}}_{2}& \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}& \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}\\ {\mathbf{A}}_{\text{s}}^{\left(1\right)}& 2{\omega }_{\text{m}}{\mathbf{I}}_{2}& {\mathbf{A}}^{\left(0\right)}& -\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}& \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}\\ 0& \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}& -\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}& {\mathbf{A}}^{\left(0\right)}& -4{\omega }_{\text{m}}{\mathbf{I}}_{2}\\ 0& \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}& \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}& 4{\omega }_{\text{m}}{\mathbf{I}}_{2}& {\mathbf{A}}^{\left(0\right)}\end{matrix}\right),\end{equation} \tag{ 29 }$$

where

$$\begin{equation}{\mathbf{A}}^{\left(0\right)}={\mathbf{A}}_{0}+{g}_{-}{\mathbf{J}}_{-}-{g}_{+}{\mathbf{J}}_{+},\quad \quad {\mathbf{A}}_{\text{c}}^{\left(1\right)}=-\frac{{g}_{-}}{\sqrt{2}}{\mathbf{J}}_{+}+\frac{{g}_{+}}{\sqrt{2}}{\mathbf{J}}_{-},\quad \quad {\mathbf{A}}_{\text{s}}^{\left(1\right)}=\frac{{g}_{-}}{\sqrt{2}}{\mathbf{G}}_{1}-\frac{{g}_{+}}{\sqrt{2}}{\mathbf{G}}_{2},\end{equation} \tag{ 30 }$$

where we introduced the matrix

$$\begin{equation}{\mathbf{G}}_{2}=\left(\begin{pmatrix}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{pmatrix}\right);\end{equation} \tag{ 31 }$$

the diffusion matrix (for a single Brillouin zone) is the same as in the case of optomechanical cooling. Under the RWA, where the time-dependent parts in equation (27) are ignored, we are left with just the Brillouin zone n = 0, i.e, the matrix A⁽⁰⁾.

To evaluate the mechanical squeezing we need the submatrix describing the mechanical covariances. The solution of the Lyapunov equation is of the form

$$\begin{equation}{\mathbf{\Gamma }}^{\left(0\right)}=\left(\begin{pmatrix}{cc}\hfill {\mathbf{\Gamma }}_{\text{cav}}\hfill & \hfill {\mathbf{\Gamma }}_{\text{C}}\hfill \\ \hfill {\mathbf{\Gamma }}_{\text{C}}^{\top }\hfill & \hfill {\mathbf{\Gamma }}_{\text{m}}\hfill \end{pmatrix}\right),\end{equation} \tag{ 32 }$$

where Γ⁽⁰⁾ is a 4 × 4 matrix. The steady-state squeezed and antisqueezed mechanical variances are

$$\begin{equation}{V}_{\text{sq}}=\mathrm{min}\;\mathrm{e}\mathrm{i}\mathrm{g}\left({\mathbf{\Gamma }}_{\text{m}}\right),\quad \quad {V}_{\text{asq}}=\mathrm{max}\;\mathrm{e}\mathrm{i}\mathrm{g}\left({\mathbf{\Gamma }}_{\text{m}}\right).\end{equation} \tag{ 33 }$$

The variances simulated by the Floquet–Lyapunov method and under the RWA are shown in figure 2. For small coupling (g_± ≪ ω_m) and a sideband resolved system (κ ≪ ω_m), there is a good agreement between the two approaches. As the coupling or the cavity decay is increased, particularly the squeezed variance is affected by the presence of the counterrotating terms. This result is not surprising—as the counterrotating terms give rise to a phase-independent backaction limit, they add a fixed amount of noise to both mechanical quadratures. As the squeezed variance is smaller, the same amount of added noise corresponds to a larger relative increase of the variance than for the antisqueezed quadrature. Similar to the case of mechanical cooling in figure 1, the Floquet solution becomes unstable for strong coupling, behaviour not visible with the RWA.

Zoom In Zoom Out Reset image size

3.2. Parametric and dissipative squeezing for a levitated particle

In the previous example, only one frequency component (oscillating at ω = 2ω_m) was present in the system. Now, we turn our attention to a system where also the next component is relevant—a levitated particle squeezed by a combination of parametric and dissipative squeezing [38]. The potential for the particle's centre-of-mass motion is defined by the laser beam used for levitation; its scattering into an empty cavity mode provides the optomechanical interaction [56, 61, 67]. To achieve strong squeezing, the optical tweezer amplitude is modulated at twice the mechanical frequency, E_tw(t) = E₀[1 + α cos(2ω_m t)], where α ∈ (0, 1) is the modulation depth, modulating both the mechanical potential and the optomechanical coupling rate. Since the potential is proportional to the tweezer intensity (i.e., the square of the amplitude) and the optomechanical coupling to the tweezer amplitude, the system dynamics are characterised by the Hamiltonian [38]

$$\begin{equation}\hat{H}=\frac{{\omega }_{\text{m}}}{2}{\hat{p}}_{2}^{2}+\frac{{\omega }_{\text{m}}}{2}{\left[1+\alpha \;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)\right]}^{2}{\hat{q}}_{2}^{2}+{\Delta}{\hat{c}}^{{\dagger}}\hat{c}-g\left[1+\alpha \;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)\right]{\hat{q}}_{2}\left({\hat{c}}^{{\dagger}}+\hat{c}\right).\end{equation} \tag{ 34 }$$

The parametric squeezing is provided by the modulation of the trapping frequency according to the second term in equation (34). To see how the dissipative squeezing arises, we set the detuning to the red mechanical sideband, Δ = ω_m, and move to the rotating frame with respect to the free oscillations, ${\hat{H}}_{0}={\omega }_{\text{m}}{\hat{c}}^{{\dagger}}\hat{c}+\frac{1}{2}{\omega }_{\text{m}}\left({\hat{q}}_{2}^{2}+{\hat{p}}_{2}^{2}\right)$; under the RWA, we then obtain the interaction-picture Hamiltonian

$$\begin{equation}{\hat{H}}_{\mathrm{i}\mathrm{n}\mathrm{t}}=\frac{{\omega }_{\text{m}}\alpha }{4}\left({\hat{b}}^{2}+{\hat{b}}^{{\dagger}2}+\alpha {\hat{b}}^{{\dagger}}\hat{b}\right)-\frac{g}{\sqrt{2}}\left(\hat{b}+\frac{\alpha }{2}{\hat{b}}^{{\dagger}}\right){\hat{c}}^{{\dagger}}-\frac{g}{\sqrt{2}}\left({\hat{b}}^{{\dagger}}+\frac{\alpha }{2}\hat{b}\right)\hat{c}.\end{equation} \tag{ 35 }$$

Introducing the mechanical Bogoliubov mode $\hat{\beta }=\left(2\hat{b}+\alpha {\hat{b}}^{{\dagger}}\right)/\sqrt{4-{\alpha }^{2}}$ (which satisfies $\left[\hat{\beta },{\hat{\beta }}^{{\dagger}}\right]=1$), we can rewrite this Hamiltonian as

$$\begin{equation}{\hat{H}}_{\mathrm{i}\mathrm{n}\mathrm{t}}=-\frac{{\omega }_{\text{m}}{\alpha }^{2}}{4}{\hat{\beta }}^{{\dagger}}\hat{\beta }+\frac{{\omega }_{\text{m}}\alpha }{4}\left({\hat{\beta }}^{2}+{\hat{\beta }}^{{\dagger}2}\right)-{\lambda }_{\text{eff}}\left({\hat{c}}^{{\dagger}}\hat{\beta }+{\hat{\beta }}^{{\dagger}}\hat{c}\right),\end{equation} \tag{ 36 }$$

where ${\lambda }_{\text{eff}}=\lambda \sqrt{\left(4-{\alpha }^{2}\right)/8}$. From this expression, we see that the Bogoliubov mode undergoes slow oscillations at frequency ω_m α² ≪ ω_m (the first term on the right-hand side) while being squeezed parametrically (the second term in the Hamiltonian); dissipative squeezing is provided by the optomechanical interaction in the last term. Remarkably, these two squeezing techniques can combine and provide stronger squeezing than either strategy alone [38].

To describe the generated squeezing beyond the RWA, we write the interaction Hamiltonian (35) including the counterrotating terms,

$$\begin{align}\hfill {\hat{H}}_{\mathrm{i}\mathrm{n}\mathrm{t}}& =\frac{{\omega }_{\text{m}}\alpha }{4}\left({\hat{b}}^{2}+{\hat{b}}^{{\dagger}2}+\alpha {\hat{b}}^{{\dagger}}\hat{b}\right)-\frac{g}{\sqrt{2}}\left(\hat{b}+\frac{\alpha }{2}{\hat{b}}^{{\dagger}}\right){\hat{c}}^{{\dagger}}-\frac{g}{\sqrt{2}}\left({\hat{b}}^{{\dagger}}+\frac{\alpha }{2}\hat{b}\right)\hat{c}\hfill \\ \hfill & \quad +\;\left[{\omega }_{\text{m}}\alpha {\hat{b}}^{{\dagger}}\hat{b}+\frac{{\omega }_{\text{m}}{\alpha }^{2}}{8}\left({\hat{b}}^{2}+{\hat{b}}^{{\dagger}2}\right)-\frac{g\alpha }{\sqrt{2}}\left({\hat{c}}^{{\dagger}}\hat{b}+{\hat{b}}^{{\dagger}}\hat{c}\right)\right]\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+\left(\frac{{\omega }_{\text{m}}{\alpha }^{2}}{16}{\hat{b}}^{2}-\frac{g}{\sqrt{2}}\hat{c}\hat{b}\right){\text{e}}^{-2\text{i}{\omega }_{\text{m}}t}\hfill \\ \hfill & \quad +\;\left(\frac{{\omega }_{\text{m}}{\alpha }^{2}}{16}{\hat{b}}^{{\dagger}2}-\frac{g}{\sqrt{2}}{\hat{c}}^{{\dagger}}{\hat{b}}^{{\dagger}}\right){\text{e}}^{2\text{i}{\omega }_{\text{m}}t}+\frac{{\omega }_{\text{m}}{\alpha }^{2}}{4}{\hat{b}}^{{\dagger}}\hat{b}\;\mathrm{cos}\left(4{\omega }_{\text{m}}t\right)+\left(\frac{{\omega }_{\text{m}}\alpha }{4}{\hat{b}}^{2}-\frac{g\alpha }{2\sqrt{2}}\hat{c}\hat{b}\right){\text{e}}^{-4\text{i}{\omega }_{\text{m}}t}\hfill \\ \hfill & \quad +\;\left(\frac{{\omega }_{\text{m}}\alpha }{4}{\hat{b}}^{{\dagger}2}-\frac{g\alpha }{2\sqrt{2}}{\hat{c}}^{{\dagger}}{\hat{b}}^{{\dagger}}\right){\text{e}}^{4\text{i}{\omega }_{\text{m}}t}+\frac{{\omega }_{\text{m}}{\alpha }^{2}}{16}\left({\hat{b}}^{{\dagger}2}{\text{e}}^{6\text{i}{\omega }_{\text{m}}t}+{\hat{b}}^{2}{\text{e}}^{-6\text{i}{\omega }_{\text{m}}t}\right).\hfill \end{align} \tag{ 37 }$$

From this Hamiltonian, we derive the equations of motion

$$\begin{align}\hfill {\dot {\hat{q}}}_{1}& =-\kappa {\hat{q}}_{1}-\frac{g}{\sqrt{2}}\left(\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{sin}\left(4{\omega }_{\text{m}}t\right)\right){\hat{q}}_{2}\hfill \\ \hfill & \quad -\frac{g}{\sqrt{2}}\left(\frac{2-\alpha }{2}-\left(1-\alpha \right)\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)-\frac{\alpha }{2}\;\mathrm{cos}\left(4{\omega }_{\text{m}}t\right)\right){\hat{p}}_{2}+\sqrt{2\kappa }\;{\hat{q}}_{1,\mathrm{i}\mathrm{n}},\hfill \end{align} \tag{ 38a }$$

$$\begin{align}\hfill {\dot {\hat{p}}}_{1}& =-\kappa {\hat{p}}_{1}+\frac{g}{\sqrt{2}}\left(\frac{2+\alpha }{2}+\left(1+\alpha \right)\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{cos}\left(4{\omega }_{\text{m}}t\right)\right){\hat{q}}_{2}\hfill \\ \hfill & \quad +\frac{g}{\sqrt{2}}\left(\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{sin}\left(4{\omega }_{\text{m}}t\right)\right){\hat{p}}_{2}+\sqrt{2\kappa }\;{\hat{p}}_{1,\mathrm{i}\mathrm{n}},\hfill \end{align} \tag{ 38b }$$

$$\begin{align}\hfill {\dot {\hat{q}}}_{2}& =-\frac{g}{\sqrt{2}}\left(\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{sin}\left(4{\omega }_{\text{m}}t\right)\right){\hat{q}}_{1}-\frac{g}{\sqrt{2}}\left(\frac{2-\alpha }{2}-\left(1-\alpha \right)\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)-\frac{\alpha }{2}\;\mathrm{cos}\left(4{\omega }_{\text{m}}t\right)\right){\hat{p}}_{1}\hfill \\ \hfill & \quad -\;\gamma {\hat{q}}_{2}+\frac{{\omega }_{\text{m}}\alpha }{2}\left(\frac{\alpha }{4}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\mathrm{sin}\left(4{\omega }_{\text{m}}t\right)+\frac{\alpha }{4}\;\mathrm{sin}\left(6{\omega }_{\text{m}}t\right)\right){\hat{q}}_{2}\hfill \\ \hfill & \quad -\;\frac{{\omega }_{\text{m}}\alpha }{4}\left(\left(2-\alpha \right)-\frac{8-3\alpha }{2}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+\left(2-\alpha \right)\mathrm{cos}\left(4{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{cos}\left(6{\omega }_{\text{m}}t\right)\right){\hat{p}}_{2}+\sqrt{2\gamma }\;{\hat{q}}_{2,\mathrm{i}\mathrm{n}},\hfill \end{align} \tag{ 38c }$$

$$\begin{align}\hfill {\dot {\hat{p}}}_{2}& =\frac{g}{\sqrt{2}}\left(\frac{2+\alpha }{2}+\left(1+\alpha \right)\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{cos}\left(4{\omega }_{\text{m}}t\right)\right){\hat{q}}_{1}+\frac{g}{\sqrt{2}}\left(\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{sin}\left(4{\omega }_{\text{m}}t\right)\right){\hat{p}}_{1}\hfill \\ \hfill & \quad -\;\frac{{\omega }_{\text{m}}\alpha }{4}\left(\left(2+\alpha \right)+\frac{8+3\alpha }{2}\;\mathrm{cos}\left(2{\omega }_{\text{m}}t\right)+\left(2+\alpha \right)\mathrm{cos}\left(4{\omega }_{\text{m}}t\right)+\frac{\alpha }{2}\;\mathrm{cos}\left(6{\omega }_{\text{m}}t\right)\right){\hat{q}}_{2}\hfill \\ \hfill & \quad -\;\gamma {\hat{p}}_{2}-\frac{{\omega }_{\text{m}}\alpha }{2}\left(\frac{\alpha }{4}\;\mathrm{sin}\left(2{\omega }_{\text{m}}t\right)+\mathrm{sin}\left(4{\omega }_{\text{m}}t\right)+\frac{\alpha }{4}\;\mathrm{sin}\left(6{\omega }_{\text{m}}t\right)\right){\hat{p}}_{2}+\sqrt{2\gamma }\;{\hat{p}}_{2,\mathrm{i}\mathrm{n}},\hfill \end{align} \tag{ 38d }$$

from which we obtain the drift matrix

$$\begin{align}\hfill {\mathbf{A}}_{\mathrm{F}}=\left(\begin{pmatrix}{ccccccc}\hfill {\mathbf{A}}^{\left(0\right)}\hfill & \hfill {\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill {\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill {\mathbf{A}}_{\text{c}}^{\left(3\right)}\hfill & \hfill {\mathbf{A}}_{\text{s}}^{\left(3\right)}\hfill \\ \hfill {\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill -2{\omega }_{\text{m}}{\mathbf{I}}_{2}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}+{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill \\ \hfill {\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill 2{\omega }_{\text{m}}{\mathbf{I}}_{2}+\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}-\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left(-{\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}-{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill -\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill \\ \hfill {\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}+{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left(-{\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill {\mathbf{A}}^{\left(0\right)}\hfill & \hfill -4{\omega }_{\text{m}}{\mathbf{I}}_{2}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill \\ \hfill {\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{s}}^{\left(1\right)}+{\mathbf{A}}_{\text{s}}^{\left(3\right)}\right)\hfill & \hfill \frac{1}{\sqrt{2}}\left({\mathbf{A}}_{\text{c}}^{\left(1\right)}-{\mathbf{A}}_{\text{c}}^{\left(3\right)}\right)\hfill & \hfill 4{\omega }_{\text{m}}{\mathbf{I}}_{2}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}\hfill & \hfill -\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill \\ \hfill {\mathbf{A}}_{\text{c}}^{\left(3\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill -\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill -\frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}\hfill & \hfill -6{\omega }_{\text{m}}{\mathbf{I}}_{2}\hfill \\ \hfill {\mathbf{A}}_{\text{s}}^{\left(3\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(2\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{s}}^{\left(1\right)}\hfill & \hfill \frac{1}{\sqrt{2}}{\mathbf{A}}_{\text{c}}^{\left(1\right)}\hfill & \hfill 6{\omega }_{\text{m}}{\mathbf{I}}_{2}\hfill & \hfill {\mathbf{A}}^{\left(0\right)}\hfill \end{pmatrix}\right),\end{align} \tag{ 39 }$$

where

$$\begin{equation}{\mathbf{A}}^{\left(0\right)}={\mathbf{A}}_{0}-\frac{g}{\sqrt{2}}{\mathbf{J}}_{-}+\frac{g\alpha }{2\sqrt{2}}{\mathbf{J}}_{+}+\frac{{\omega }_{\text{m}}{\alpha }^{2}}{4}{\mathbf{M}}_{-}-\frac{{\omega }_{\text{m}}\alpha }{2}{\mathbf{M}}_{+},\end{equation} \tag{ 40a }$$

$$\begin{equation}{\mathbf{A}}_{\text{c}}^{\left(1\right)}=\frac{g}{2}{\mathbf{J}}_{+}-\frac{g\alpha }{2}{\mathbf{J}}_{-}+\frac{{\omega }_{\text{m}}\alpha }{\sqrt{2}}{\mathbf{M}}_{-}-\frac{3{\omega }_{\text{m}}{\alpha }^{2}}{8\sqrt{2}}{\mathbf{M}}_{+},\quad \quad {\mathbf{A}}_{\text{s}}^{\left(1\right)}=-\frac{g}{2}{\mathbf{G}}_{1}+\frac{{\omega }_{\text{m}}{\alpha }^{2}}{8\sqrt{2}}{\mathbf{M}}_{0},\end{equation} \tag{ 40b }$$

$$\begin{equation}{\mathbf{A}}_{\text{c}}^{\left(2\right)}=\frac{g\alpha }{4}{\mathbf{J}}_{+}-\frac{{\omega }_{\text{m}}\alpha }{2\sqrt{2}}{\mathbf{M}}_{+}+\frac{{\omega }_{\text{m}}{\alpha }^{2}}{4\sqrt{2}}{\mathbf{M}}_{-},\quad \quad {\mathbf{A}}_{\text{s}}^{\left(2\right)}=-\frac{g\alpha }{4}{\mathbf{G}}_{1}+\frac{{\omega }_{\text{m}}\alpha }{2\sqrt{2}}{\mathbf{M}}_{0},\end{equation} \tag{ 40c }$$

$$\begin{equation}{\mathbf{A}}_{\text{c}}^{\left(3\right)}=-\frac{{\omega }_{\text{m}}{\alpha }^{2}}{8\sqrt{2}}{\mathbf{M}}_{+},\quad \quad {\mathbf{A}}_{\text{s}}^{\left(3\right)}=\frac{{\omega }_{\text{m}}{\alpha }^{2}}{8\sqrt{2}}{\mathbf{M}}_{0}\end{equation} \tag{ 40d }$$

with

$$\begin{equation}{\mathbf{M}}_{{\pm}}=\left(\begin{pmatrix}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\pm}1\hfill & \hfill 0\hfill \end{pmatrix}\right),\quad \quad {\mathbf{M}}_{0}=\mathrm{diag}\left(0,0,1,-1\right).\end{equation} \tag{ 41 }$$

Again, the diffusion matrix is the same as before and the RWA corresponds to working with the zeroth Brillouin zone only.

We compare the results of the Floquet–Lyapunov approach and the RWA in figure 3. Again, the squeezed quadrature is more sensitive to the validity of the RWA than the antisqueezed quadrature; generally, this effect seems more pronounced now owing to the larger number of counterrotating terms present in the system. Surprisingly, there exist regimes for the levitated particle where the antisqueezed variance is smaller with the counterrotating terms than under the RWA (see the red lines in figure 3(f)). This result points to a nontrivial role that the counterrotating terms play in this situation, which might be harnessed by optimising the modulation phase in the tweezer amplitude, E_tw(t) = E₀[1 + α cos(2ω_m t + ϕ)] (we considered ϕ = 0 for simplicity here).

Zoom In Zoom Out Reset image size

The dynamical stability of the system is now more complicated than in the case of dissipative squeezing via a two-tone drive. Similar to the previous case, the system becomes unstable for large coupling, which is not captured by the RWA calculation. Additionally, the system is unstable (both under RWA and with the Floquet–Lyapunov approach) also for large modulation depth, where the parametric squeezing effect becomes too strong, and for weak optomechanical coupling. We can understand the latter result by noting that, for g → 0, only the parametric squeezing effect persists; for the modulation depth considered in figure 3, this would correspond to more than 3 dB of squeezing, making the system unstable.

In this manuscript, we presented an efficient method for analysing dynamics of Gaussian systems with periodic Hamiltonians. Our approach is based on deriving equations of motion for the canonical operators in the Floquet space from which a Lyapunov equation for the covariance matrix of the system's Wigner function can be formulated. In this way, we arrive at a time-independent algebraic equation for the steady state of the system starting from a time-dependent one; the price we pay for this—increased size of the Hilbert space—is not prohibitive for numerical simulations. We expect the technique to find wide use in cavity optomechanics and electromechanics where time-dependent coupling rates are commonly used for dissipative state preparation. In these systems, large thermal phonon occupations of the mechanical reservoirs preclude direct simulations of the master equation in the Fock basis.

We demonstrated the use of our techniques on dissipative generation of mechanical squeezing in conventional dispersive and novel levitated optomechanical systems and performed a detailed analysis of the role counterrotating terms play in these systems. We showed the intuitive result that the squeezed quadrature is more sensitive to the validity of the RWA than the antisqueezed quadrature owing to its overall lower amount of noise. Importantly, even though counterrotating terms have stronger effect in the recent proposal of combined parametric and dissipative squeezing for levitated particles [38], sub-vacuum squeezing is still possible with state-of-the-art systems; our simulations show that about 2 dB of squeezing (with about 11 dB of antisqueezing) can be generated with parameters similar to reference [62]. Finally, we note that the presence of counterrotating terms can also affect the dynamical stability of the system (particularly when approaching the ultrastrong coupling regime) which cannot be seen from simple analysis under the RWA.

Our approach can be readily adapted for various other state preparation schemes that rely on multitone driving, such as generation of Gaussian entanglement between two mechanical resonators [17, 18]. Our technique allows a straightforward analysis of dissipative state preparation beyond the RWA and a direct evaluation of approaches to mitigate the undesired effects of the counterrotating terms. One such possible technique is using cavity fields with squeezed input noise. This strategy has already been used to suppress the heating associated with counterrotating terms in sideband cooling [68, 69]; it would be interesting to see whether its benefits persist also in more complex dissipative dynamics such as generation of one- and two-mode mechanical squeezing.

Our strategy could also be used to analyse conditional dynamics of optomechanical systems under backaction-evading measurements [34, 36]. In these schemes, driving on both mechanical sidebands is used to perform a quantum nondemolition measurement of one of the mechanical quadratures; similar to dissipative state preparation, the effects of counterrotating terms are often neglected under the RWA. Since the measurement at the output is homodyne (projecting the system onto a Gaussian state), the system can still be fully characterised using the covariance matrix formalism; the important distinction is that the dynamics of the conditional state obey an algebraic equation of the Riccati type [22]. The Riccati equation is closely related to the Lyapunov equation—it operates with the same drift and diffusion matrices (which account for the unconditional part of the evolution) with additional, nonlinear terms stemming from the weak measurement. Since the measurement is typically time-independent in backaction-evading measurements, extension of our formalism to Riccati-type equations should be straightforward.

All in all, dissipative preparation of nonclassical mechanical states and their tomography using quantum nondemolition measurements are important techniques for modern optomechanical and electromechanical systems. As we have discussed, they are often based on the use of periodic Hamiltonians in which resonant parts give rise to the desired dynamics and fast oscillating terms are neglected under the RWA. The RWA can be invalidated by two factors: nonzero sideband ratio (between the cavity damping and mechanical frequency), which stems from practical limitations in realistic experimental systems, and ultrastrong coupling, associated with improved experimental design. With our techniques, the dynamics of these systems can be efficiently analysed including the effects of the fast oscillating terms, allowing the new physics associated with these processes to be uncovered.

Our work was supported by the project LTAUSA19099 from the Czech Ministry of Education, Youth and Sports (MEYS ČR), European Union's Horizon 2020 (2014–2020) research and innovation framework programme under Grant Agreement No. 731473 (project 8C18003 TheBlinQC) (IP and RF), and the project CZ.02.1.01/0.0/0.0/16_026/0008460 of MEYS ČR (OČ and RF). RF was also supported by the Project 20-16577S of the Czech Science Foundation. Project TheBlinQC has received funding from the QuantERA ERA-NET Cofund in Quantum Technologies implemented within the European Union's Horizon 2020 Programme.

According to the Floquet theory, the solution to a linear system $\dot {\mathbf{y}}=\mathbf{B}\left(t\right)\mathbf{y}$, where B(t) is periodic, can be written in the form ${\dot {\mathbf{y}}}_{\mathrm{F}}={\mathbf{B}}_{\mathrm{F}}{\mathbf{y}}_{\mathrm{F}}$, where B_F is time independent. The vectors transform as y = Q(t)y_F, where Q(t) is periodic; we will show in this section that the time-independent matrix B_F can be found from

$$\begin{equation}{\mathbf{Q}\mathbf{B}}_{\mathrm{F}}=\mathbf{B}\mathbf{Q}-\dot {\mathbf{Q}}.\end{equation} \tag{ A1 }$$

Let us consider a time-dependent Hamiltonian that is τ-periodic, Ĥ(t) = Ĥ(t + τ), and quadratic in the canonical operators. Representing all modes with creation (${\hat{a}}_{i}^{{\dagger}}$) and annihilation (â_i ) operators, we can write the time development of their expectation values ${\alpha }_{i}=\left\langle {\hat{a}}_{i}\right\rangle ,{\alpha }_{i}^{{\ast}}=\left\langle {\hat{a}}_{i}^{{\dagger}}\right\rangle$ as

$$\begin{equation}\dot {\mathbf{y}}=\mathbf{B}\left(t\right)\mathbf{y},\end{equation} \tag{ A2 }$$

where we defined the vector $\mathbf{y}={\left({\alpha }_{1},{\alpha }_{1}^{{\ast}},\dots ,\;{\alpha }_{N},{\alpha }_{N}^{{\ast}}\right)}^{\top }\;$. Since the linear system (A2) is τ-periodic with the corresponding frequency ω = 2π/τ, we can decompose y and B(t) in their frequency components according to

$$\begin{equation}{\alpha }_{i}=\sum _{n=-\infty }^{\infty }{\alpha }_{i}^{\left(n\right)}\;{\text{e}}^{\text{i}n\omega \mathrm{t}},\quad {\alpha }_{i}^{{\ast}}=\sum _{n=-\infty }^{\infty }{\alpha }_{i}^{{\ast}\left(n\right)}\;{\text{e}}^{\text{i}n\omega \mathrm{t}},\quad \mathbf{B}\left(t\right)=\sum _{n=-\infty }^{\infty }{\mathbf{B}}^{\left(n\right)}\;{\text{e}}^{\text{i}n\omega \mathrm{t}}.\end{equation} \tag{ A3 }$$

Defining the vector of frequency components ${\mathbf{y}}^{\left(n\right)}={\left({\alpha }_{1}^{\left(n\right)},{\alpha }_{1}^{\left(n\right){\ast}},\dots ,\;{\alpha }_{N}^{\left(n\right)},{\alpha }_{N}^{\left(n\right){\ast}}\right)}^{\top }$ and the Floquet-space vector ${\mathbf{y}}_{\mathrm{F}}={\left(\dots ,\;{\mathbf{y}}^{\left(n\right)},{\mathbf{y}}^{\left(n+1\right)},\dots \right)}^{\top }$, we can write the relation y = Q(t)y_F with the 2N × ∞ matrix

$$\begin{equation}\mathbf{Q}\left(t\right)=\left(\begin{pmatrix}{c}\hfill \dots ,\;{\text{e}}^{\text{i}n\omega \mathrm{t}}{\mathbf{I}}_{N},{\text{e}}^{\text{i}\left(n+1\right)\omega \mathrm{t}}{\mathbf{I}}_{N},\dots \hfill \end{pmatrix}\right),\end{equation} \tag{ A4 }$$

where I_N is the N-mode (and hence 2N-dimensional) identity matrix. With the frequency decomposition of B(t) and the matrix Q(t), we can evaluate equation (A1); a straightforward calculation reveals that individual frequency components obey the differential equation

$$\begin{equation}{\dot {\mathbf{y}}}^{\left(n\right)}=\sum _{m}\left({\mathbf{B}}^{\left(m\right)}-\text{i}n\omega {\delta }_{m0}\right){\mathbf{y}}^{\left(n-m\right)}.\end{equation} \tag{ A5 }$$

Transformation from the frequency components of the creation and annihilation operators (defined in terms of the complex exponential e^inωt ) to the sine and cosine components of the quadratures (as introduced in equation (8)) is possible by introducing the definitions

$$\begin{equation}{q}_{i}^{\left(0\right)}=\frac{{\alpha }_{i}^{\left(0\right)}+{\alpha }_{i}^{{\ast}\left(0\right)}}{\sqrt{2}},\quad \quad {p}_{i}^{\left(0\right)}=-\text{i}\frac{{\alpha }_{i}^{\left(0\right)}-{\alpha }_{i}^{{\ast}\left(0\right)}}{\sqrt{2}},\end{equation} \tag{ A6a }$$

$$\begin{equation}{q}_{i,\mathrm{c}}^{\left(n\right)}=\frac{{\alpha }_{i}^{\left(n\right)}+{\alpha }_{i}^{\left(-n\right)}+{\alpha }_{i}^{{\ast}\left(n\right)}+{\alpha }_{i}^{{\ast}\left(-n\right)}}{2},\quad \quad {p}_{i,\mathrm{c}}^{\left(n\right)}=-\text{i}\frac{{\alpha }_{i}^{\left(n\right)}+{\alpha }_{i}^{\left(-n\right)}-{\alpha }_{i}^{{\ast}\left(n\right)}-{\alpha }_{i}^{{\ast}\left(-n\right)}}{2},\end{equation} \tag{ A6b }$$

$$\begin{equation}{q}_{i,\mathrm{s}}^{\left(n\right)}=\text{i}\frac{{\alpha }_{i}^{\left(n\right)}-{\alpha }_{i}^{\left(-n\right)}+{\alpha }_{i}^{{\ast}\left(n\right)}-{\alpha }_{i}^{{\ast}\left(-n\right)}}{2},\quad \quad {p}_{i,\mathrm{s}}^{\left(n\right)}=\frac{{\alpha }_{i}^{\left(n\right)}-{\alpha }_{i}^{\left(-n\right)}-{\alpha }_{i}^{{\ast}\left(n\right)}+{\alpha }_{i}^{{\ast}\left(-n\right)}}{2}.\end{equation} \tag{ A6c }$$

Note that these quantities are real since, owing to the definitions of the frequency components in equation (A3), we have ${\left({\alpha }_{i}^{\left(n\right)}\right)}^{{\ast}}={\alpha }_{i}^{{\ast}\left(-n\right)}$. Moreover, the normalisation of the operators is chosen such that the basis change is unitary; this choice also fixes the prefactor in front of the sum in the Fourier transform

$$\begin{equation}f\left(t\right)={f}^{\left(0\right)}+\sqrt{2}\sum _{n=1}^{\infty }\left[{f}_{\text{c}}^{\left(n\right)}\;\mathrm{cos}\left(n\omega t\right)+{f}_{\text{s}}^{\left(n\right)}\;\mathrm{sin}\left(n\omega t\right)\right].\end{equation} \tag{ A7 }$$

Indeed, introducing the vectors ${\mathbf{r}}^{\left(0\right)}={\left({q}_{1}^{\left(0\right)},{p}_{1}^{\left(0\right)},\dots ,\;{q}_{N}^{\left(0\right)},{p}_{N}^{\left(0\right)}\right)}^{\top }$, ${\mathbf{r}}_{\text{c}}^{\left(n\right)}={\left({q}_{1,\mathrm{c}}^{\left(n\right)},{p}_{1,\mathrm{c}}^{\left(n\right)},\dots ,\;{q}_{N,\mathrm{c}}^{\left(n\right)},{p}_{N,\mathrm{c}}^{\left(n\right)}\right)}^{\top }$, and ${\mathbf{r}}_{\text{s}}^{\left(n\right)}={\left({q}_{1,\mathrm{s}}^{\left(n\right)},{p}_{1,\mathrm{s}}^{\left(n\right)},\dots ,\;{q}_{N,\mathrm{s}}^{\left(n\right)},{p}_{N,\mathrm{s}}^{\left(n\right)}\right)}^{\top }$, as well as the Floquet vector ${\mathbf{r}}_{\mathrm{F}}={\left({\mathbf{r}}^{\left(0\right)},{\mathbf{r}}_{\text{c}}^{\left(1\right)},{\mathbf{r}}_{\text{s}}^{\left(1\right)},\dots \right)}^{\top }$, we can write

$$\begin{equation}{\mathbf{r}}_{\mathrm{F}}={\mathbf{R}\mathbf{y}}_{\mathrm{F}},\end{equation} \tag{ A8 }$$

where we introduced the unitary matrix

$$\begin{equation}\mathbf{R}=\left(\begin{pmatrix}{cccccc}\hfill \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_{0}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_{\text{c}}\hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_{\text{c}}\hfill & \hfill 0\hfill \\ \hfill \hfill & \hfill 0\hfill & \hfill -{\mathbf{R}}_{\text{s}}\hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_{\text{s}}\hfill & \hfill 0\hfill \\ \hfill \hfill & \hfill {\mathbf{R}}_{\text{c}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_{\text{c}}\hfill \\ \hfill \hfill & \hfill -{\mathbf{R}}_{\text{s}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_{\text{s}}\hfill \\ \hfill \unicode{8944} \hfill & \hfill \hfill & \hfill \hfill & \hfill {\vdots}\hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{pmatrix}\right),\end{equation} \tag{ A9a }$$

$$\begin{align}\hfill {\mathbf{R}}_{0}& =\frac{1}{\sqrt{2}}\left(\begin{pmatrix}{cc}\hfill 1\hfill & \hfill 1\hfill \\ \hfill -\text{i}\hfill & \hfill \text{i}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill \hfill & \hfill \hfill & \hfill -\text{i}\hfill & \hfill \text{i}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{pmatrix}\right),\quad \quad {\mathbf{R}}_{\text{c}}=\frac{1}{2}\left(\begin{pmatrix}{cc}\hfill 1\hfill & \hfill 1\hfill \\ \hfill -\text{i}\hfill & \hfill \text{i}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill \hfill & \hfill \hfill & \hfill -\text{i}\hfill & \hfill \text{i}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{pmatrix}\right)=\frac{1}{\sqrt{2}}{\mathbf{R}}_{0},\hfill \\ \hfill & \hfill \quad {\mathbf{R}}_{\text{s}}=\frac{1}{2}\left(\begin{pmatrix}{cc}\hfill \text{i}\hfill & \hfill \text{i}\hfill \\ \hfill 1\hfill & \hfill -1\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \text{i}\hfill & \hfill \text{i}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill -1\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{pmatrix}\right)=\frac{\mathrm{i}}{\sqrt{2}}{\mathbf{R}}_{0}.\end{align} \tag{ A9b }$$

We can now transform the Langevin equation in the Floquet space, ${\dot {\mathbf{y}}}_{\mathrm{F}}={\mathbf{B}}_{\mathrm{F}}{\mathbf{y}}_{\mathrm{F}}$, into

$$\begin{equation}{\dot {\mathbf{r}}}_{\mathrm{F}}=\mathbf{R}{\dot {\mathbf{y}}}_{\mathrm{F}}={\mathbf{R}\mathbf{B}}_{\mathrm{F}}{\mathbf{R}}^{{\dagger}}{\mathbf{R}\mathbf{y}}_{\mathrm{F}}={\mathbf{A}}_{\mathrm{F}}{\mathbf{r}}_{\mathrm{F}},\end{equation} \tag{ A10 }$$

where A_F is given by equation (11) in the main text.

So far we assumed that the system is homogeneous, which is valid only for the classical expectation values ${\alpha }_{i}=\left\langle {\hat{a}}_{i}\right\rangle$, ${q}_{i}=\left\langle {\hat{q}}_{i}\right\rangle$, ${p}_{i}=\left\langle {\hat{p}}_{i}\right\rangle$. To include quantum fluctuations, we append the original linear system (A2) with a vector of noise operators ${\hat{\xi }}_{y}={\left({\hat{a}}_{1,\mathrm{i}\mathrm{n}},{\hat{a}}_{1,\mathrm{i}\mathrm{n}}^{{\dagger}},\dots ,\;{\hat{a}}_{N,\mathrm{i}\mathrm{n}},{\hat{a}}_{N,\mathrm{i}\mathrm{n}}^{{\dagger}}\right)}^{\top }$,

$$\begin{equation}\dot {\hat{\mathbf{y}}}=\mathbf{B}\left(t\right)\hat{\mathbf{y}}+{\hat{\xi }}_{y},\end{equation} \tag{ A11 }$$

where we now work with the vector of operators $\hat{\mathbf{y}}={\left({\hat{a}}_{1},{\hat{a}}_{1}^{{\dagger}},\dots ,\;{\hat{a}}_{N},{\hat{a}}_{N}^{{\dagger}}\right)}^{\top }$ instead of their expectation values $\mathbf{y}=\left\langle \hat{\mathbf{y}}\right\rangle$. Since the system is Markovian, the noise operators have flat spectra and can thus be separated into their individual frequency components without changing their correlation properties. This allows us to proceed in full analogy to the classical expectation values and write a quantum version of equation (A5) in the form

$$\begin{align}\hfill {\dot {\hat{\mathbf{y}}}}^{\left(n\right)}=\sum _{m}\left({\mathbf{B}}^{\left(m\right)}-\text{i}n\omega {\delta }_{m0}\right){\hat{\mathbf{y}}}^{\left(n-m\right)}+{\hat{\xi }}_{y}^{\left(n\right)}.\end{align} \tag{ A12 }$$

Transformation from annihilation and creation operators to canonical operators then gives the time-independent Langevin equation

$$\begin{align}\hfill {\dot {\hat{\mathbf{r}}}}_{\mathrm{F}}={\mathbf{A}}_{\mathrm{F}}{\hat{\mathbf{r}}}_{\mathrm{F}}+{\hat{\xi }}_{\mathrm{F}}\end{align} \tag{ A13 }$$

with ${\hat{\xi }}_{\mathrm{F}}$ defined in full analogy with ${\hat{\mathbf{r}}}_{\mathrm{F}}$. Crucially, Markovianity ensures that the individual frequency components of the noise are uncorrelated so we can write

$$\begin{equation}{\mathbf{N}}^{\left(0\right)}=\left\langle {\hat{\xi }}^{\left(0\right)}{\hat{\xi }}^{\left(0\right)\top }\right\rangle ={\mathbf{N}}_{\text{k}}^{\left(n\right)}=\left\langle {\hat{\xi }}_{\text{k}}^{\left(n\right)}{\hat{\xi }}_{\text{k}}^{\left(n\right)\top }\right\rangle =\left\langle \hat{\xi }{\hat{\xi }}^{\top }\right\rangle =\mathbf{N},\end{equation} \tag{ A14 }$$

where k ∈ {c, s} and $\hat{\xi }$ is the input quadrature noise of the original system, and N_F = diag(..., N, N, ...). The Floquet-space covariance matrix Γ_F then obeys the Lyapunov equation

$$\begin{equation}{\dot {\mathbf{\Gamma }}}_{\mathrm{F}}={\mathbf{A}}_{\mathrm{F}}{\mathbf{\Gamma }}_{\mathrm{F}}+{\mathbf{\Gamma }}_{\mathrm{F}}{\mathbf{A}}_{\mathrm{F}}^{\top }+{\mathbf{N}}_{\mathrm{F}}.\end{equation} \tag{ A15 }$$

The Floquet approach reproduces the full time-dependent system exactly only when infinitely many frequency components are included in the simulation. In practice, a suitable cutoff frequency has to be chosen in such a way that the inclusion of additional Brillouin zones does not change the result. For the simulations in this manuscript, we confirm the speed of this convergence with respect to the number of Brillouin zones by repeating the calculations with different number of Brillouin zones. The results of this approach can be found in figure 4 where we plot the squeezed variances for dissipative squeezing via a two-tone drive (panel (a)) and for the levitated particle (panel (b)) with different numbers of Brillouin zones.

Zoom In Zoom Out Reset image size

Considering only the first Brillouin zone (the frequency component n = 0; the dotted black lines in figure 4) corresponds to the RWA. For the two-tone driven dissipative squeezing (figure 4(a)), adding the next two Brillouin zones (corresponding to the sine and cosine components for the frequency component n = 1) already seems sufficient as this limit agrees perfectly with the next level of the approximation (i.e., adding the two components for the frequency n = 2). For the simulations in section 3.1, we therefore use five Brillouin zones (frequency components n ⩽ 2) as shown explicitly in the drift matrix (29).

With the joint dissipative and parametric generation for a levitated particle, more Brillouin zones have to be considered as the interaction Hamiltonian includes higher frequency components; here, agreement is reached one step later in the approximation between the cases of n ⩽ 2 and n ⩽ 3 as shown in figure 4(b). For the simulations in section 3.2, we use the first seven Brillouin zones (frequencies n ⩽ 3), corresponding to the drift matrix (39).