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

In this work we consider the two-mode Hamiltonian

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{main:time:independent:Hamiltonian:to:decouple:dimensionful} \hat {H} = \hat{H}_0+ \hbar \,\mathcal{D}_1(t)\,\left(\hat{b}+ \hat{b}{}^{\dagger}\right) + \hbar \, \mathcal{D}_2(t)\left(\hat{b}+ \hat{b}{}^{\dagger} \right)^2 - \hbar\,\mathcal{G}(t) \hat{a}^{\dagger}\hat{a} \, \left(\hat{b}+ \hat{b}{}^{\dagger}\right), \nonumber \end{align} \tag{ 1 }$$

where $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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.

Zoom In Zoom Out Reset image size

The Hamiltonian (1) reduces to the standard optomechanical Hamiltonian when $\mathcal{D}_1=\mathcal{D}_2=0$⁹. 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].

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

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{main:time:independent:Hamiltonian:to:decouple} \frac{\hat {H}}{\hbar \, \omega_{\mathrm{m}}} = \, \Omega_{\mathrm{c}}\,\hat{a}^{\dagger}\hat{a}+\hat{b}^{\dagger}\hat{b}+ \tilde{\mathcal{D}}_1(\tau)\,\left(\hat{b}+ \hat{b}^{\dagger}\right)+ \tilde{\mathcal{D}}_2(\tau)\left(\hat{b}+ \hat{b}^{\dagger} \right)^2 - \tilde{\mathcal{G}}(\tau)\,\hat{a}^{\dagger}\hat{a} \,\left(\hat{b}+ \hat{b}^{\dagger}\right), \nonumber \end{align} \tag{ 2 }$$

which is the Hamiltonian that we will be working with.

When solving the dynamics, we employ methods from the continuous variable formalism [3, 43]. 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 [3, 43], 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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\{\hat{a}_n,\hat{a}^{{\dagger}}_n\}$ evolves to a final state $\hat{\rho}_f$ through the standard Schrödinger equation $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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:={{\rm Tr}}(\hat{\mathcal{A}}\,\hat{\rho})$.

In this language, the canonical commutation relations read $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}[\hat{X}_n,\hat{X}_m^{{\dagger}}]={\rm i}\,\Omega_{nm}$, where the $2N\times2N$ 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\times2N$ symplectic matrix $\boldsymbol{S}$ that preserves the symplectic form, i.e. $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\boldsymbol{S}^{{\dagger}}\,\boldsymbol{\Omega}\,\boldsymbol{S}=\boldsymbol{S}\,\boldsymbol{\Omega}\,\boldsymbol{S}^{{\dagger}}=\boldsymbol{\Omega}$¹⁰. 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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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\times2N$ matrix multiplication. Here, the indices i and f  denote the initial and final state, respectively. Finally, Williamson's theorem guarantees that any $2N\times2N$ Hermitian matrix, such as the covariance matrix $\boldsymbol{\sigma}$, can be decomposed as $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\boldsymbol{\sigma}=\boldsymbol{S}^{{\dagger}}\,\boldsymbol{\nu}_{\oplus}\,\boldsymbol{S}$, where $\boldsymbol{S}$ is an appropriate symplectic matrix. The diagonal matrix $\boldsymbol{\nu}_{\oplus}={\rm 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_{\rm B}\,T})\geqslant1$ (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}_{\rm th}$ of a N-mode bosonic system is simply given by its Williamson form, i.e. $\boldsymbol{\sigma}_{\rm th}=\boldsymbol{\nu}_{\oplus}$.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{general:time:evolution:operator} \hat{U}(t)=\overset{\leftarrow}{\mathcal{T}}\,\exp\left[-\frac{{\rm i}}{\hbar}\int_0^{t} {\rm d}t'\,\hat{H}(t')\right], \nonumber \end{align} \tag{ 3 }$$

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 $\newcommand{\e}{{\rm e}} \hat{U}(t)=\exp[-\frac{{\rm 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]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{decoupled:time:evolution:operator} \hat{U}(t)=\prod_n\hat{U}_n(t)= \prod_n \exp[-{\rm i}\,F_n(t)\,\hat{G}_n], \nonumber \end{align} \tag{ 4 }$$

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 appendix 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]¹¹. 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\times2N$ 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{general:symplectic:decomposition} \boldsymbol{S}=\prod_{n=1}^{N\,(2N+1)}\boldsymbol{S}_n, \nonumber \end{align} \tag{ 5 }$$

where the symplectic matrices $\boldsymbol{S}_n$ are given by $\newcommand{\e}{{\rm e}} \boldsymbol{S}_n:=\exp[F_n(t)\,\boldsymbol{\Omega}\,\boldsymbol{G}_n]$ and the matrices $\boldsymbol{G}_n$ can be obtained through $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 appendix A.

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

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{basis:operator:Lie:algebra} \hat{N}^2_a &:= (\hat{a}^{\dagger} \hat{a})^2 ~~~~\qquad \hat{N}_a := \hat{a}^{\dagger} \hat{a} \qquad\qquad~ \hat{N}_b := \hat{b}^{\dagger} \hat{b} \nonumber \\ \hat{B}_+ &:= \hat{b}^{\dagger} +\hat{b} ~~~~\qquad \hat{B}_- := {\rm i}\,(\hat{b}^{\dagger} -\hat{b}) \nonumber \\ \hat{B}^{(2)}_+ &:= \hat{b}^{\dagger2}+\hat{b}^2~~\qquad \hat{B}^{(2)}_- := {\rm i}\,(\hat{b}^{\dagger2}-\hat{b}^2) \nonumber \\ \hat{N}_a\,\hat{B}_+ &:= \hat{N}_a\,(\hat{b}^{\dagger}+\hat{b}) \quad \hat{N}_a\,\hat{B}_- := \hat{N}_a\,{\rm i}\,(\hat{b}^{\dagger}-\hat{b}), \nonumber \end{align} \tag{ 6 }$$

which form the smallest set of operators in the Lie algebra that generate the Hamiltonian (2)¹².

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 (6) 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.

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:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{U} \hat{U}(\tau)&:={\rm e}^{-{\rm i}\,\Omega_{\rm c} \hat{N}_a\,\tau}\,\hat{U}_{\mathrm{sq}}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a}\,\hat{N}_a}\,{\rm e}^{-{\rm i}\,F_{\hat{N}^2_a}\,\hat{N}^2_a}\,{\rm e}^{-{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{B}_-}\,\hat{B}_-}\nonumber \\ &\quad\times\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_-}\,\hat{N}_a\,\hat{B}_-}, \nonumber \end{align} \tag{ 7 }$$

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

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{decoupling:form:to:be:used} \hat{U}_{\mathrm{sq}} = \overleftarrow{T} \exp\biggl[ - {\rm i} \int^\tau_0 {\rm d}\,\tau' \, 2\,\left(\frac{1}{2} + \tilde{\mathcal{D}}_2(\tau')\right)\,\hat{N}_b+ \tilde{\mathcal{D}}_2(\tau')\,\hat{B}^{(2)}_+ \biggr]. \nonumber \end{align} \tag{ 8 }$$

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 appendix B. We obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sub:algebra:decoupling:solution} F_{\hat{N}_a}&= -2\,\int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau')\int_0^{\tau'}{\rm d}\tau''\,\tilde{\mathcal{G}}(\tau'')\,\Re\xi(\tau'')\, \nonumber \\ &\quad -2 \,\int^\tau_0\,{\rm d}\tau' \,\tilde{\mathcal{G}}(\tau')\, \Im \xi \, \int^{\tau'}_0 \,{\rm d}\tau''\, \tilde{\mathcal{D}}_1(\tau'') \, \Re \xi(\tau'') \, , \nonumber \\ F_{\hat{N}^2_a}&= 2\,\int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{G}}(\tau')\,\Im\xi(\tau')\int_0^{\tau'}{\rm d}\tau''\,\tilde{\mathcal{G}}(\tau'')\,\Re\xi(\tau'') \, , \nonumber \\ F_{\hat{B}_+}&=\int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{D}}_1(\tau')\,\Re\xi(\tau') \, , \nonumber \\ F_{\hat{B}_-}&=- \int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau') \, , \nonumber \\ F_{\hat{N}_a\,\hat{B}_+}&=- \int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{G}}(\tau')\,\Re\xi(\tau') \, , \nonumber \\ F_{\hat{N}_a\,\hat{B}_-}&=\int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{G}}(\tau')\,\Im\xi(\tau'), \nonumber \end{align} \tag{ 9 }$$

where we have introduced the function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:def:xi} \xi(\tau) = P_{11}(\tau)-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,P_{22}(\tau), \nonumber \end{align} \tag{ 10 }$$

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.

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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{b}^{\dagger}$. First of all, we note that a Bogoliubov transformation of a single mode operator always has the general expression $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 appendix B we show that the Bogoliubov coefficients $\alpha(\tau)$ and $\beta(\tau)$ can be obtained through

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:bogoliubov:coefficients} \alpha(\tau)&=\frac{1}{2}\biggl[P_{11}(\tau)+P_{22}(\tau)-{\rm i} \int_0^\tau {\rm d}\tau' P_{22}(\tau') -{\rm i} \int_0^\tau {\rm d}\tau' (1+4\,\tilde{\mathcal{D}}_2(\tau'))\,P_{11}(\tau')\biggr] \, , \nonumber \\ \beta(\tau) &=\frac{1}{2}\biggl[P_{11}(\tau)-P_{22}(\tau)+{\rm i} \int_0^\tau {\rm d}\tau' P_{22}(\tau') -{\rm i} \int_0^\tau {\rm d}\tau' (1+4\,\tilde{\mathcal{D}}_2(\tau'))\,P_{11}(\tau')\biggr] \, , \nonumber \end{align} \tag{ 11 }$$

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) + {\rm i} \dot{\xi}(\tau))$ and $\beta(\tau) = \frac{1}{2}(\xi^*(\tau) + {\rm i} \dot{\xi}^*(\tau))$, where dotted functions imply differentiation with respect to $\tau$.

The two functions P₁₁ and P₂₂ are determined by the two following uncoupled differential equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{differential:equation:written:in:paper} \ddot{P}_{11}+(1+4\,\tilde{\mathcal{D}}_2(\tau))\,P_{11}&=0\nonumber \\ \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}&=0, \nonumber \end{align} \tag{ 12 }$$

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 (12) can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:IP22} \ddot{I}_{P_{22}} + ( 1 + 4 \, \tilde{\mathcal{D}}_2 (\tau)) I_{P_{22}} = 0 \, , \nonumber \end{align} \tag{ 13 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_{P_{22}} = \int^\tau_0 {\rm d}\tau' P_{22} ( \tau') \, . \nonumber \end{align} \tag{ 14 }$$

The solutions to P₁₁ and P₂₂ (or $I_{P_{22}}$) can then be used in the expressions (9)–(11) 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.

In this work, we assume that both the optical and mechanical modes are initially in a coherent states, namely $\renewcommand{\ket}[1]{|#1\rangle} \ket{\mu_{\mathrm{c}}}$ and $\renewcommand{\ket}[1]{|#1\rangle} \ket{\mu_{\mathrm{m}}}$ respectively, defined as the eigenstates of the annihilation operators, i.e. $\renewcommand{\ket}[1]{|#1\rangle} \hat{a} \ket{\mu_{\mathrm{c}}} = \mu_{\mathrm{c}} \ket{\mu_{\mathrm{c}}}$ and $\renewcommand{\ket}[1]{|#1\rangle} \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

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \displaystyle \label{initial:state:two} |\Psi(0)\rangle = \ket{\mu_{\mathrm{c}}}\otimes \ket{\mu_{\mathrm{m}}}. \nonumber \end{align} \tag{ 15 }$$

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

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

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \displaystyle \label{non:linear:state:evolution} \ket{\Psi(\tau)} &= \, {\rm e}^{- {\rm i} \,\left( F_{\hat{B}_+}\, F_{\hat{B}_-}+\Im\left(K \, \mu_{{\rm m}} \right)\right) } \,{\rm e}^{- \frac{1}{2}|\mu_{{\rm c}}|^2}\,\sum_n \biggl[ \frac{\mu_{{\rm c}}^n}{\sqrt{n!}} \, {\rm e}^{-{\rm i} \,\left(F_{\hat{N}^2_a}+ F_{\hat{N}_a \, \hat{B}_+} \, F_{\hat{N}_a \, \hat{B}_-} \right)\, n^2}\nonumber \\ &\,\,\quad{\rm e}^{- {\rm i} \, \left( \Omega_{\rm 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] \, , \nonumber \end{align} \tag{ 16 }$$

where we have defined $K:=F_{\hat{B}_-} + {\rm i} F_{\hat{B}_+}$ and $K_{\hat{N}_a}:=F_{\hat{N}_a \, \hat{B}_- }+ {\rm i} F_{\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 (16) 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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{U}_{{\rm sq}} = {\rm e}^{-{\rm i} \tau \hat{b}^{\dagger} \hat{b}}$ only if $\tilde{\mathcal{D}}_2 \neq 0$.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle \label{non:linear:reduced:mechanical:state:evolution} \hat{\rho}_{{\rm Mech.}}(\tau) = {\rm e}^{- |\mu_{{\rm c}}|^2}\,\sum_n \frac{|\mu_{{\rm c}}|^{2\,n}}{n!} \ket{\phi_n, \tilde{\mathcal{D}}_2(\tau)} \bra{\phi_n, \tilde{\mathcal{D}}_2(\tau)}. \nonumber \end{align} \tag{ 17 }$$

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

The exact expression for $\delta(\tau)$ is long and cumbersome due to the complex expressions of the covariance matrix elements (C.7). 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:araki:lieb} |S(\hat{\rho}_{A})-S(\hat{\rho}_{B})|\leqslant S(\hat{\rho}_{AB})\leqslant S(\hat{\rho}_{A})+S(\hat{\rho}_{B}) \, , \nonumber \end{align} \tag{ 20 }$$

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 appendix 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 = {\rm e}^{- |K_{\hat{N}_a}|^2}$ (see appendix C for a definition of $E_{\hat{B}_+ \hat{B}_-}$ and its appearance in the first and second moments) allow us to find

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{sympelctic:eigenvalue:reduced:state} \nu_{{\rm Op}}^2 &= 1 +4 \, |\mu_{\mathrm{c}}|^2 \left( 1 - {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} \,{\rm e}^{- |K_{\hat{N}_a}|^2} \right) + 4 \, |\mu_{\mathrm{c}}|^4 \biggl( 1 - 2 \, {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} {\rm e}^{- |K_{\hat{N}_a}|^2 }\nonumber \\ &\quad - {\rm e}^{- 4 |\mu_{{\rm c}}|^2 \sin^2\theta} \, {\rm e}^{-4 |K_{\hat{N}_a}|^2 } 2 \, {\rm e}^{- 3|K_{\hat{N}_a}|^2} \, \Re \left\{{\rm e}^{{\rm i} \theta} \, {\rm e}^{|\mu_{\mathrm{c}}|^2 ( {\rm e}^{2{\rm i} \theta} - 1)} \, {\rm e}^{2 |\mu_{{\rm c}}|^2 ({\rm e}^{-{\rm i} \theta} - 1)} \right\} \biggr) \, , \nonumber \\ \nu_{\mathrm{Me}}^2&=1+4\,|K_{\hat{N}_a}|^2 |\mu_{\mathrm{c}}|^2 \, , \nonumber \end{align} \tag{ 21 }$$

where we recall that $K_{\hat{N}_a}:=F_{\hat{N}_a \, \hat{B}_- } +{\rm 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 (21) is bounded by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu_{{\rm Op}} < \sqrt{1 + 4|\mu_{{\rm c}}|^2 + 4 |\mu_{{\rm c}}|^4 } \, , \nonumber \end{align} \tag{ 22 }$$

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 (21) 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 \gg1$, and when $|K_{\hat{N}_a}|^2 \gg1$, which requires $\tilde{g}_0 \gg1$ 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sympelctic:eigenvalue:reduced:state:maximum:value} \nu_{\mathrm{Op,max}}\sim&1+2\,|\mu_{\mathrm{c}}|^2 \, ,\nonumber \\ \nu_{\mathrm{Me,max}}\sim&2\, |K_{\hat{N}_a}| \,|\mu_{\mathrm{c}}| \, . \nonumber \end{align} \tag{ 23 }$$

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} |\ll2|\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

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:many:photons} S(\hat \rho_{\mathrm{G}}(\tau))\sim s_V(1+2\,|\mu_{\mathrm{c}}|^2). \nonumber \end{align} \tag{ 24 }$$
• (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

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:nG:large:KNa} S(\hat \rho_{\mathrm{G}}(\tau))\sim s_V(2\,|K_{\hat{N}_a}|\,|\mu_{\mathrm{c}}|). \nonumber \end{align} \tag{ 25 }$$
• (iii)  
  Finally, when $|K_{\hat{N}_a}|\sim2|\mu_{\mathrm{c}}|$ and $|\mu_{\mathrm{c}}|\gg1$, 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 (9). 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.

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 appendix D. We begin by deriving analytic expressions for the coefficients in (9) given this choice of parameters.

5.2.1. Decoupled dynamics.

We use the methods discussed in section 3 to start by solving the differential equation (12). 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 (11)):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(\tau) &= \frac{1}{2} \left( 2 \cos{\zeta \tau} - \frac{{\rm i}}{\zeta} \left( 1 + \zeta^2 \right)\sin{\zeta \tau} \right), \nonumber \\ \beta(\tau) &= - 2\,{\rm i}\,\frac{\tilde{d}_2 }{\zeta} \sin{\zeta \tau}. \nonumber \end{align} \tag{ 26 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{f:functions:constant:squeezing} F_{\hat{N}_a^2} &= -\frac{\tilde{g}_0^2}{\zeta^2} \left( 1 - {\rm sinc}(2\,\zeta\,\tau) \right)\,\tau, \nonumber \\ F_{\hat{N} _a \, \hat{B}_+} &= - \frac{\tilde{g}_0}{\zeta} \sin{\zeta \tau}, \quad\quad\quad F_{\hat{N} _a \, \hat{B}_-} = \frac{\tilde{g}_0}{\zeta^2}( \cos{\zeta \tau}-1) \, , \nonumber \end{align} \tag{ 27 }$$

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

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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{x}_1 = (\hat{a}^{\dagger} + \hat{a})/\sqrt{2}$ and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{p}_1 = {\rm i} (\hat{a}^{\dagger} - \hat{a})/\sqrt{2}$ and would correspond to classical trajectories in phase space. The full expression for the expectation values $\renewcommand{\bra}[1]{\langle#1|} \renewcommand{\braket}[1]{\langle#1\rangle} \braket{\hat x_1}$ and $\renewcommand{\bra}[1]{\langle#1|} \renewcommand{\braket}[1]{\langle#1\rangle} \braket{\hat p_1}$ can be found in (C.5) in appendix C.

Zoom In Zoom Out Reset image size

In figures 2(a)–(d), 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 2(e)–(h), 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 3(d)) and the lower and upper bounds $\delta_{{\rm min}}$ and $\delta_{{\rm max}}$ provided by the Araki–Lieb inequality in figures 3(e) and (f).

Zoom In Zoom Out Reset image size

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 regimes, as discussed in section 5.1. Using the expression (27) we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |K_{\hat{N}_a}|^2 = \frac{\tilde{g}_0^2}{\zeta^4} \left[ \left( \zeta^2 + 1 \right) \sin^2(\zeta \, \tau) + \cos( 2\, \zeta\, \tau) - 2 \, \cos(\zeta \, \tau) + 1 \right] \, . \nonumber \end{align} \tag{ 28 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \theta(\tau) = 2 \, \frac{\tilde{g}_0^2}{\zeta^3}\left( \sin(\zeta \, \tau) - \zeta \, \tau \right) \, . \nonumber \end{align} \tag{ 29 }$$

We find that $\lim_{\tilde{d}_2 \rightarrow \infty} \theta(\tau) = 0$. We then look at the symplectic eigenvalues (21) 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.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:modulated:squeezing} \tilde{\mathcal{D}}_2(\tau) = \,\tilde{d}_2\, \cos(\Omega_0\, \tau) \, , \nonumber \end{align} \tag{ 30 }$$

where $\tilde{d}_2 = d_2/\omega_{\mathrm{m}}$ is the amplitude of the squeezing and $\Omega_0$ denotes the time-scale of squeezing¹³.

The differential equations in (12) 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₁₁, which we reprint here for convenience,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:reprint:P11} \ddot{P}_{11} + \left( 1 + 4 \, \tilde{d}_2 \cos(\Omega_0 \, \tau) \right) P_{11} = 0 \, . \nonumber \end{align} \tag{ 31 }$$

Equation (31) is that of a parametric oscillator, which is used elsewhere in physics to describe, for example, a driven pendulum. As shown in appendix B, the equation for the integral of P₂₂ (B.19) takes the same form.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Mathieu} \frac{{\rm d}^2 y}{{\rm d}x^2} + \left[ a - 2 \, q \, \cos(2 \, x) \right] y= 0, \nonumber \end{align} \tag{ 32 }$$

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' = \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'$, we rewrite the equation for P₁₁ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\Omega_0^2}{4} \frac{{\rm d}^2 P_{11}}{{\rm d} \tau^{\prime 2}} + \left( 1 + 4 \, \tilde{d}_2 \cos( 2\, \tau') \right) P_{11} = 0 \, , \nonumber \end{align} \tag{ 33 }$$

where we identify the variables $a = 4 /\Omega_0^2$, and $q = -8 \, \tilde{d}_2 / \Omega_0^2$. The boundary conditions P₁₁(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 appendix E), which means that a  =  1 and $q = -2\tilde{d}_2$.

The Mathieu equations are notoriously difficult to evaluate numerically. Instead, we use a two-scale method to derive perturbative solutions to P₁₁ 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 appendix E for the full derivation, where we also show that these approximate solutions correspond exactly to the more physically intuitive RWA when $\tau \gg1$. 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₁₁ and $I_{P_{22}}$ (the integral of P₂₂) are given by, respectively,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{11} &= \cos (\tau) \, \cosh (\tilde{d}_2 \, \tau) - \sin (\tau) \, \sinh ( \tilde{d}_2 \, \tau) \, \nonumber , \nonumber \\ I_{P_{22}} &= - \frac{1}{1- \tilde{d}_2}\left(\cos (\tau) \, \sinh (\tilde{d}_2 \, \tau) - \sin (\tau) \, \cosh (\tilde{d}_2 \, \tau) \right) \, . \nonumber \end{align} \tag{ 34 }$$

We then compute $\xi(\tau)$ in (E.15). We assume that $\tilde{d}_2\tau\ll 1$ to find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:approx:xi} \xi(\tau) \approx \, {\rm e}^{-{\rm i} \tau}\left( 1 + \frac{\tilde{d}_2^2 \tau^2}{2}\right) + {\rm i} \, {\rm e}^{{\rm i} \, \tau} \, \tilde{d}_2 \, \tau \, , \nonumber \end{align} \tag{ 35 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |\alpha(\tau)|^2 - |\beta(\tau)|^2 \approx 1 + \mathcal{O}[(\tilde{d}_2 \tau)^4] \, . \nonumber \end{align} \tag{ 36 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:F:coefficients:modulated:squeezing} F_{\hat{N}_a^2} &= \, \tilde{g}_0^2 \tau\, \left( 1 - \tilde{d}_2 \right) (\mathrm{sinc} (2 \, \tau) - 1) + \frac{1}{2} \tilde{g}_0^2 \, \tilde{d}_2^2 \left(\left(2\,\tau^2 - 3\right) \sin (2 \, \tau) + 2 \, \tau + 4\,\tau \,\cos(2\,\tau)\right)\, ,\nonumber \\ F_{\hat{N}_a \, \hat{B}_+} &= - \tilde{g}_0 \, \sin (\tau) - \tilde{g}_0 \, \tilde{d}_2 \, (\tau \,\cos(\tau) - \sin(\tau)) - \frac{1}{2}\,\tilde{g}_0 \, \tilde{d}_2^2 \left[\left(\tau^2 - 2\right) \sin (\tau) + 2\,\tau \,\cos(\tau)\right] ,\nonumber \\ F_{\hat{N}_a \, \hat{B}_-} &= - 2\,\tilde{g}_0 \,\sin^2(\tau/2) + \tilde{g}_0 \, \tilde{d}_2 \, ( \tau\, \sin(\tau) - 2\,\sin^2(\tau/2)) \nonumber \\ &\quad + \frac{1}{2} \tilde{g}_0 \, \tilde{d}_2^2 \left[\left(\tau^2 - 2\right) \cos (\tau) - 2\,\tau\,\sin(\tau) + 2\right]\, , \nonumber \end{align} \tag{ 37 }$$

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 4(a) 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 4(c) 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.

Zoom In Zoom Out Reset image size

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 (37), we find that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:modulated:K} |K_{\hat{N}_a}|^2 &= \, 4 \, \tilde{g}_0^2 \, \sin ^2(\tau/2) + \tilde{g}_0^2 \, \tilde{d}_2^2 \left(\tau^2 - 2\left(2 - \tau^2\right) \sin ^2(\tau/2)\right) \nonumber \\ &\quad- 2\tilde{g}_0^2 \, \tilde{d}_2 \, \left( \tau \, (\sin (\tau) - \sin (2\,\tau)) + (\cos (\tau) - \cos(2\,\tau)) - 2 \sin ^2(\tau/2) \right), \nonumber \end{align} \tag{ 38 }$$

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.

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 (12), 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 equation (12) 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 (9), which must be evaluated numerically. For each value of $\tau$, the integral will be evaluated from 0 to the final $\tau'$, and then from 0 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.

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

The time evolution operator $\hat{U}(t)$ induced by a Hamiltonian $\hat{H}(t)$ reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{U}(t)=\overset{\leftarrow}{\mathcal{T}}\,\exp\left[-\frac{{\rm i}}{\hbar}\int_0^{t} {\rm d}t'\,\hat{H}(t')\right]. \nonumber \end{align} \tag{ A.1 }$$

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 functions. The choice of $\hat{G}_n$ need not be unique.

It has been shown [28, 47] that it is always possible to obtain the decoupling

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{decoupled:time:evolution:operator:appendix} \hat{U}(t)=\prod_n\hat{U}_n(t), \nonumber \end{align} \tag{ A.2 }$$

where we have defined $\newcommand{\e}{{\rm e}} \hat{U}_n:=\exp[-{\rm i}\,F_n(t)\,\hat{G}_n]$ and the real, time-dependent functions F_n(t), and the ordering of the operators is $\hat{U}_1\hat{U}_2\ldots$.

The functions F_n(t) are uniquely determined by the coupled, nonlinear, first order differential equations

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{operator:differential:equations} \frac{1}{\hbar}\hat{H}=\dot{F}_1\,\hat{G}_1+\dot{F}_2\,\hat{U}_1\,\hat{G}_2\,\hat{U}_1^{{\dagger}}+\dot{F}_3\,\hat{U}_1\,\hat{U}_2\,\hat{G}_3\,\hat{U}_2^{{\dagger}}\,\hat{U}_1^{{\dagger}}+\dot{F}_4\,\hat{U}_1\,\hat{U}_2\,\hat{U}_3\,\hat{G}_4\,\hat{U}_3^{{\dagger}}\,\hat{U}_2^{{\dagger}}\,\hat{U}_1^{{\dagger}}+\dots, \nonumber \end{align} \tag{ A.3 }$$

where $\hat{H}$ and all F_n are taken at time t. This is the general method we will be employing in the following sections.

If the Hamiltonian is quadratic in the mode operators the techniques described in the previous subsection have a more powerful representation. As explained in the main text, the time evolution operator has a symplectic representation $\boldsymbol{S}(t)$ and the decoupled ansatz (A.2) has the form $\boldsymbol{S}(t)=\prod_n\boldsymbol{S}_n(t)$, where we have introduced $\newcommand{\e}{{\rm e}} \boldsymbol{S}_n:=\exp[F_n(t)\,\boldsymbol{\Omega}\,\boldsymbol{G}_n]$ and the real, time-dependent functions F_n(t) are the same as those obtained with the technique above.

The real, time-dependent functions F_n(t) can be obtained by solving the following system of coupled nonlinear first order differential equations

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{matrix:differential:equations} \frac{1}{\hbar}\boldsymbol{H} =\dot{F}_1\,\boldsymbol{G}_1+\dot{F}_2\,\boldsymbol{S}_1^{{\dagger}}\,\boldsymbol{G}_2\,\boldsymbol{S}_1+\dot{F}_3\,\boldsymbol{S}_1^{{\dagger}}\,\boldsymbol{S}_2^{{\dagger}}\,\boldsymbol{G}_3\,\boldsymbol{S}_2\,\boldsymbol{S}_1+\dot{F}_4\,\boldsymbol{S}_1^{{\dagger}}\,\boldsymbol{S}_2^{{\dagger}}\,\boldsymbol{S}_3^{{\dagger}}\,\boldsymbol{G}_4\,\boldsymbol{S}_3\,\boldsymbol{S}_2\,\boldsymbol{S}_1+\dots, \nonumber \end{align} \tag{ A.4 }$$

where the matrix $\boldsymbol{H}$ can be obtained by $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{H}(t)=\frac{1}{2}\mathbb{X}^{{\dagger}}\,\boldsymbol{H}\,\mathbb{X}$ and the summation is over $N\,(2N+1)$ terms [47]. This is the matrix version of the operator differential equation (A.3) for quadratic Hamiltonians, which reduces the problem of operator algebra to matrix multiplication.

We start from the dimensionless Hamiltonian (2), which we reprint here

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \hat {H} = \, \Omega_{\mathrm{c}}\,\hat{a}^{\dagger}\hat{a}+\hat{b}^{\dagger}\hat{b}+ \tilde{\mathcal{D}}_1(\tau)\,\left(\hat{b}+ \hat{b}^{\dagger}\right)+ \tilde{\mathcal{D}}_2(\tau)\left(\hat{b}+ \hat{b}^{\dagger} \right)^2 - \tilde{\mathcal{G}}(\tau)\,\hat{a}^{\dagger}\hat{a} \,\left(\hat{b}+ \hat{b}^{\dagger}\right). \nonumber \end{align} \tag{ B.1 }$$

This Hamiltonian can be conveniently re-written as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{main:time:independent:Hamiltonian:to:decouple:reprinted} \hat {H}&= \, \Omega_{\mathrm{c}} \hat{a}^{\dagger}\hat{a}+ \tilde{\mathcal{D}}_1(\tau) \left(\hat{b}+ \hat{b}^{\dagger}\right)- \tilde{\mathcal{G}}(\tau) \hat{a}^{\dagger}\hat{a} \left(\hat{b}+ \hat{b}^{\dagger}\right) \nonumber \\ &\quad+ 2 \left(\frac{1}{2}+\tilde{\mathcal{D}}_2(\tau)\right) \hat{b}^{\dagger}\hat{b} +\tilde{\mathcal{D}}_2(\tau) \left(\hat{b}^2+ \hat{b}^{\dagger2} \right). \nonumber \end{align} \tag{ B.2 }$$

Our goal now is to separate the Hamiltonian into a linear and a quadratic contribution. We then study the effect of the quadratic Hamiltonian on the nonlinear part, which allows us to focus fully on decoupling the nonlinear dynamics.

We therefore rewrite the time evolution operator (3) in the alternative form

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{formal:solution:to:the:decoupling:easier} \hat{\tilde{U}}(\tau):={\rm e}^{-{\rm i}\,\Omega_\mathrm{c} \hat{a}^{\dagger} \hat{a}\,\tau}\,\hat{U}_{\mathrm{sq}}\,\overset{\leftarrow}{\mathcal{T}}\,\exp\left[-\frac{{\rm i}}{\hbar}\,\int_0^\tau\,{\rm d}\tau'\,\hat{U}_{\mathrm{sq}}^{\dagger}(\tau')\,\hat {H}_1(\tau')\,\hat{U}_{\mathrm{sq}}(\tau')\right], \nonumber \end{align} \tag{ B.3 }$$

where we have introduced

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{app:eq:definition:of:Hsq} \hat {H}_{\mathrm{sq}}(\tau)&:=2\,\left(\frac{1}{2}+\tilde{\mathcal{D}}_2(\tau)\right)\,\hat{b}^{\dagger}\hat{b}+\tilde{\mathcal{D}}_2(\tau)\,\left(\hat{b}^2+ \hat{b}^{\dagger2} \right)\, , \nonumber \\ \hat {H}_1(\tau)&:=\tilde{\mathcal{D}}_1(\tau)\,\left(\hat{b}+ \hat{b}^{\dagger}\right)- \tilde{\mathcal{G}}(\tau)\,\hat{a}^{\dagger}\hat{a} \,\left(\hat{b}+ \hat{b}^{\dagger}\right)\, , \nonumber \\ \hat{U}_{\mathrm{sq}}(\tau)&:=\overset{\leftarrow}{\mathcal{T}}\,\exp\left[-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\hat {H}_{\mathrm{sq}}(\tau')\right]. \nonumber \end{align} \tag{ B.4 }$$

The transition to the expression (B.3) is similar to moving from the Heisenberg (or Schrödinger) picture to the interaction picture.

Our first goal is to solve the time-evolution of the quadratic part $\hat{U}_{{\rm sq}}$. We define the basis vector $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbb{X}:=(\hat{b}, \hat{b}^{\dagger}){}^{\mathrm{T}}$. From standard symplectic (i.e. Bogoliubov) theory we know that

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:Bogoliubov:transform} \mathbb{X}' = \hat{U}_{\mathrm{sq}}^{\dagger} \,\mathbb{X}\, \hat{U}_{\mathrm{sq}} = \left(\begin{array}{@{}c@{}} \hat{U}_{\mathrm{sq}}^{\dagger}\,\hat{b}\,\hat{U}_{\mathrm{sq}} \nonumber \\ \hat{U}_{\mathrm{sq}}^{\dagger}\,\hat{b}^{\dagger}\,\hat{U}_{\mathrm{sq}} \end{array}\right) = \boldsymbol{S}_{\mathrm{sq}}(\tau)\, \mathbb{X}, \nonumber \end{align} \tag{ B.5 }$$

where the $2\times2$ symplectic matrix $\boldsymbol{S}_{\mathrm{sq}}(\tau)$ is the symplectic representation of $\hat {H}_{{\rm sq}}(\tau)$ and satisfies $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\boldsymbol{S}_{\mathrm{sq}}^{\dagger}(\tau)\,\boldsymbol{\Omega}\,\boldsymbol{S}_{\mathrm{sq}}(\tau)=\boldsymbol{\Omega}$. Here $\boldsymbol{\Omega}={{\rm diag}}(-{\rm i},{\rm i})$ is the symplectic form.

The matrix $\boldsymbol{S}_{\mathrm{sq}}(\tau)$ therefore has the expression $\newcommand{\e}{{\rm e}} \boldsymbol{S}_{\mathrm{sq}}(\tau)=\overset{\leftarrow}{\mathcal{T}}\,\exp[\boldsymbol{\, \Omega}\,\int_0^\tau\,{\rm d}\tau'\,\tilde{\boldsymbol{H}}_{{\rm sq}}(\tau')]$. Here $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat {H}_{{\rm sq}}=\frac{1}{2}\mathbb{X}^{\dagger}\tilde{\boldsymbol{H}}_{{\rm sq}} \mathbb{X}$ and one has that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\boldsymbol{H}}_{{\rm sq}} = \left(\begin{array}{@{}cc@{}} 1+2\,\tilde{\mathcal{D}}_2(\tau) & 2\,\tilde{\mathcal{D}}_2(\tau) \nonumber \\ 2\,\tilde{\mathcal{D}}_2(\tau) & 1+2\,\tilde{\mathcal{D}}_2(\tau) \end{array}\right). \nonumber \end{align} \tag{ B.6 }$$

Therefore, we have that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{main:general:expression:time:dependent:squeezing} \boldsymbol{S}_{\mathrm{sq}}(\tau)=\overset{\leftarrow}{\mathcal{T}}\,\exp\left[-{\rm i}\,\left(\begin{array}{@{}cc@{}} 1 & 0 \nonumber \\ 0 & -1 \end{array}\right)\,\int_0^\tau\,{\rm d}\tau'\,\left(\begin{array}{@{}cc@{}} 1+2\,\tilde{\mathcal{D}}_2(\tau') & 2\,\tilde{\mathcal{D}}_2(\tau') \nonumber \\ 2\,\tilde{\mathcal{D}}_2(\tau') & 1+2\,\tilde{\mathcal{D}}_2(\tau') \end{array}\right)\right]. \nonumber \end{align} \tag{ B.7 }$$

It can be shown that the time independent orthogonal matrix $\boldsymbol{M}_{{\rm ort}}$ with expression

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{M}_{{\rm ort}} =\frac{1}{\sqrt{2}} \left(\begin{array}{@{}cc@{}} 1 & 1 \nonumber \\ -1 & 1 \end{array}\right)\, , \nonumber \end{align} \tag{ B.8 }$$

puts the Hamiltonian matrix $\tilde{\boldsymbol{H}}_s$ in diagonal form through

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(\begin{array}{@{}cc@{}} 1+2\,\tilde{\mathcal{D}}_2(\tau) & 2\,\tilde{\mathcal{D}}_2(\tau) \nonumber \\ 2\,\tilde{\mathcal{D}}_2(\tau) & 1+2\,\tilde{\mathcal{D}}_2(\tau) \end{array}\right) = \frac{1}{2} \left(\begin{array}{@{}cc@{}} 1 & -1 \nonumber \\ 1 & 1 \end{array}\right) \left(\begin{array}{@{}cc@{}} 1+4\,\tilde{\mathcal{D}}_2(\tau) & 0 \nonumber \\ 0 &1 \end{array}\right) \left(\begin{array}{@{}cc@{}} 1 & 1 \nonumber \\ -1 & 1 \end{array}\right) . \nonumber \end{align} \tag{ B.9 }$$

This allows us to manipulate (B.7) and find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{main:general:expression:time:dependent:squeezing:manipulated} \boldsymbol{S}_{\mathrm{sq}}(\tau)= \boldsymbol{M}_{{\rm ort}}^{{\rm T}}\,\overset{\leftarrow}{\mathcal{T}}\,\exp\left[{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ 1+4\,\tilde{\mathcal{D}}_2(\tau') & 0 \end{array}\right)\right] \,\boldsymbol{M}_{{\rm ort}}. \nonumber \end{align} \tag{ B.10 }$$

This means that we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbb{X}'= \boldsymbol{M}_{{\rm ort}}^{{\rm T}}\,\overset{\leftarrow}{\mathcal{T}}\,\exp\left[{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ 1+4\,\tilde{\mathcal{D}}_2(\tau') & 0 \end{array}\right)\right] \,\boldsymbol{M}_{{\rm ort}}\, \mathbb{X}. \nonumber \end{align} \tag{ B.11 }$$

We introduce the new vector $\mathbb{Y}:=\boldsymbol{M}_{{\rm ort}}\,\mathbb{X}$, which is just a rotation of the operators. Then we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{matrix:time:ordered:exponential} \mathbb{Y}'= \overset{\leftarrow}{\mathcal{T}}\,\exp\left[{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ 1+4\,\tilde{\mathcal{D}}_2(\tau') & 0 \end{array}\right)\right] \,\mathbb{Y}. \nonumber \end{align} \tag{ B.12 }$$

Here we seek a formal expression for (B.12). We write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{initial:useful:conditions} \overset{\leftarrow}{\mathcal{T}}\,\exp\left[{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ 1+4\,\tilde{\mathcal{D}}_2(\tau') & 0 \end{array}\right)\right]=\boldsymbol{P}+{\rm i}\,\int_0^\tau {\rm d}\tau'\,\boldsymbol{K}\,\boldsymbol{P}, \nonumber \end{align} \tag{ B.13 }$$

in terms of the matrix $\boldsymbol{K}$ defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{matrix:time:ordered:exponential:explicit} \boldsymbol{K}:= \left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ 1+4\,\tilde{\mathcal{D}}_2(\tau') & 0 \end{array}\right), \nonumber \end{align} \tag{ B.14 }$$

and the matrix $\boldsymbol{P}$, which we will determine and which is diagonal. This follows from the fact that the matrix on the left-hand side of (B.13) is diagonal when squared, and therefore any even powers in the expansion will be diagonal. We use the fact that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\rm d}{{\rm d}\tau}\overset{\leftarrow}{\mathcal{T}}\,\exp\left[{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\boldsymbol{K}(\tau')\right] = {\rm i}\,\boldsymbol{K}\, \overset{\leftarrow}{\mathcal{T}}\,\exp\left[{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\boldsymbol{K}(\tau')\right]\, , \nonumber \end{align} \tag{ B.15 }$$

to find the equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -\boldsymbol{K}\,\int_0^\tau {\rm d}\tau'\,\boldsymbol{K}\,\boldsymbol{P}=\dot{\boldsymbol{P}}. \nonumber \end{align} \tag{ B.16 }$$

Since $\boldsymbol{K}$ is invertible, we manipulate this equation and obtain, after some algebra,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{differential:equation:for:matrix} \ddot{\boldsymbol{P}}- 4\,\frac{\dot{\tilde{\mathcal{D}}}_2(\tau)}{1 + 4\,\tilde{\mathcal{D}}_2(\tau)} \left(\begin{array}{@{}cc@{}} 0 & 0 \nonumber \\ 0 & 1 \end{array}\right) \dot{\boldsymbol{P}} +(1+4\,\tilde{\mathcal{D}}_2(\tau))\,\boldsymbol{P}=0. \nonumber \end{align} \tag{ B.17 }$$

We can now solve the four differential equations contained in (B.17) which read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{differential:equation:written:down} \ddot{P}_{11}+(1+4\,\tilde{\mathcal{D}}_2(\tau))\,P_{11}&= \, 0\nonumber \\ P_{12}=P_{21}&= \, 0\nonumber \\ \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}& = \, 0. \nonumber \end{align} \tag{ B.18 }$$

The differential equation (B.18) must be supplemented by initial conditions. We note that since the left hand side of (B.13) is the identity matrix for $\tau=0$ we have that $\boldsymbol{P}(0)=\mathbb{1}$ which implies $P_{11}(0)=P_{22}(0)=1$. In addition, taking the time derivative of both sides of (B.13) and setting $\tau=0$ implies $\dot{P}_{11}(0)=\dot{P}_{22}(0)=0$.

By introducing the integral $I_{P_{22}} = \int^\tau_0 {\rm d} \tau' P_{22} (\tau')$, one can also rewrite the second equation as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:eq:IP22} \ddot{I}_{P_{22}} + ( 1 + 4 \, \tilde{\mathcal{D}}_2 (\tau)) I_{P_{22}} = 0 \, , \nonumber \end{align} \tag{ B.19 }$$

so that it becomes the same as that for P₁₁. The boundary conditions are now $I_{P_{22}}(0) = 0$ and $\dot{I}_{P_{22}} = 1$.

We were not able to find a general solution to the differential equation for P₁₁ (B.18) and the equation for $I_{P_{22}}$ (B.19), but they can be integrated numerically when an explicit form of $\tilde{\mathcal{D}}_2(\tau)$ is given. For specific choices of $\tilde{\mathcal{D}}_2(\tau)$, which we discuss in the main text of this work, the equations become the well-studied Mathieu equation.

Proceeding with the solution to the quadratic time-evolution, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overset{\leftarrow}{\mathcal{T}}\,\exp\left[{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ 1+4\,\tilde{\mathcal{D}}_2(\tau') & 0 \end{array}\right)\right]&=\boldsymbol{P}+{\rm i}\,\int_0^\tau {\rm d}\tau'\,\boldsymbol{K}\,\boldsymbol{P} \nonumber \\ &= \left(\begin{array}{@{}cc@{}} P_{11} & {\rm i}\,\int_0^\tau\,{\rm d}\tau'\,P_{22} \nonumber \\ {\rm i}\,\int_0^\tau\,{\rm d}\tau'\,(1+4\,\tilde{\mathcal{D}}_2(\tau'))\,P_{11} & P_{22} \end{array}\right). \nonumber \end{align} \tag{ B.20 }$$

In turn, this allows us to get

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{S}_{\mathrm{sq}}(\tau) &=\left(\begin{array}{@{}cc@{}} \alpha(\tau) & \beta(\tau) \nonumber \\ \beta^*(\tau) & \alpha^*(\tau) \end{array}\right) \nonumber \\ &=\boldsymbol{M}_{{\rm ort}}^{{\rm T}}\, \left(\begin{array}{@{}cc@{}} P_{11} & {\rm i}\,\int_0^\tau\,{\rm d}\tau'\,P_{22} \nonumber \\ {\rm i}\,\int_0^\tau\,{\rm d}\tau'\,(1+4\,\tilde{\mathcal{D}}_2(\tau'))\,P_{11} & P_{22} \end{array}\right)\, \boldsymbol{M}_{{\rm ort}}, \nonumber \end{align} \tag{ B.21 }$$

where we have introduced the Bogoliubov matrix $\boldsymbol{S}_{\mathrm{sq}}(\tau)$ with coefficients $\alpha(\tau)$ and $\beta(\tau)$.

This gives us, after a little algebra,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(\tau)&=\frac{1}{2}\,\left[P_{11}(\tau)+P_{22}(\tau)-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,P_{22}(\tau')-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,(1+4\,\tilde{\mathcal{D}}_2(\tau'))\,P_{11}(\tau')\right] \, , \nonumber \\ \beta(\tau) &=\frac{1}{2}\,\left[P_{11}(\tau)-P_{22}(\tau)+{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,P_{22}(\tau')-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,(1+4\,\tilde{\mathcal{D}}_2(\tau'))\,P_{11}(\tau')\right]. \nonumber \end{align} \tag{ B.22 }$$

These quantities can also be written in terms of $I_{P_{22}}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(\tau) &= \frac{1}{2}\,\left[ P_{11} - {\rm i} I_{P_{22}} + {\rm i}\frac{\rm d}{{\rm d}\tau} ( P_{11} - {\rm i} I_{P_{22}})\right] \, ,\nonumber \\ \beta(\tau) &= \frac{1}{2}\,\left[ P_{11} + {\rm i} I_{P_{22}} + {\rm i}\frac{\rm d}{{\rm d}\tau} ( P_{11} + {\rm i} I_{P_{22}})\right]. \nonumber \end{align} \tag{ B.23 }$$

This means that the basis vector $\mathbb{X}$ (B.5) transforms as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:Bogoliubov:transform:solved} \mathbb{X}'= \hat{U}_{\mathrm{sq}}^{\dagger} \,\mathbb{X}\, \hat{U}_{\mathrm{sq}} = \left(\begin{array}{@{}cc@{}} \alpha(\tau) & \beta(\tau) \nonumber \\ \beta^*(\tau) & \alpha^*(\tau) \end{array}\right) \,\mathbb{X}. \nonumber \end{align} \tag{ B.24 }$$

As a side remark, the Bogoliubov (symplectic) identities $|\alpha(t)|^2-|\beta(t)|^2=1$ read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{11}\,P_{22}+\left(\int_0^\tau\,{\rm d}\tau'\,P_{22}\right)\,\left(\int_0^\tau\,{\rm d}\tau'\,(1+4\,\tilde{\mathcal{D}}_2(\tau'))\,P_{11}\right)=1. \nonumber \end{align} \tag{ B.25 }$$

We can now go back to the time evolution operator (B.3) which we reprint here

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \hat{\tilde{U}}(\tau)={\rm e}^{-{\rm i}\,\Omega_{\rm c} \hat{a}^{\dagger} \hat{a}\,\tau}\,\hat{U}_{\mathrm{sq}}(\tau)\,\overset{\leftarrow}{\mathcal{T}}\,\exp\left[-\frac{{\rm i}}{\hbar}\,\int_0^\tau\,{\rm d}\tau'\,\hat{U}_{\mathrm{sq}}^{\dagger}(\tau')\,\hat {H}_1(\tau')\,\hat{U}_{\mathrm{sq}}(\tau')\right]. \nonumber \end{align} \tag{ B.26 }$$

Our work above allows to obtain

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \hat{\tilde{U}}(\tau)={\rm e}^{-{\rm i}\,\Omega_{\rm c} \hat{a}^{\dagger} \hat{a}\,\tau}\,\hat{U}_{\mathrm{sq}}(\tau)\,\overset{\leftarrow}{\mathcal{T}}\,&\exp\biggl[-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\bigl(\tilde{\mathcal{D}}_1(\tau') \bigl(\xi(\tau')\,\hat{b}+ \xi^*(\tau')\,\hat{b}{}^{\dagger}\bigr) \nonumber \\ &\quad\quad\quad\quad\quad\quad\quad\quad-\tilde{\mathcal{G}}(\tau') \, \hat{a}^{\dagger}\hat{a}\,\bigl(\xi(\tau')\,\hat{b}+ \xi^*(\tau')\,\hat{b}{}^{\dagger}\bigr)\bigr)\biggr], \nonumber \end{align} \tag{ B.27 }$$

which can be conveniently rewritten as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{almost:final:time:evolution:operator} \hat{\tilde{U}}(t)&={\rm e}^{-{\rm i}\,\Omega_{\rm c} \hat{a}^{\dagger} \hat{a}\,\tau}\,\hat{U}_{\mathrm{sq}}\,\overset{\leftarrow}{\mathcal{T}}\,\exp\biggl[-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\biggl(\tilde{\mathcal{D}}_1(\tau') \Re\xi(\tau')\,\hat{B}_+-{\rm i} \, \tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau')\,\hat{B}_- \nonumber \\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad-\tilde{\mathcal{G}}(\tau') \,\Re\xi(\tau')\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_++{\rm i} \, \tilde{\mathcal{G}}(\tau') \,\Im\xi(\tau')\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_-\biggr)\biggr]. \nonumber \end{align} \tag{ B.28 }$$

Here we have introduced

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:eq:simplified:xi} \xi(\tau):=\alpha(\tau)+\beta^*(\tau)=P_{11}-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,P_{22} \nonumber \end{align} \tag{ B.29 }$$

for conveniency of presentation, where $\Re \xi(\tau)$ and $\Im \xi(\tau)$ are the real and imaginary parts of $\xi(\tau)$. Our definition of $\xi(\tau)$ also implies that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:eq:xi:alpha:beta:relation} \alpha(\tau) = \frac{1}{2}(\xi(\tau) + {\rm i} \dot{\xi}(\tau)) \qquad{{\rm and}} \qquad\beta(\tau) = \frac{1}{2}(\xi^*(\tau) + {\rm i} \dot{\xi}^*(\tau)) . \nonumber \end{align} \tag{ B.30 }$$

We now note that the operators $\hat{N}_a,\hat{N}_a^2, \hat{B}_+,\hat{B}_-,\hat{N}_a\,\hat{B}_+,\hat{N}_a\,\hat{B}_-$ form a closed Lie sub-algebra of the full algebra of our operators. Therefore, we can apply the decoupling techniques described above.

We proceed to decouple the remaining part of the operator (B.28), which reads

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \hat{U}_{\rm remain}(t):=\overset{\leftarrow}{\mathcal{T}}\,&\exp\biggl[-{\rm i}\,\int_0^\tau\,{\rm d}\tau'\,\biggl(\tilde{\mathcal{D}}_1(\tau') \Re\xi(\tau')\,\hat{B}_+-\tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau')\,\hat{B}_- \nonumber \\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad-\tilde{\mathcal{G}}(\tau') \,\Re\xi(\tau')\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_++\tilde{\mathcal{G}}(\tau') \,\Im\xi(\tau')\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_-\biggr)\biggr]. \nonumber \end{align} \tag{ B.31 }$$

We make the ansatz

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{decoupling:ansatz} \hat{U}_{\rm remain}(t) ={\rm e}^{-{\rm i}\,F_{\hat{N}_a}\,\hat{N}_a}\,{\rm e}^{-{\rm i}\,F_{\hat{N}^2_a}\,\hat{N}^2_a}\,{\rm e}^{-{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{B}_-}\,\hat{B}_-}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_-}\,\hat{N}_a\,\hat{B}_-}. \nonumber \end{align} \tag{ B.32 }$$

Taking the time derivative on both sides and arranging in a similar fashion to (A.3) we find

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{main:docupling:Hamiltonian:technique:simplified} &\tilde{\mathcal{D}}_1(\tau) \Re\xi(\tau)\,\hat{B}_+-\tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau)\,\hat{B}_--\tilde{\mathcal{G}}(\tau) \,\Re\xi(\tau)\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_++\tilde{\mathcal{G}}(\tau) \,\Im\xi(\tau)\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_-\nonumber \\ & = \dot F_{\hat{N}_a}\,\hat{N}_a +\dot F_{\hat{N}^2_a}\,\hat{N}^2_a + \dot F_{\hat{B}_+}\,\hat{B}_+ + \dot F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+\nonumber \\ &\quad+ \dot F_{\hat{B}_-}\,{\rm e}^{-{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\, \hat{B}_- \,{\rm e}^{{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\,{\rm e}^{{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\nonumber \\ &\quad+ \dot F_{\hat{N}_a\,\hat{B}_-}\,{\rm e}^{-{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\, \hat{N}_a\,\hat{B}_- \,{\rm e}^{{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\,{\rm e}^{{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\nonumber \\ & = \dot F_{\hat{N}_a}\,\hat{N}_a +\dot F_{\hat{N}^2_a}\,\hat{N}^2_a + \dot F_{\hat{B}_+}\,\hat{B}_+ + \dot F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+\nonumber \\ &\quad+ (\dot F_{\hat{B}_-}+\dot F_{\hat{N}_a\,\hat{B}_-}\,\hat{N}_a)\,{\rm e}^{-{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\, \hat{B}_- \,{\rm e}^{{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\,{\rm e}^{{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\nonumber \\ & = \dot F_{\hat{N}_a}\,\hat{N}_a +\dot F_{\hat{N}^2_a}\,\hat{N}^2_a+ \dot F_{\hat{B}_+}\,\hat{B}_+ + \dot F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+\nonumber \\ &\quad+ (\dot F_{\hat{B}_-} +\dot F_{\hat{N}_a\,\hat{B}_-}\,\hat{N}_a)\,(\hat{B}_-+2\, F_{\hat{B}_+}+2\, F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a). \nonumber \end{align} \tag{ B.33 }$$

Therefore our main differential equations can be obtained by equating the coefficient of the different, linearly independent operators of the Lie algebra from the equation below

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{main:docupling:Hamiltonian:technique:simplifieda} &\tilde{\mathcal{D}}_1(\tau) \Re\xi(\tau)\,\hat{B}_+-\tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau)\,\hat{B}_--\tilde{\mathcal{G}}(\tau) \,\Re\xi(\tau)\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_++\tilde{\mathcal{G}}(\tau) \,\Im\xi(\tau)\,\hat{a}^{\dagger}\hat{a}\,\hat{B}_-\nonumber \\ & = \dot F_{\hat{N}_a}\,\hat{N}_a +\dot F_{\hat{N}^2_a}\,\hat{N}^2_a + \dot F_{\hat{B}_+}\,\hat{B}_+ + \dot F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+\nonumber \\ &\quad+ (\dot F_{\hat{B}_-}+\dot F_{\hat{N}_a\,\hat{B}_-}\,\hat{N}_a)\,(\hat{B}_-+2\, F_{\hat{B}_+}+2\, F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a). \nonumber \end{align} \tag{ B.34 }$$

The solutions with operators proportional to $\hat{B}_\pm$ can be independently solved. However, the solutions for $F_{\hat{N}_a}$, $F_{\hat{N}_a^2}$, and $F_{\hat{N}_a \, \hat{B}_\pm}$ are less straight-forward. We identify the following four coupled differential equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:diff:eqs:explicit} \dot{F}_{\hat{N}_a} &= - 2 \, \dot{F}_{\hat{B}_-} \, F_{\hat{N}_a \, \hat{B}_+} - 2 \, F_{\hat{B}_+} \, \dot{F}_{\hat{N}_a \, \hat{B}_-} \, , \nonumber \\ \dot{F}_{\hat{N}_a^2} &= - 2 \, \dot{F}_{\hat{N}_a \, \hat{B}_-} \, F_{\hat{N}_a \, \hat{B}_+} \, , \nonumber \\ \dot{F}_{\hat{N}_a \, \hat{B}_+} &=- \tilde{\mathcal{G}} ( \tau) \, \Re \xi(\tau) \, , \nonumber \\ \dot{F}_{\hat{N}_a \, \hat{B}_-} &= \tilde{\mathcal{G}} \, \Im \xi(\tau) \, . \nonumber \end{align} \tag{ B.35 }$$

By first solving the equations for $\dot{F}_{\hat{N}_a \, \hat{B}_\pm}$ and $\dot{F}_{\hat{B}_\pm}$, it is then possible to insert the solutions into the expressions for $\dot{F}_{\hat{N}_a}$ and $\dot{F}_{\hat{N}_a^2}$. We find the following key expression for this work

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{appendix:sub:algebra:decoupling:solution} F_{\hat{N}_a}(\tau)&= -2 \,\int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau')\int_0^{\tau'}{\rm d}\tau''\, \tilde{\mathcal{G}}(\tau'')\,\Re\xi(\tau'')\, \nonumber \\ &\quad-2 \, \int^\tau_0\,{\rm d}\tau' \,\tilde{\mathcal{G}}(\tau')\, \Im \xi(\tau') \, \int^{\tau'}_0 \,{\rm d}\tau''\, \tilde{\mathcal{D}}_1(\tau'') \, \Re \xi(\tau'') \, , \nonumber \\ F_{\hat{N}^2_a}(\tau)&= 2 \,\int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{G}}(\tau')\,\Im\xi(\tau')\int_0^{\tau'}{\rm d}\tau''\,\tilde{\mathcal{G}}(\tau'')\,\Re\xi(\tau'')\, , \nonumber \\ F_{\hat{B}_+}(\tau)&= \int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{\mathcal{D}}}_1(\tau')\,\Re\xi(\tau')\, , \nonumber \\ F_{\hat{B}_-}(\tau)&=- \int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{D}}_1(\tau')\,\Im\xi(\tau')\, , \nonumber \\ F_{\hat{N}_a\,\hat{B}_+}(\tau)&=- \int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{G}}(\tau')\,\Re\xi(\tau')\, , \nonumber \\ F_{\hat{N}_a\,\hat{B}_-}(\tau)&=\int_0^\tau\,{\rm d}\tau'\,\tilde{\mathcal{G}}(\tau')\,\Im\xi(\tau') \, , \nonumber \end{align} \tag{ B.36 }$$

where $\Re\xi$ denotes the real part of $\xi$ and $\Im\xi$ denotes the imaginary part of $\xi$. This result concludes the decoupling part of our work. The expressions (B.36), together with the decoupling form (B.32), can be used in the expression for $\hat{U}(\tau)$ (B.3) to obtain an explicit (up to a formal solution for $\xi(\tau)$) time-evolved expression for the observables of the system.

Finally, let us once more write down the final expression of the time-evolution operator

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:app:evolution:operator} \hat{U}(\tau)&= \, {\rm e}^{-{\rm i}\,\Omega_\mathrm{c} \hat{a}^{\dagger} \hat{a}\,\tau}\,\hat{U}_{\mathrm{sq}}(\tau)\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a}\,\hat{N}_a}\,{\rm e}^{-{\rm i}\,F_{\hat{N}^2_a}\,\hat{N}^2_a}\,{\rm e}^{-{\rm i}\,F_{\hat{B}_+}\,\hat{B}_+}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_+}\,\hat{N}_a\,\hat{B}_+}\,\nonumber \\ &\quad\times {\rm e}^{-{\rm i}\,F_{\hat{B}_-}\,\hat{B}_-}\,{\rm e}^{-{\rm i}\,F_{\hat{N}_a\,\hat{B}_-}\,\hat{N}_a\,\hat{B}_-}, \nonumber \end{align} \tag{ B.37 }$$

to be complemented with the functions listed in (B.36).

In this appendix, we use the covariance matrix elements to compute the symplectic eigenvalues of the optical and mechanical subsystems.

C.1.1. The optical symplectic eigenvalue.

We wish to compute the symplectic eigenvalue of the optical subsystem. It is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu_{\rm Op}^2 = \sigma_{11}^2 - |\sigma_{13}|^2 \, . \nonumber \end{align} \tag{ C.8 }$$

From the covariance matrix elements in (C.7) we find

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{11}^2 &= 1 + 4\, |\mu_{\mathrm{c}}|^2 \left( 1 - {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} \,{\rm e}^{- |K_{\hat{N}_a}|^2} \right)+ 4 \, |\mu_{\mathrm{c}}|^4 \left( 1 - {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} \,{\rm e}^{- |K_{\hat{N}_a}|^2} \right) ^2 \, ,\nonumber \\ |\sigma_{31}|^2 &= 4\, |\mu_{\mathrm{c}}|^4 \, |E_{\hat{B}_+ \hat{B}_-}|^4 \left( {\rm e}^{{\rm i} \theta} \, {\rm e}^{|\mu_{\mathrm{c}}|^2 ( {\rm e}^{2{\rm i} \theta} - 1)} {\rm e}^{-|K_{\hat{N}_a}|^2} \, - {\rm e}^{2|\mu_{\mathrm{c}}|^2 ( {\rm e}^{{\rm i} \theta} - 1)} \right) \, \nonumber \\ &\quad\quad\quad\quad\quad\quad\quad\quad\times \left( {\rm e}^{-{\rm i} \theta} \, {\rm e}^{|\mu_{\mathrm{c}}|^2 ( {\rm e}^{- 2{\rm i} \theta} - 1)} {\rm e}^{-|K_{\hat{N}_a}|^2} \, - {\rm e}^{2|\mu_{\mathrm{c}}|^2 ( {\rm e}^{-{\rm i} \theta} - 1)} \right)\nonumber \\ &= 4\, |\mu_{\mathrm{c}}|^4 \, |E_{\hat{B}_+ \hat{B}_-}|^4 \biggl( {\rm e}^{-2|K_{\hat{N}_a}|^2}{\rm e}^{|\mu_c|^2 \left( {\rm e}^{2{\rm i} \theta} + {\rm e}^{- 2 {\rm i} \theta} -2\right)} + {\rm e}^{2|\mu_{\mathrm{c}}|^2 ( {\rm e}^{{\rm i} \theta} + {\rm e}^{-{\rm i} \theta} - 2)} \nonumber \\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad- 2 \, \Re \left[ {\rm e}^{{\rm i} \theta} \, {\rm e}^{|\mu_{\mathrm{c}}|^2 ( {\rm e}^{2{\rm i} \theta} - 1)} {\rm e}^{-|K_{\hat{N}_a}|^2} {\rm e}^{2 |\mu_c|^2 ({\rm e}^{-{\rm i} \theta} - 1)} \right] \biggr) . \nonumber \end{align} \tag{ C.9 }$$

By using the following trigonometric identities,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm e}^{2{\rm i} \theta} + {\rm e}^{-2 {\rm i} \theta } - 2 &= 2 \cos{2\theta} - 2= - 4 \sin^2{\theta} \, , \nonumber \\ {\rm e}^{{\rm i} \theta} + {\rm e}^{-{\rm i} \theta } - 2 &= 2( \cos\theta - 1) = - 4 \sin^2{\theta/2} \, , \nonumber \end{align} \tag{ C.10 }$$

and from the fact that $|E_{\hat{B}_+ \hat{B}_-}|^2 = {\rm e}^{- |K_{\hat{N}_a}|^2 }$, we obtain

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{11}^2 &= 1 + 4\, |\mu_{\mathrm{c}}|^2 \left( 1 - {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} \,{\rm e}^{- |K_{\hat{N}_a}|^2} \right)+ 4 \, |\mu_{\mathrm{c}}|^4 \left( 1 - {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} \,{\rm e}^{- |K_{\hat{N}_a}|^2} \right) ^2 \, ,\nonumber \\ |\sigma_{31}|^2 &= 4\, |\mu_{\mathrm{c}}|^4 \, {\rm e}^{-2 |K_{\hat{N}_a}|^2 } \biggl( {\rm e}^{-2|K_{\hat{N}_a}|^2}{\rm e}^{- 4|\mu_c|^2 \sin^2{\theta}} + {\rm e}^{-8 |\mu_{\mathrm{c}}|^2 \sin^2\theta/2} \nonumber \\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad- 2 \, \Re \left[ {\rm e}^{{\rm i} \theta} \, {\rm e}^{|\mu_{\mathrm{c}}|^2 ( {\rm e}^{2{\rm i} \theta} - 1)} {\rm e}^{-|K_{\hat{N}_a}|^2} {\rm e}^{2 |\mu_c|^2 ({\rm e}^{-{\rm i} \theta} - 1)} \right] \biggr) \, . \nonumber \end{align} \tag{ C.11 }$$

Putting them together, we find

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \nu_{\rm Op}^2 &= 1 +4 \, |\mu_{\mathrm{c}}|^2 \left( 1 - {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} \,{\rm e}^{- |K_{\hat{N}_a}|^2} \right) \nonumber \\ &+ 4 \, |\mu_c|^4 \biggl( 1 - 2 \, {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} {\rm e}^{- |K_{\hat{N}_a}|^2 } + {\rm e}^{-8|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} {\rm e}^{- 2|K_{\hat{N}_a}|^2 } \nonumber \\ &\quad\quad\quad\quad\quad- {\rm e}^{-4 |K_{\hat{N}_a}|^2 } {\rm e}^{- 4 |\mu_c|^2 \sin^2\theta} - {\rm e}^{- 2 |K_{\hat{N}_a}|^2 } {\rm e}^{- 8 |\mu_{{\rm c}}|^2 \sin^2\theta/2} \nonumber \\ &\quad\quad\quad\quad\quad+ 2 \, {\rm e}^{- 3|K_{\hat{N}_a}|^2} \, \Re \left[ {\rm e}^{{\rm i} \theta} \, {\rm e}^{|\mu_{\mathrm{c}}|^2 ( {\rm e}^{2{\rm i} \theta} - 1)} {\rm e}^{2 |\mu_c|^2 ({\rm e}^{-{\rm i} \theta} - 1)} \right] \biggr) \, , \nonumber \end{align} \tag{ C.12 }$$

which can be simplified into

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \nu_{\rm Op}^2 &= 1 +4 \, |\mu_{\mathrm{c}}|^2 \left( 1 - {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} \,{\rm e}^{- |K_{\hat{N}_a}|^2} \right) \nonumber \\ &+ 4 |\mu_c|^4 \biggl( 1 - 2 {\rm e}^{-4|\mu_{\mathrm{c}}|^2 \sin^2{\theta/2}} {\rm e}^{- |K_{\hat{N}_a}|^2 } - {\rm e}^{-4 |K_{\hat{N}_a}|^2 } {\rm e}^{- 4 |\mu_c|^2 \sin^2\theta} \nonumber \\ &\quad\quad\quad\quad\quad\quad+ 2 \, {\rm e}^{- 3|K_{\hat{N}_a}|^2} \, \Re \left[ {\rm e}^{{\rm i} \theta} \, {\rm e}^{|\mu_{\mathrm{c}}|^2 ( {\rm e}^{2{\rm i} \theta} - 1)} \, {\rm e}^{2 |\mu_c|^2 ({\rm e}^{-{\rm i} \theta} - 1)} \right] \biggr) \, . \nonumber \end{align} \tag{ C.13 }$$

Next, we compute the mechanical symplectic eigenvalue.

C.1.2. The mechanical symplectic eigenvalue.

We first recall the Bogoliubov identities, which are $|\alpha(\tau)|^2 - |\beta(\tau)|^2 = 1$ and $\alpha(\tau) \beta^*(\tau) - \alpha^*(\tau) \beta(\tau) = 0$.

The mechanical eigenvalue is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu_{{\rm Me}}^2 = \sigma_{22}^2 - |\sigma_{42}|^2 \, . \nonumber \end{align} \tag{ C.14 }$$

Given the covariance matrix elements in (C.7), we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{22}^2 &= 1 + 4 |\beta|^2 + 4 |\beta|^4+ 4 |\Delta|^2 |\mu_{{\rm c}}|^2 + 8 |\beta|^2 \, |\Delta|^2 |\mu_{{\rm c}}|^2 + 4 |\Delta|^4 |\mu_{{\rm c}}|^4 \, ,\nonumber \\ |\sigma_{42}|^2 &= 4 |\alpha|^2 |\beta|^2 + 4 \alpha^* \beta^* \Delta^2 |\mu_{{\rm c}}|^2 + 4 \alpha \beta \Delta^{*2} |\mu_{{\rm c}}|^2 + 4 |\Delta|^4 |\mu_{{\rm c}}|^4 \, . \nonumber \end{align} \tag{ C.15 }$$

This allows us to write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu_{{\rm Me}}^2 = \sigma_{22}^2 - |\sigma_{42}|^2 = 1 + 4\left[ \left( 1 + 2 |\beta|^2 \right) |\Delta|^2 - 2 \, \Re \left( \alpha \, \beta \Delta^{*2} \right) \right] |\mu_{{\rm c}}|^2 \, , \nonumber \end{align} \tag{ C.16 }$$

where we have suppressed the dependence of $\tau$ for notational clarity.

We wish to simplify this expression by examining each term in turn and using the Bogoliubov conditions. We recall that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta = \left( \alpha + \beta \right) F_{\hat{N}_a \, \hat{B}_-} - {\rm i} \, \left( \alpha - \beta \right) F_{\hat{N}_a \,\hat{B}_+} \, . \nonumber \end{align} \tag{ C.17 }$$

We can now use the Bogoliubov identities to show that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |\Delta|^2 = \left( 1 + 2|\beta|^2 \right) \, |K_{\hat{N}_a}|^2 - ( \alpha \beta^* + \alpha^* \beta) \left( F_{\hat{N}_a \, \hat{B}_+}^2 - F_{\hat{N}_a \, \hat{B}_-}^2 \right) \, , \nonumber \end{align} \tag{ C.18 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta^{*2} &= \, ( \alpha^* + \beta^*)^2 F_{\hat{N}_a \, \hat{B}_-}^2 - ( \alpha^* - \beta^*)^2 F_{\hat{N}_a \, \hat{B}_+}^2 + 2 \, i \, ( \alpha^{*2} - \beta^{*2}) \, F_{\hat{N}_a \, \hat{B}_+} \, F_{\hat{N}_a \, \hat{B}_-} \nonumber \\ &= \, ( \alpha^{*2} + \beta^{*2} + 2 \alpha^* \beta^*) F_{\hat{N}_a \, \hat{B}_-}^2 - ( \alpha^{*2} + \beta^{*2} - 2 \, \alpha^* \beta^*) \, F_{\hat{N}_a \, \hat{B}_+}^2 \nonumber \\ &\quad+ 2 \, i \, ( \alpha^{*2} - \beta^{*2}) F_{\hat{N}_a \, \hat{B}_+} F_{\hat{N}_a \, \hat{B}_-} \nonumber \\ &= \, \left( \alpha^{*2} + \beta^{*2} \right) \left( F_{\hat{N}_a \, \hat{B}_-}^2 - F_{\hat{N}_a \, \hat{B}_+}^2 \right) + 2\, \alpha^* \beta^* |K_{\hat{N}_a}|^2 \nonumber \\ &\quad+ 2 {\rm i} \left( \alpha^{*2} - \beta^{*2} \right) F_{\hat{N}_a \, \hat{B}_+} \, F_{\hat{N}_a \, \hat{B}_-} \, , \nonumber \end{align} \tag{ C.19 }$$

where we recall that $|K_{\hat{N}_a}|^2 = F_{\hat{N}_a \, \hat{B}_+}^2 + F_{\hat{N}_a \, \hat{B}_-}^2$. Furthermore, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha \beta \, \Delta^{*2} &= \, \left( |\alpha|^2 \alpha^* \beta + \alpha \beta^* |\beta|^2 \right) \left( F_{\hat{N}_a \, \hat{B}_-}^2 - F_{\hat{N}_a \, \hat{B}_+} ^2\right) + 2\, |\alpha|^2 |\beta|^2 |K_{\hat{N}_a}|^2 \nonumber \\ &+ 2 {\rm i} \left( |\alpha|^2 \alpha^* \beta - \alpha \beta^* |\beta|^2 \right) F_{\hat{N}_a \,\hat{B}_+} F_{\hat{N}_a \, \hat{B}_-} \nonumber \\ &= \, \alpha^* \beta \left( F_{\hat{N}_a \, \hat{B}_-}^2 - F_{\hat{N}_a \, \hat{B}_+}^2 \right) + |\beta|^2 \left( \alpha^* \beta + \alpha \beta^* \right) \left( F_{\hat{N}_a \, \hat{B}_-}^2 - F_{\hat{N}_a \, \hat{B}_+}^2 \right) \nonumber \\ &+ 2 \left( 1 + |\beta|^2 \right) |\beta|^2 |K_{\hat{N}_a}|^2 + 2 {\rm i} \alpha^* \beta F_{\hat{N}_a \, \hat{B}_+} F_{\hat{N}_a \, \hat{B}_-} \nonumber \\ &+ 2 {\rm i} |\beta|^2 \left( \alpha^* \beta - \alpha \beta^* \right) F_{\hat{N}_a \, \hat{B}_+} F_{\hat{N}_a \, \hat{B}_-} \, , \nonumber \end{align} \tag{ C.20 }$$

where we used $|\alpha|^2 = |\beta|^2 + 1$ everywhere. We now note that the last term disappears because $\alpha^* \beta - \alpha \beta^* = 0$. We also note that $\alpha^* \beta = \frac{1}{2} ( \alpha ^* \beta + \alpha \beta^*)$ is real, which follows from the Bogoliubov condition $\alpha^* \beta = \alpha \beta^*$. When we take the real part, the second-to-last term disappears as well because it has an additional i, meaning that we are left with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Re \left( \alpha \beta \Delta^{*2} \right) &= \frac{1}{2} \left( \alpha^* \beta + \alpha \beta^* \right) \left( F_{\hat{N}_a \, \hat{B}_-}^2 - F_{\hat{N}_a \, \hat{B}_+}^2 \right) + |\beta|^2 \left( \alpha^* \beta + \alpha \beta^* \right) \left( F_{\hat{N}_a \, \hat{B}_-}^2 - F_{\hat{N}_a \, \hat{B}_+}^2 \right) \nonumber \\ &+ 2\left( 1 + |\beta|^2 \right) |\beta|^2 |K_{\hat{N}_a}|^2 \, . \nonumber \end{align} \tag{ C.21 }$$

We turn again to the symplectic eigenvalue, which we can now simplify as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \nu_{{\rm Me}}^2 &= 1 + 4\left[ \left( 1 + 2 |\beta|^2 \right) |\Delta|^2 - 2 \, \Re \left( \alpha \, \beta \Delta^{*2} \right) \right] |\mu_{{\rm c}}|^2 \nonumber \\ &= 1 + 4 \biggl[ ( 1 + 2 |\beta|^2)^2 |K_{\hat{N}_a}|^2 - ( 1 + 2 |\beta|^2) \left( \alpha \beta^* + \alpha^* \beta\right) \left( F_{\hat{N}_a \, \hat{B}_+}^2 - F_{\hat{N}_a \, \hat{B}_-}^2 \right) \nonumber \\ &\quad\quad\quad- \left( \alpha^* \beta + \alpha \beta^* \right) \left( F_{\hat{N}_a \, \hat{B}_-}^2 - F_{\hat{N}_a \, \hat{B}_+}^2 \right) \nonumber \\ &\quad\quad\quad- 2 |\beta|^2 (\alpha^* \beta + \alpha \beta^*) \left( F_{\hat{N}_a \, \hat{B}_-} ^2 - F_{\hat{N}_a \, \hat{B}_+}^2 \right) \nonumber \\ &\quad\quad\quad- 4 (1 + |\beta|^2) |\beta|^2 |K_{\hat{N}_a}|^2 \biggr] |\mu_{{\rm c}}|^2 \nonumber \\ &= 1 + 4 \left[ \left( 1 + 4 |\beta|^2 + 4|\beta|^4 \right) |K_{\hat{N}_a}|^2 - 4 \left( |\beta|^2 + |\beta|^4 \right) |K_{\hat{N}_a}|^2 \right] |\mu_{{\rm c}}|^2 \nonumber \\ &= 1 + 4 |K_{\hat{N}_a}|^2 |\mu_{{\rm c}}|^2 \, . \nonumber \end{align} \tag{ C.22 }$$

Our goal is to obtain approximate solutions to the differential equations (B.18) and (B.19). We will do so by following the derivation in [53] with some modifications.

The general form of Mathieu's equation is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:Mathieu} \frac{{\rm d}^2 y}{{\rm d}x^2} + \left[ a - 2 \, q \, \cos(2 \, x) \right] \, y= 0 \, . \nonumber \end{align} \tag{ E.1 }$$

We will use the general notation in this appendix and then compare with the system in the main text.

We begin by defining a slow time scale X  =  qx. We then assume that the solutions y  depend on both scales, such that $y(x, X)$. This means that we can treat x and X as independent variables and the absolute derivative ${\rm d}/{\rm d}x$ in (E.1) can be split in two:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\rm d}{{\rm d}x} = \partial_x + q \, \partial_X \, . \nonumber \end{align} \tag{ E.2 }$$

Mathieu's equation (E.1) therefore becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \partial_x + q \, \partial_X \right)^2 y(x, X) + ( a - 2 q \cos(2 x)) \, y(x, X) = 0 \, . \nonumber \end{align} \tag{ E.3 }$$

We then expand the solution $y(x, X)$ for small q as $y(x,X) = y_0(x,X) + q\, y_1 (x,X) + \mathcal{O}(q^2)$ and insert this into the differential equation above. Our goal is to obtain a solution for y ₀ which incorporates a number of restrictions from the differential equation for y ₁.

To zeroth order, we recover the regular harmonic differential equation for y ₀, which is the limiting case as $q \rightarrow 0$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \partial^2_x y_0 + a \, y_0 = 0 \, , \nonumber \end{align} \tag{ E.4 }$$

where we know that the solutions are sinusoidal, while the coefficients must depend on X. We choose the following trial solution:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y_0(x,X) = A(X) \, {\rm e}^{{\rm i} \sqrt{a} \, x} +A^*(X) \, {\rm e}^{-{\rm i} \sqrt{a} \, x} \, . \nonumber \end{align} \tag{ E.5 }$$

Our goal is now to find explicit solutions to the complex function $A(X)$. We continue with the equation for y ₁. We discard all terms of order q² to find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q \, \partial^2_x y_1 + 2\, q \, \partial_x \partial _X y_0 + a \, q \, y_1 -2 \, q \cos(2 x) y_0 = 0 \, . \nonumber \end{align} \tag{ E.6 }$$

We divide by q and insert our solution for y ₀ to find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\partial^2_x y_1 + a \, y_1 + 2 \, {\rm i} \, \sqrt{a} \, \left(\frac{\partial A(X)}{\partial X} \, {\rm e}^{{\rm i} \sqrt{a} x} - \frac{\partial A^*(X)}{\partial X} \, {\rm e}^{-{\rm i} \sqrt{a} x} \right) \nonumber \\ &\quad- 2 \, \cos(2 x) \left( A(X) \, {\rm e}^{{\rm i} \sqrt{a} x} + A^*(X) \, {\rm e}^{-{\rm i} \sqrt{a} x} \right) = 0 \, . \nonumber \end{align} \tag{ E.7 }$$

At this point, we specialise to a  =  1, which corresponds to setting $\Omega_0 = 2$ in the main text. We combine the exponentials to find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\partial^2_x \, y_1 + a\, y_1 + \left( 2 {\rm i} \frac{\partial A(X)}{\partial X} - \, A^*(X) \right) \, {\rm e}^{{\rm i} x} + \left( 2 {\rm i} \frac{\partial A^*(X)}{\partial X} + \, A(X) \right) \, {\rm e}^{-{\rm i} x} \nonumber \\ &\quad - A(X) \, {\rm e}^{3 {\rm i} x} - A^*(X) \, {\rm e}^{- 3 {\rm i} x} = 0 \, . \nonumber \end{align} \tag{ E.8 }$$

In order for the solution to be stable, we require that secular terms such as resonant terms ${\rm e}^{{\rm i} x}$ vanish. If these do not vanish, the solution will grow exponentially [53]. We also neglect terms that oscillate much faster, such as ${\rm e}^{3 {\rm i} x}$. This leaves us with the condition that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:resonance:condition} \left( 2 {\rm i} \frac{\partial A^*(X)}{\partial X} + \, A(X) \right) = 0 \, , \nonumber \end{align} \tag{ E.9 }$$

which can be differentiated again and solved with the trial solution $A(X) = (c_1- {\rm i} \, c_2)$${\rm e}^{X/2} + (c_3- {\rm i}\, c_4) \, {\rm e}^{- X/2}$ for the parameters $c_1,c_2,c_3$ and c₄. From the requirement in (E.9), it is now possible to fix two of the coefficients in (E.10). We differentiate $A(X)$ and use (E.9) to find that the conditions $c_1= c_2$ and $c_3 = -c_4$ must always be fulfilled.

We then recall that X  =  qx and after combining some exponentials, we obtain the full trial solution for the zeroth order term y ₀:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:y0:almost:final:solution} y_0(x) &= A(qx) \, {\rm e}^{{\rm i} \, x} +A^*(qx) \, {\rm e}^{-{\rm i} x} \, \nonumber \\ &= 2 \, \left( g_0 c_1 \, {\rm e}^{qx /2} + c_3 \, {\rm e}^{- qx/2} \right) \, \cos(x) + 2 \, \left( g_0c_1 \, {\rm e}^{qx/2} - c_3 \, {\rm e}^{- qx/2} \right) \sin(x) \, . \nonumber \end{align} \tag{ E.10 }$$

We now proceed to compare this solution with the parameters and initial conditions given for P₁₁ in (B.18) and $I_{P_{22}}$ in (B.19) in the main text, which are both solved by the Mathieu equation.

First, we note that $q = - 2\, \tilde{d}_2$ and that $x = \tau$. Then we consider the boundary conditions for P₁₁, which are P₁₁(0)  =  1 and $\dot{P}_{11} (0) = 0$. From these conditions, we find that $g_0 c_1= c_3 = 1/4$, and the the approximate solution to P₁₁ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:eq:P11:approx} P_{11}(\tau) =\cos (\tau) \, \cosh (\tilde{d}_2 \, \tau) - \sin (\tau) \, \sinh ( \tilde{d}_2 \, \tau) \, . \nonumber \end{align} \tag{ E.11 }$$

The equation for $I_{P_{22}}$ has the opposite initial conditions $I_{P_{22}}(0) = 0$ and $\dot{I}_{P_{22}} = 1$. For this case, we find that $g_0c_1= -c_3 = 1/(4 (1 - \tilde{d}_2))$. The full solution to $I_{P_{22}}$ is therefore

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:eq:IP22:approx} I_{P_{22}}(\tau) = - \frac{1}{1- \tilde{d}_2}\left(\cos (\tau) \, \sinh (\tilde{d}_2 \, \tau) - \sin (\tau) \, \cosh (\tilde{d}_2 \, \tau) \right) \, , \nonumber \end{align} \tag{ E.12 }$$

and thus

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{22} = \cos (\tau) \cosh (\tilde{d}_2 \, \tau)-\frac{\tilde{d}_2 + 1}{\tilde{d}_2-1} \sin (\tau) \sinh (\tilde{d}_2 \, \tau) \, . \nonumber \end{align} \tag{ E.13 }$$

Both solutions reduce to the correct solutions for the zero-squeezing case as $\tilde{d}_2 \rightarrow 0$.

The validity of the pertubative approach can be determined as follows. By inserting the trial solution for $P_{11}(\tau)$ into the Mathieu equation, we are left with terms that are multiplied by $\tilde{d}_2$. The leading term is in fact $\tilde{d}_2 \, \cosh(\tilde{d}_2 \, \tau)$. These terms must be approximately zero to solve the Mathieu equation, which means that we require $\tilde{d}_2 \, \cosh(\tilde{d}_2 \, \tau) \ll1$. This means that while $\tilde{d}_2 \ll1$, we can allow $\tilde{d}_2 \, \tau \sim 1$, which means that $\tau$ can be large provided that $\tilde{d}_2$ is sufficiently small.

From the expression for $\xi(\tau)$ in (B.29) we then find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi(\tau) &=\cos (\tau) \, \cosh (\tilde{d}_2 \, \tau) - \sin (\tau) \, \sinh (\tilde{d}_2 \, \tau) \nonumber \\ &- \frac{{\rm i}}{1- \tilde{d}_2} \, \left( \sin (\tau) \cosh (\tilde{d}_2 \, \tau) -\cos (\tau) \sinh (\tilde{d}_2 \, \tau)\right) \, . \nonumber \end{align} \tag{ E.14 }$$

For very small $\tilde{d}_2 \ll 1$, which was the condition for deriving the approximate solutions in the first place, we can approximate the fraction as unity and we find the compact expression

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:eq:RWA} \xi(\tau) = \, {\rm e}^{-{\rm i} \, \tau} \, \cosh (\tilde{d}_2 \, \tau) + {\rm i} \, {\rm e}^{{\rm i} \, \tau} \, \sinh (\tilde{d}_2 \, \tau) \, . \nonumber \end{align} \tag{ E.15 }$$

To better understand what this approximation entails physically, we compare it with the RWA, which has a well-known physical interpretation.

There is another solution which more explicitly demonstrates how the resonance conditions helps constrain the solution. We write the solution to the differential equation for P₁₁ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{11}(\tau) = Q_{c}(\tau) \cos(\tau + \pi/4) + Q_{s}(\tau) \sin(\tau + \pi/4). \nonumber \end{align} \tag{ E.16 }$$

Then, the differential equation $\ddot P_{11} + (1+f(\tau)) P_{11} = 0$ is solved approximately, considering only terms of first order in $\tilde{d}_2$ and neglecting terms rotating with frequency 3 off resonantly, by the following set of differential equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \dot Q_c(\tau) = \tilde{d}_2 Q_c(\tau) \quad\mathrm{and}\quad\dot Q_s(\tau) = - \tilde{d}_2 Q_s(\tau)\, . \nonumber \end{align} \tag{ E.17 }$$

These can be solved as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q_c(\tau) = {\rm e}^{\tilde{d}_2 \tau} Q_c(0) \quad\mathrm{and}\quad Q_s(\tau) = {\rm e}^{- \tilde{d}_2 \tau} Q_s(0). \nonumber \end{align} \tag{ E.18 }$$

As the initial conditions for P₁₁ are P₁₁(0)  =  1 and $\dot P_{11}(0)=0$, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q_{c}(0) + Q_{s}(0) = \sqrt{2} \, , \quad{\rm and}\quad \left(1 - \tilde{d}_2\right) \left( Q_{c}(0) - Q_{s}(0) \right) = 0 \, , \nonumber \end{align} \tag{ E.19 }$$

which implies $Q_{c}(0) = Q_{s}(0) = 1/\sqrt{2}$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{11}(\tau) = \frac{1}{\sqrt{2}}\left( {\rm e}^{\tilde{d}_2\tau} \cos(\tau + \pi/4) + {\rm e}^{-\tilde{d}_2\tau} \sin(\tau + \pi/4)\right). \nonumber \end{align} \tag{ E.20 }$$

The same steps as above can be applied to find an approximate solution for $I_{P_{22}}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_{P_{22}}(\tau) = \bar Q_{c}(\tau) \cos(\tau + \pi/4) + \bar Q_{s}(\tau) \sin(\tau + \pi/4), \nonumber \end{align} \tag{ E.21 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar Q_c(\tau) = {\rm e}^{\tilde{d}_2\tau} \bar Q_c(0) \quad\mathrm{and}\quad \bar Q_s(\tau) = {\rm e}^{-\tilde{d}_2\tau} \bar Q_s(0)\, \nonumber \end{align} \tag{ E.22 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar Q_{c}(0) + \bar Q_{s}(0) = 0 \, , \quad{\rm and}\quad - \left(1 - \tilde d_2\right) \left( \bar Q_{c}(0) - \bar Q_{s}(0) \right) = \sqrt{2} \, , \nonumber \end{align} \tag{ E.23 }$$

which leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar Q_c(0)= -\frac{1}{\sqrt{2}\left(1 - \tilde d_2\right)} \, ,\quad \mathrm{and} \quad \bar Q_s(0)= \frac{1}{\sqrt{2}\left(1 - \tilde d_2\right)} \, , \nonumber \end{align} \tag{ E.24 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_{P_{22}}(\tau) = -\frac{1}{\sqrt{2}\left(1 - \tilde d_2\right)}\left( {\rm e}^{\tilde{d}_2\tau} \cos\left(\tau + \frac{\pi}{4}\right) - {\rm e}^{-\tilde{d}_2\tau} \sin\left(\tau + \frac{\pi}{4}\right)\right). \nonumber \end{align} \tag{ E.25 }$$

We find that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi(\tau) = \frac{1}{\sqrt{2}}\left(\left( 1 + {\rm i}\frac{1}{\left(1 - \tilde d_2\right)}\right) {\rm e}^{\tilde{d}_2\tau} \cos\left(\tau + \frac{\pi}{4}\right) + \left( 1 - {\rm i}\frac{1}{\left(1 - \tilde d_2\right)}\right) {\rm e}^{-\tilde{d}_2\tau} \sin\left(\tau + \frac{\pi}{4}\right)\right) \, \, , \nonumber \end{align} \tag{ E.26 }$$

which exactly coincides with (E.15).

Here we compare the approximate resonance solution for P₁₁ in (E.11) and $I_{P_{22}}$ in (E.12) with the RWA, which is obtained as an approximation to the Hamiltonian in (1) when $\tau \gg 1$. In the main text, we separated the Hamiltonian (1) into a free term and a squeezing term (B.4), which for our specific choice of $\tilde{\mathcal{D}}_2(\tau) = \tilde{d}_2 \, \cos(\Omega_0 \, \tau)$ becomes

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \hat{{H}}_{{\rm sq}} = \hat{b}^{\dagger} \hat{b} + \tilde{d}_2 \cos(\Omega_0 \tau) \left( \hat{b}^{\dagger} + \hat{b} \right)^2 \, . \nonumber \end{align} \tag{ E.27 }$$

We now define the free evolution Hamiltonian $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{H}_0 = \hat{b}^{\dagger} \hat{b}$ and the squeezing term $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{H}'_{{\rm sq}} = \tilde{d}_2 \cos(\Omega_0 \, \tau) ( \hat{b}^{{\dagger}} + \hat{b}){}^2$ as separate operators. We transform into a frame rotating with $\newcommand{\re}{{\rm Re}} \newcommand{\e}{{\rm e}} \renewcommand{\dag}{\dagger}\exp[-{\rm i} \, \tau \, \hat{b}^{\dagger} \hat{b} ]$, which means that the squeezing term transforms into

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:rotating:frame:more:exponentials} {\rm e}^{{\rm i} \hat{H}_0 \, \tau} \hat{H}_{{\rm sq}}' {\rm e}^{- {\rm i} \hat{H}_0 , \tau} = \, \tilde{d}_2 \, \cos(\Omega_0 \, \tau) \left( {\rm e}^{2 \, {\rm i} \, \tau} \, b^{2\dagger} + {\rm e}^{- 2 \, {\rm i} \, \tau} \, \hat{b}^2 + 2 \, \hat{b}^{\dagger} \hat{b} + 1\right) \, . \nonumber \end{align} \tag{ E.28 }$$

For this specific case, the system becomes resonant when $\Omega_0 = \omega_0/\omega_{\mathrm{m}} = 2$. We can see this by expanding the cosine in terms of exponentials to obtain

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\rm e}^{{\rm i} \hat{H}_0 \, \tau} \hat{H}_{{\rm sq}}' {\rm e}^{- {\rm i} \hat{H}_0 , \tau} &= \frac{1}{2} \tilde{d}_2 \left[ \left( {\rm e}^{{\rm i} \,(2 + \Omega_0) \, \tau} + {\rm e}^{{\rm i} \, (2 - \Omega_0)\, \tau} \right) \hat{b}^{2\dagger} + \left( {\rm e}^{-{\rm i} \, (2 + \Omega_0) \, \tau} + {\rm e}^{-{\rm i} \, ( 2 - \Omega_0)\, \tau} \right) \hat{b}^2 \right] \, \nonumber \\ &+ \tilde{d}_2 \, \cos(\Omega_0 \, \tau) \left( 2 \, \hat{b}^{\dagger} \hat{b} + 1 \right) \, . \nonumber \end{align} \tag{ E.29 }$$

When $\Omega_0 = 2$, two of the time-dependent terms will cancel. We then perform the RWA, i.e. we neglect all remaining time-dependent terms. This approximation is only valid for $\tau \gg1$. In the interaction frame, we find

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{app:rotating:frame:approx} \hat H_{{\rm sq}, I} = {\rm e}^{{\rm i} \hat H_0\, \tau} \hat{H}_{{\rm sq}}' {\rm e}^{- {\rm i} \hat H_0 \tau} \approx \frac{1}{2} \tilde{d}_2 \, ( \hat{b}^{{\dagger} 2} + \hat{b}^2) \, . \nonumber \end{align} \tag{ E.30 }$$

In the symplectic basis $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat {\vec r} = (\hat{b}, \hat{b} ^{\dagger}){}^{{\rm T}}$ the corresponding symplectic operator, given by $\boldsymbol{S}_{{\rm sq}} = {\rm e}^{\boldsymbol{\Omega} \, \boldsymbol{H}_{{\rm sq}}}$, where $\boldsymbol{\Omega} = {\rm i} \, \mathrm{diag}(-1, 1)$ and $\boldsymbol{H}_{{\rm sq}}$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{H}_{{\rm sq}} = \tilde{d}_2 \, \left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ 1 & 0 \end{array}\right) \, . \nonumber \end{align} \tag{ E.31 }$$

The symplectic representation of the squeezing operator (B.7) in the lab frame reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app:lab:frame:evolution} \boldsymbol{S}_{\mathrm{sq}}(\tau) =\boldsymbol{S}_0(\tau) \,\left(\begin{array}{@{}cc@{}} \cosh (\tilde{d}_2 \,\tau) & -{\rm i}\,\sinh(\tilde{d}_2 \,\tau) \nonumber \\ {\rm i}\,\sinh(\tilde{d}_2 \,\tau) & \cosh (\tilde{d}_2 \,\tau) \end{array}\right), \nonumber \end{align} \tag{ E.32 }$$

where $\boldsymbol{S}_0 = {\rm e}^{-{\rm i} \, \tau}$ encodes the evolution from the Hamiltonian $\hat{H}_0$. We therefore find the Bogoliubov coefficients

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(\tau)&= \, {\rm e}^{-{\rm i} \, \tau}\,\cosh (\tilde{d}_2 \,\tau) \, ,\nonumber \\ \beta(\tau)&=-{\rm i}\, {\rm e}^{-{\rm i} \, \tau}\,\sinh (\tilde{d}_2 \,\tau) \, , \nonumber \end{align} \tag{ E.33 }$$

which evidently satisfy the Bogoliubov conditions, and obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi = \alpha(\tau) + \beta^*(\tau) = {\rm e}^{-{\rm i} \, \tau} \, \cosh (\tilde{d}_2 \, \tau) + {\rm i} \, {\rm e}^{{\rm i} \, \tau} \, \sinh (\tilde{d}_2 \, \tau) \, . \nonumber \end{align} \tag{ E.34 }$$

This expression exactly matches the one we derived as a perturbative solution to the Mathieu equations in (E.15). However, the requirement for the validity of the RWA is that $\tau \gg1$, while the approximate solutions are still valid for small $\tau$. We conclude that the approximate solutions only coincide with the RWA for large $\tau$, while this interpretation cannot be used when $\tau \sim1$.