Time evolution of coupled multimode and multiresonator optomechanical systems

TL;DR

This paper derives analytical solutions for the time evolution of complex multimode optomechanical systems, enabling detailed analysis of quantum coherence, mixedness, and the effects of nonlinear interactions versus linear approximations.

Contribution

It provides the first analytical expressions for the full time evolution of coupled multimode optomechanical systems with arbitrary interactions, including nonlinear effects.

Findings

Analytical expressions for system evolution and expectation values.

Identification of regimes with significant nonlinear effects.

Comparison showing differences between nonlinear and linearized models.

Abstract

We study the time evolution of bosonic systems where multiple driven bosonic modes of light interact with multiple mechanical resonators through arbitrary, time-dependent, optomechanical-like interactions. We find the analytical expression for the full time evolution of the system and compute the expectation value of relevant quantities of interest. Among the most interesting ones, we are able to compute the first-order quantum bipartite coherence between pairs of subsystems, and the analytical expression for the mixedness induced by the nonlinear interaction in the reduced state of the mechanical oscillators. Finally, we also compare our results with a linearised version of the system, and we find a regime where there are qualitative and quantitative differences in the behavior of some measurable quantities. Our results can therefore be used to describe the full time-evolution of the system, to characterise its nonlinear character and explore the validity of the linearisation approach.

Equations

$${\hat{H}}_{\text{full}}=$$

$${\hat{B}}_p^{(+)}=$$

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

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

$${\hat{H}}_p=$$

$$\hat{H}:=$$

$$\hat{U}(t)=$$

$$F_b=$$

$$F_n=$$

$$-2{\int }_0^{t}d{t}^{\prime }\left[g_n^{(+)}\sin ({\omega }_{\mathrm{m}}t)+g_n^{(-)}\cos ({\omega }_{\mathrm{m}}{t}^{\prime })\right]{\int }_0^{t^{\prime }}d{t}^{\prime \prime }\left[{\lambda }^{(+)}\cos ({\omega }_{\mathrm{m}}{t}^{\prime \prime })-{\lambda }^{(-)}\text{s}(t^{\prime \prime })\right]$$

$$F_{nm}=$$

$$F_{+}=$$

$$F_{-}=$$

$$F_n^{(+)}=$$

$$F_n^{(-)}=$$

$${\hat{U}}_{\text{full}}(t)=$$

$$\times e^{-i{\sum }_p{F}_{-}^{(p)}{\hat{B}}^{(p,-)}}{e}^{-i{\sum }_{np}{F}_n^{(p,-)}{\hat{N}}_n{\hat{B}}^{(p,-)}},$$

$$F_{\text{m}}^{(p)}$$

$${\tilde{F}}_{\text{c},n}^{(p)}$$

$$-2{\int }_0^{t}d{t}^{\prime }\left[g_{np}^{(+)}\sin ({\omega }_{\mathrm{m},p}{t}^{\prime })+g_{np}^{(-)}\cos ({\omega }_{\mathrm{m},p}{t}^{\prime })\right]{\int }_0^{t^{\prime }}d{t}^{\prime \prime }\left[{\lambda }^{(+)}\cos ({\omega }_{\mathrm{m},p}{t}^{\prime \prime })-{\lambda }^{(-)}\sin ({\omega }_{\mathrm{m},p}{t}^{\prime \prime })\right]$$

$$F_{nm}^{(p)}$$

$$F_{+}^{(p)}$$

$$F_{-}^{(p)}$$

$$F_n^{(p,+)}$$

$$F_n^{(p,-)}$$

$$\hat{H}=$$

$$\hat{U}(t)=$$

$${\hat{a}}_k(t)=$$

$${\hat{b}}_p(t)=$$

$${\hat{\rho }}_{\text{m}}(t)=\sum _{\stackrel{\{n_k\}}{k\in \mathcal{I}}}{p}_{\{n_k\}}{\hat{D}}_{\{n_k\}}{\hat{\rho }}_{\text{m}}(0){\hat{D}}_{\{n_k\}}^{†},$$

$$p_{\{n_k\}}:=$$

$${\hat{D}}_{\{n_k\}}:=$$

$${\hat{a}}_k^{†}(t){\hat{a}}_k(t)=$$

$${\hat{b}}_k^{†}(t){\hat{b}}_k(t)=$$

$$+|F^{(k)}{|}^2+\sum _{nm}{F}_n^{(k)}{F}_m^{(k)\ast }{\hat{N}}_n{\hat{N}}_m.$$

$$g_{mn}^{(1)}(t):=\frac{|⟨{\hat{d}}_m^{†}(t){\hat{d}}_n(t)⟩|}{\sqrt{⟨{\hat{d}}_m^{†}(t){\hat{d}}_m(t)⟩⟨{\hat{d}}_n^{†}(t){\hat{d}}_n(t)⟩}}.$$

$$S_N=$$

$$=$$

$$\hat{\rho }(0)=$$

Full text

I Introduction

Physical systems with many constituents are typically difficult to investigate in depth. Arguably, the most imposing challenge faced when studying such systems is the ability to obtain analytical insight into their dynamics, starting from fundamental equations of the theory that describes them. Classical and quantum many-body physics, statistical mechanics and thermodynamics have been developed to successfully tackle systems with large numbers of constituents. Approaches developed using tools and concepts from these areas of research rarely provide full analytical descriptions of the dynamics of complex systems, but are able to provide important coarse-grained information of key aspects of these systems. Regardless of the success of such approaches, a complete, analytical understanding of any system is highly desirable in order to explain existing features, and predict new ones.

In this work we study a quantum system composed of an arbitrary number of bosons, which we conveniently separate into field modes and mechanical oscillators. We assume that these two sets of bosonic modes interact through a nonlinear Hamiltonian, which can be used to model different physical implementations, such as Fabry–Pérot cavities with a moving-end mirror Favero and Karrai (2009), levitated nano-diamonds Barker and Shneider (2010); Yin et al. (2013), membrane-in-the-middle configurations Jayich et al. (2008) and optomechanical crystals Eichenfield et al. (2009); Safavi-Naeini et al. (2014).

We employ techniques developed to decouple the time-evolution operator analytically Puri (2001); Bruschi, Lee, and Fuentes (2013); Bruschi and Xuereb (2018),111An alternative attempt to study analytically the time evolution of such systems was developed in parallel Brown et al. (2013). and we provide a full analytical solution to the time evolution of the system. Our results apply to an arbitrary number of interacting systems with arbitrary time-dependent couplings, and are free from approximations. We employ the decoupling of the time-evolution operator to compute analytically the expectation value of physically interesting quantities, such as the average number of field- and resonator-excitations, and we compute the coherence that is induced between different bosonic modes. Furthermore, we are able to provide an analytical expression for the mixedness of the reduced state of the resonators, which can contribute to the understanding of the nonlinear nature of the dynamics. We also comment on the possibility of adding optical drive (i.e., for laser cooling purposes) to the main Hamiltonian, which leads to a system that cannot be treated with the techniques developed in this work

Finally, we apply our results to a set of simple setup, where one cavity mode interacts with an arbitrary number of operators through a standard optomechanical Hamiltonian. Our analytical results are able to provide us with an insight of the quantum behavior of the resonators in this context, and on the coherence introduced between the different subsystems. In particular, we find a simple expression for the mixedness of the resonators which is always nonzero if the nonlinear coupling is present. We are also able to compare our results with those obtained through a “linearised” version of the system, an approach typically taken in most approaches to optomechanics Aspelmeyer, Kippenberg, and Marquardt (2014). We find qualitative and quantitative differences, which we are able to quantify. We conclude discussing the validity and scope of the techniques.

This work is organised as follows. In section II we introduce the necessary tools and the core Hamiltonian. In section III we decouple the time-evolution operator analytically. In section IV we compute the time evolution of quantities of interest. In section V and section VI we specialise to a physically-relevant initial state and apply our techniques to a few examples.

II Tools

We start this work by introducing the necessary tools. The tools that we employ and develop here are not tied to a particular physical system. The only assumption is that the Hamiltonian of the system corresponds to the Hamiltonian presented below, which can be used to model, for example, a cavity optomechanical system Aspelmeyer, Kippenberg, and Marquardt (2014). Therefore, we emphasize that the approach and results of this work do not depend on the specific implementation chosen. For this reason, we choose to refer to the system as a cavity optomechanical system for simplicity of presentations, and without loss of generality.

Cavity optomechanics studies the interaction of light confined and matter Aspelmeyer, Kippenberg, and Marquardt (2014). A typical implementation is that of a cavity with a semitransparent mirror-wall, from which an electromagnetic beam can enter. The beam is in resonance with at least one cavity mode, therefore being “trapped”. On the other end, a mirror-membrane acts as the second wall of the cavity, fully reflecting light. The membrane can vibrate, which affects the fundamental frequency of the mode. The action of the membrane is then modelled effectively as the coupling of the position of a harmonic oscillator to one or more cavity modes. More membranes can be included, such as semitransparent membranes placed at the antinodes of the stationary cavity modes. When one such membrane is present, the system is known to be in a “membrane-in-the-middle” configuration Jayich et al. (2008); Aspelmeyer, Kippenberg, and Marquardt (2014).

II.1 Optomechanical Hamiltonians

In this work, we consider a bosonic system (which we refer to as an optomechanical system) with an arbitrary number of cavity field modes $\{\hat{a}_{n},\hat{a}^{\dagger}_{n}\}$ and an arbitrary number of mechanical modes $\{\hat{b}_{p},\hat{b}^{\dagger}_{p}\}$, where the creation and annihilation operators satisfy the canonical commutation relations $[\hat{a}_{n},\hat{a}^{\dagger}_{n^{\prime}}]=\delta_{nn^{\prime}}$ and $[\hat{b}_{p},\hat{b}^{\dagger}_{p^{\prime}}]=\delta_{pp^{\prime}}$, and all other vanish. We then assume that these modes interact through the generalised nonlinear optomechanical Hamiltonian

[TABLE]

where $\hat{H}_{0}:=\sum_{n}\hbar\,\omega_{\textrm{c},n}\hat{a}_{n}^{\dagger}\hat{a}_{n}+\sum_{p}\hbar\,\omega_{\textrm{m},p}\,\hat{b}_{p}^{\dagger}\hat{b}_{p}$ is the free Hamiltonian.

Here we have the cavity mode frequencies $\omega_{\textrm{c},n}$, the mechanical resonator frequencies $\omega_{\textrm{m},p}$, the time dependent couplings $\lambda^{(\pm)}_{p}(t)$ and $g_{np}^{(\pm)}(t)$, and the Hermitian operators

[TABLE]

The operators defined in (2) can be cast in a more conventional form by noting that they are proportional to the quadrature operators $\hat{x}_{\mathrm{m},p},\hat{p}_{\mathrm{m},p}$ of the mechanical resonators, i.e., $\hat{B}^{(+)}_{p}\propto\hat{x}_{\mathrm{m},p}$ and $\hat{B}^{(-)}_{p}\propto\hat{p}_{\mathrm{m},p}$. We retain our choice because the decomposition in terms of creation and annihilation operators is natural to this work. Figure (1) illustrates the general scheme, that is, it is not a representation of any particular implementation but pictorial description of the system only (1) only.

II.2 Tackling time evolution of bosonic systems

Given a set of $N$ bosonic modes and an arbitrary time-dependent Hamiltonian $\hat{H}(t)$, the unitary time-evolution operator reads

[TABLE]

where $\overset{\leftarrow}{\mathcal{T}}$ is the time ordering operator Bruschi, Lee, and Fuentes (2013). This expression simplifies dramatically when the Hamiltonian $\hat{H}$ is time independent, in which case one simply has $\hat{U}(t)=\exp[-\frac{i}{\hbar}\,\hat{H}\,t]$. However, we are here interested in general time evolution induced by the time-dependent Hamiltonian (1).

It is the main interest of this work to seek an expression for (3) of the form

[TABLE]

where the $\hat{G}_{n}$ are time-independent Hermitian operators and the $F_{n}(t)$ are real time-dependent functions. If the time evolution operator (3) can be cast in the form (4) we say that it has been decoupled. The solution to the formal expression (3) in terms of its decoupled form (4) can be found by employing specifically developed techniques Bruschi, Lee, and Fuentes (2013), which we outline in A. The results are presented in the next section.

III Decoupling the time evolution of the system

The aim of this work is to show that a decoupling (4) exists for the time evolution operator induced by the Hamiltonian (1), and to find an explicit expression for the functions $F_{n}(t)$. This can be then used to compute the time evolution of any quantity of interest, and of the state of the system. To achieve our goal, we will employ the techniques developed in Puri (2001); Bruschi, Lee, and Fuentes (2013) to decouple the time-evolution operator $\hat{U}$ induced the Hamiltonian $\hat{H}$ in (1). More details can be found in A.

We start our task by noting that we can rearrange conveniently the contributing terms to $\hat{H}(t)$ as $\hat{H}(t):=\sum_{n}\hbar\,\omega_{\textrm{c},n}\hat{a}_{n}^{\dagger}\hat{a}_{n}+\sum_{p}\,\hat{H}_{p}$, where we have introduced the multimode Hamiltonians $\hat{H}_{p}$ for fixed resonator $p$ as

[TABLE]

An important property of such Hamiltonians $\hat{H}_{p}$ is that they commute with each other at all times, i.e., $[\hat{H}_{p}(t),\hat{H}_{p^{\prime}}(t^{\prime})]=0$ for all $p,p^{\prime}$ and all $t,t^{\prime}$. This property is important for the next step.

III.1 Choice of core Hamiltonian

The first step is to employ the fact that $[\hat{H}_{p},\hat{H}_{p^{\prime}}]=0$ for all $p,p^{\prime}$. This stimulates us to look at an individual Hamiltonian $\hat{H}_{p}$ and we therefore drop the label $p$ for this first part of the work. This leaves us the Hamiltonian

[TABLE]

where we have re-defined the operators $\hat{N}_{n}:=\hat{a}^{\dagger}_{n}\hat{a}_{n}$, $\hat{N}_{b}:=\hat{b}^{\dagger}\hat{b}$, $\hat{B}^{(+)}:=\hat{b}{}^{\dagger}+\hat{b}$ and $\hat{B}^{(-)}:=i\,[\hat{b}{}^{\dagger}-\hat{b}]$ for notational convenience.

Decoupling of the Hamiltonian (6) will allow us to decouple the full Hamiltonian (1), as we will show later. The Hamiltonian (6) is an extension of the standard optomechanical Hamiltonian which, in its simplest form, is obtained from (6) by setting $g_{n}^{(-)}=\lambda^{(\pm)}=0$, considering only one cavity mode and ignoring the usual external field drive Aspelmeyer, Kippenberg, and Marquardt (2014).

III.2 Decoupled form of the time-evolution operator

We argued that an important achievement that would allow us to understand and control the time evolution of our quantum system is to be able to obtain a decoupled form (4) of the time-evolution operator induced by the Hamiltonian (6). We are able to show that such decoupling exists and has the formal solution

[TABLE]

and we have introduced the Hermitian operator $\hat{N}_{nm}:=\hat{a}^{\dagger}_{n}\hat{a}_{n}\hat{a}^{\dagger}_{m}\hat{a}_{m}$ for convenience.

The time dependent and real functions $F$ can be found explicitly Bruschi, Lee, and Fuentes (2013); Bruschi and Xuereb (2018), and the full calculations necessary to achieve them are located in A. We find

[TABLE]

We would like to emphasise that the existence of an analytical expression for the $F$ functions is remarkable in its own right. Notwithstanding the fact that the set of necessary operators (A) is infinite, and therefore the number of terms to be expected in the decoupled form (4) is also infinite, it is possible to obtain an explicit expression of the time evolution through (III.2) given an explicit expression for the time dependent couplings $\lambda_{\pm}$ and $g_{n}^{(\pm)}$.

This achievement should also shed light on the reasons why it has been difficult, so far, to fully characterise how such a nonlinear system evolves in time. A naive application of the Baker-Campbell-Hausdorff formula would have not been able to provide us in an efficient way with the full set (A) of necessary operators, which in turn allowed us to compute the functions (III.2). It should therefore be clear now why the techniques used and developed here can lead to a deeper understanding of the dynamics of nonlinear quantum systems.

III.3 Full multiresonator nonlinear decoupled solution

We have decoupled analytically the time-evolution operator induced by the Hamiltonian (6). This has led us to the analytical expressions (III.2). We recall that there is a Hamiltonian (6) for each mechanical mode $\hat{b}_{p}$. Therefore, we can reconstruct the decoupled time-evolution operator $\hat{U}_{\textrm{full}}$ induced by the Hamiltonian $\hat{H}_{\textrm{full}}$ in (1) by recalling that $[\hat{H}_{p},\hat{H}_{p^{\prime}}]=0$ for all $p,p^{\prime}$. Some algebra leads us to

[TABLE]

where the updated functions read

[TABLE]

As noted before, we find it remarkable that the time evolution for this system can be decoupled exactly and analytically. Our ability to find such solution is a consequence of the structure of the full operator algebra, which allows for the use of many “tricks” to our advantage Puri (2001); Bruschi, Lee, and Fuentes (2013). Therefore, we believe that this approach has shown a clear and marked advantage with respect to a more brute-force application of the Baker-Campbell-Hausdorff formula.

III.4 Time evolution of optically driven multimode and multiresonator systems

We would like to extend the results obtained above to tackle the full time evolution $\hat{U}(t)$ induced by an optomechanical Hamiltonian $\hat{H}$ that includes an external optical drive (potentially for each mode) in addition the interactions already considered in (1).222It is a common practice to employ such a drive in physical implementations, for purposes such as cooling. We therefore would like to consider the Hamiltonian

[TABLE]

Here we have introduced the operators $\hat{A}^{(+)}_{n}=\hat{a}_{n}{}^{\dagger}+\hat{a}_{n}$ and $\hat{A}^{(-)}_{n}=i\,\left[\hat{a}_{n}{}^{\dagger}-\hat{a}_{n}\right]$, and the time-dependent couplings $\xi^{(\pm)}_{n}(t)$. Note that the operators $\hat{A}^{(\pm)}_{n}$ are proportional to the quadrature operators $\hat{x}_{\mathrm{c},n},\hat{p}_{\mathrm{c},n}$ of the cavity modes, i.e., $\hat{A}^{(+)}_{n}\propto\hat{x}_{\mathrm{c},n}$, $\hat{A}^{(-)}_{n}\propto\hat{p}_{\mathrm{c},n}$.

The Hamiltonian (11) has the form $\hat{H}=\hat{H}_{\textrm{full}}+\hat{H}_{\textrm{Dr}}$, where we have decided to define $\hat{H}_{\textrm{Dr}}=\sum_{n}\left[\hbar\,\xi^{(+)}_{n}\,\hat{A}^{(+)}_{n}+\hbar\,\xi^{(-)}_{n}\,\hat{A}^{(-)}_{n}\right]$ for convenience. The time-evolution operator $\hat{U}(t)$ induced by $\hat{H}$ has the expression (3), which can be manipulated to obtain the equivalent expression

[TABLE]

where $\hat{U}_{\textrm{full}}(t):=\overset{\leftarrow}{\mathcal{T}}\,\exp[-\frac{i}{\hbar}\int_{0}^{t}dt^{\prime}\,\hat{H}_{\textrm{full}}(t^{\prime})]$. Since $[\hat{H}_{p},\hat{H}_{p^{\prime}}]=0$, we obtain

[TABLE]

where we have defined the complex coupling $\xi_{n}:=\xi_{n}^{(+)}+i\,\xi_{n}^{(-)}$.

Unfortunately, we cannot proceed any further with simplifications of (13). The explicit expression in the time-ordered exponential leads to a time ordered exponential of exponential operators, which cannot be treated with the tools described here. We leave it to further work to study this important case in more detail.

IV Time evolution of quantities of interest

The decoupling achieved above allows us to obtain analytical control on the time evolution of the system. This includes our ability to compute explicitly the expectation value of many quantities of interest, which we proceed to do below.

IV.1 Time evolution of quantities of interest: mode operators

We start by computing the time evolution of the mode operators $\hat{a}_{k}$ and $\hat{b}_{p}$, through which one can compute the time evolution of most quantities of interest. Their time dependence is induced by the nonlinear Hamiltonian $\hat{H}_{\textrm{full}}$ and is obtained through the standard Heisenberg equation as $\hat{a}_{k}(t):=\hat{U}_{\textrm{full}}^{\dagger}\,\hat{a}_{k}\,\hat{U}_{\textrm{full}}$ and $\hat{b}_{p}(t):=\hat{U}_{\textrm{full}}^{\dagger}\,\hat{b}_{p}\,\hat{U}_{\textrm{full}}$ respectively. After some algebra we find

[TABLE]

where we have defined $F^{(p)}(t):=F^{(p)}_{-}(t)+i\,F^{(p)}_{+}(t)$ and $F^{(p)}_{n}(t):=F^{(p,+)}_{n}(t)+i\,F^{(p,-)}_{n}(t)$, and the notation $F^{(p)}_{\{nm\}}:=\frac{1}{2}[F^{(p)}_{nm}+F^{(p)}_{mn}]$.

IV.2 Final reduced state of the mechanical resonators

We continue by computing the final reduced state $\hat{\rho}_{\textrm{m}}(t)$ of the mechanical resonators. The reduced state is defined by $\hat{\rho}_{\textrm{m}}(t):=\textrm{Tr}_{\textrm{Phot}}(\hat{\rho}_{\textrm{NL}}(t))$, that is, by tracing over all of the cavity modes. In C we provide all of the detailed computations for this part. We assume that the initial state $\hat{\rho}_{\mathrm{c}}(0)$ of the cavity modes is separable from the initial state $\hat{\rho}_{\mathrm{m}}(0)$ of the mechanical modes. This implies that the full initial state is $\hat{\rho}_{0}=\hat{\rho}_{\mathrm{c}}(0)\otimes\hat{\rho}_{\mathrm{m}}(0)$. It is not difficult to show that the reduced state $\hat{\rho}_{\textrm{m}}(t)$ at time $t$ has the form

[TABLE]

where we have introduced $\sum_{\overset{\{n_{k}\}}{k\in\mathcal{I}}}:=\sum_{n_{1},n_{2},...,n_{N}}$ for $N$ modes that belong to the set of all possible combinations of excitations $\mathcal{I}$, while $\sum_{k\in\mathcal{I}}\,J_{k}:=J_{1}+J_{2}+...+J_{N}$ for any $k$-dependent quantities $J_{k}$. We have also introduced

[TABLE]

Note that we have $\textrm{Tr}(\hat{\rho}_{\textrm{m}}(t))=\sum_{\overset{\{n_{k}\}}{k\in\mathcal{I}}}p_{\{n_{k}\}}=1$ as expected.

IV.3 Mode population

The time evolution of the modes allows us to immediately compute the operators that “count” the number of excitations at any moment in time, namely $\hat{a}_{k}^{\dagger}(t)\hat{a}_{k}(t)$ and $\hat{b}_{k}^{\dagger}(t)\hat{b}_{k}(t)$. They read

[TABLE]

The number operator $\hat{a}_{k}^{\dagger}(t)\hat{a}_{k}(t)$ of each field mode is a conserved quantity as can be immediately seen from the Hamiltonian (1) (that is, it commutes with the whole Hamiltonian). However, the number operators $\hat{b}_{k}^{\dagger}(t)\hat{b}_{k}(t)$ of the resonators are not, and they depend on the nonlinear coupling through the functions $F^{(k)}$ and $F^{(k)}_{n}$.

We note here that, if the nonlinear couplings are small, i.e., they are proportional to $\epsilon\ll 1$, we have that $F^{(k)}\sim\epsilon$, $F^{(k)}_{n}\sim\epsilon$ and $F^{(k)}_{n}\,F^{(k)*}_{m}\sim\epsilon^{2}$. Therefore, the last term contributing to the number operator $\hat{b}_{k}^{\dagger}(t)\hat{b}_{k}(t)$ in (IV.3) is a negligible contribution in this regime. When this occurs, it is easy to check that the result is equivalent to what would be obtained through first-order perturbation theory, as expected. In this sense, the operators (IV.3) inform us on the full evolution of the resonator’s population and contain a signature of the full nonlinear character of the system.

IV.4 First-order bipartite quantum coherence

Given two modes $m$ and $n$, we call the correlation $\langle\hat{d}^{\dagger}_{m}\hat{d}_{n}\rangle$ the (first-order) bipartite coherence, sometimes denoted by $G^{(1)}_{mn}$ in optics Walls and Milburn (2008); Kalaga et al. (2016). This definition applies in general to any state. This measure corresponds to a simple interferometric setup, where we collect the photons in the modes $m$ and $n$, add a phase difference between their paths, and let them interfere. We will witness the formation of an interference pattern only if the quantity $\langle\hat{d}^{\dagger}_{m}\hat{d}_{n}\rangle$ is non-zero. This quantity can be normalized by the power in each mode, and in this case we recover the standard definition of first-order amplitude correlation function $g^{(1)}_{mn}$ from quantum optics applied to modes $m$ and $n$, namely

[TABLE]

It is not difficult to employ our results and compute $g^{(1)}_{mn}$ for pairs of cavity modes $\hat{a}_{k}$ and $\hat{a}_{k^{\prime}}$, for pairs of resonator modes $\hat{b}_{k}$ and $\hat{b}_{k^{\prime}}$, or for pairs of cavity and resonator modes $\hat{a}_{k}$ and $\hat{b}_{k^{\prime}}$. The results can be obtained analytically but are not illuminating and we omit to print the general formulas. Instead, we will give a few explicit results later on, when the initial state has been chosen.

IV.5 Mixedness and linear entropy

The final quantity that we are interested in computing is the amount of mixedness induced in the two main subsystems, i.e., the optical and mechanical one, due to the interaction between these systems. In particular, we would like to be able to isolate the contribution to any such induced mixedness and the nonlinear part of the interaction.

A measure of the mixedness of a state is the linear entropy $S_{N}$, defined as $S_{N}=1-\textrm{Tr}(\hat{\rho}^{2})$ and vanishes for pure states Peters, Wei, and Kwiat (2004).333Any pure state $\hat{\rho}$ satisfies $\hat{\rho}^{2}=\hat{\rho}$. Since $\textrm{Tr}(\hat{\rho})=1$ this implies that $\textrm{Tr}(\hat{\rho}^{2})=1$. This motivates the definition of the linear entropy. We also note that, contrary to the case of finite-dimensional systems with dimension $d$, where the state $\rho$ with maximal mixedness $S_{N}=1/d$ is the diagonal state with uniform eigenvalues $\lambda=1/d$, i.e., $\hat{\rho}=\sum_{n}\lambda|n\rangle\langle n|$ for an complete orthonormal basis $|n\rangle$, the case of infinite dimensional systems is more subtle. Clearly, there cannot exist a diagonal state with uniform eigenvalues $\lambda=1/d$, since d is infinite. However, one could have diagonal states with finite amount of uniform, non-zero eigenvalues $\lambda=1/k$, were $k$ can be arbitrary. In this case, one would still have $S_{N}=1/k$, which can be made arbitrarily small by increasing $k$. The time evolution of the system is expected to induce coherence between the photonic part and the mechanical part and the coherence will induce mixedness between these parts, which we can quantify using the linear entropy $S_{N}$. We focus on the reduced state $\hat{\rho}_{\textrm{m}}(t)$ of the mechanical subsystem obtain above (15). In our case we have

[TABLE]

where it is easy to check that $\hat{D}_{\{n_{k}\}}^{\dagger}\,\hat{D}_{\{m_{k}\}}=\prod_{p}e^{i\theta_{p}^{\prime}}\,\exp\left[i\,\Delta^{(p)}_{\{n_{k},m_{k}\}}\,\hat{b}^{\dagger}_{p}+\textrm{h.c.}\right]$. Here we have introduced the very useful functions $F^{(p)}_{k}:=F^{(p,+)}_{k}+i\,F^{(p,-)}_{k}$ and $\Delta^{(p)}_{\{n_{k},m_{k}\}}:=\sum_{k\in\mathcal{I}}\left(n_{k}-m_{k}\right)F^{(p)}_{k}$. The exact expression of the phase $e^{i\,\theta^{\prime}}$ is irrelevant since it clearly cancels out in (IV.5).

We know that when all the nonlinear couplings vanish, i.e., $g^{(\pm)}_{np}=0$, we have $F^{(p)}_{k}=0$. As an immediate consequence of this is that the expression (IV.5) reduces to $S_{N}(t)=1-\textrm{Tr}(\hat{\rho}^{2}_{\textrm{m}}(0))=S_{N}(0)$. Therefore, the mixedness introduced in the reduced state $\hat{\rho}_{\textrm{m}}$ of the resonators is a direct and only consequence of the nonlinear interaction. This analytical insight is, of course, perfectly in line with what is expected.

V Application to initial coherent state of the cavity modes and thermal state of the mechanical modes

Here we apply our results to a more concrete setup. We assume that there are a limited amount of modes $k\in\mathcal{I}$ that are initially in a coherent state $|\mu_{k}\rangle$ with parameter $\mu_{k}$ and defined by $\hat{a}_{k}|\mu_{k}\rangle=\mu_{k}\,|\mu_{k}\rangle$, while the cavity modes $s\notin\mathcal{I}$ are each in their respective vacuum state $|0\rangle_{s}$.

We also assume that the mechanical modes $\hat{b}_{p}$ are initially in a thermal state $\hat{\rho}_{\textrm{m}}(0)=\prod_{p\notin\mathcal{I}}\sum_{j_{p}}\frac{\tanh^{2\,j_{p}}(r_{p})}{\cosh^{2}(r_{p})}|j_{p}\rangle\langle j_{p}|$, with temperature $T$ and parameter $r_{p}$ defined by $\tanh(r_{p}):=\exp[-\frac{\hbar\,\omega_{\textrm{m},p}}{2\,k_{\mathrm{B}}\,T}]$. This is the standard initial setup in most applications, such as those with mechanical oscillators Teufel et al. (2011) and with levitated nano-objects Jain et al. (2016). Note that the set $\mathcal{I}$ might include any number $N$ of modes with $N\geq 1$, and that $N_{\textrm{i},p}\equiv\sinh^{2}(r_{p})$ is the initial population of thermal mechanical phonons.444The temperatures can be lowered to values that allow to reduce the number $N_{\textrm{i},p}$ of initial thermal phonons to $N_{\textrm{i},p}\sim 0.34$, in the case of mechanical oscillators Teufel et al. (2011). In the case of levitated nano-objects Jain et al. (2016), one can reach temperatures that give rise to an average population of $N_{\textrm{i},p}\sim 60$.

The initial state of the system $\hat{\rho}(0)$ is then separable in the mode/resonator bipartition, and has the expression

[TABLE]

Note that the state is mixed due to the initial temperature in the mechanical modes.

V.1 Final reduced state of the mechanical resonators

The final reduced state $\hat{\rho}_{\mathrm{m}}(t)$ of the mechanical modes has been computed for the general case and reads (15). In the specific case we are studying here, we have to use $p_{\{n_{k}\}}=\prod_{k\in\mathcal{I}}\frac{|\mu_{k}|^{2\,n_{k}}}{n_{k}!}\,e^{-|\mu_{k}|^{2}}$ and we will obtain the final expression. This expression is cumbersome and we avoid printing it here.

V.2 First-order bipartite coherence

In this case, we can provide some explicit formulas given that we have specified the initial state of the system. We use the expression (18) and we provide expressions for the nominator and denominator separately.

To compute the first-order coherence between two resonators $\hat{b}_{k}$ and $\hat{b}_{k^{\prime}}$ or a cavity mode $\hat{a}_{k}$ and resonator $\hat{b}_{k^{\prime}}$ we need the following on-diagonal expressions

[TABLE]

and off-diagonal expressions

[TABLE]

all of which have been computed using the quantities in B. In particular, to obtain the expressions (V.2) and (V.2) we faced no conceptual hurdles but only lengthy algebra. We do not provide all the steps of the computations here for the sake of clarity of presentation.

V.3 Mixedness of the final reduced state of the mechanical resonators

Given our chosen initial state $\hat{\rho}(0)$, and the final reduced state of the mechanical resonators $\hat{\rho}_{\textrm{m}}(t)$, we can compute the mixedness induced by the full evolution $\hat{U}_{\textrm{full}}$. This gives us the linear entropy for finite temperature, which has a simple and analytical expression

[TABLE]

We can now look at the contribution when the nonlinearity is switched off, i.e., $\left|\Delta^{(p)}_{\{n_{k},m_{k}\}}\right|=0$. This implies that (23) yields the mixedness $S^{\textrm{in}}_{N}$ of the initial state, which simply reads

[TABLE]

This allows us to express the full mixedness (23) as

[TABLE]

which is one of our main results.

The zero temperature $T=0$ case is simply obtained by setting $r_{p}=0$ for all $p$. When the temperature becomes increasingly high, it tends to inhibit the generation of mixedness, i.e., correlations between the two systems. This is in line with previous results that studied the competition between initial mixedness (due to temperature) and the coherent generation of excitations Bruschi et al. (2013).

It is now clear from (25) that the last fraction in the expression is the direct and full contribution of the nonlinearity to the mixedness, since it vanishes (together with the whole expression) for vanishing nonlinear coupling.

VI Applications: optomechanics

Now that we have obtained our results for a general setup of many cavity modes interacting with many mechanical resonators, we can apply them to simple scenarios that model cases of physical interest. In particular, we focus on Hamiltonians that describe more in particular standard closed optomechanical systems without external drive (i.e., we still consider the evolution through the nonlinear unitary operator $\xi_{p}^{(\pm)}=0$ for all $p$).

Let us assume that $\lambda_{p}^{(\pm)}=g^{(-)}_{np}=0$ and call $g^{(+)}_{np}(t)\equiv g_{np}(t)$, that is, we will consider the Hamiltonian (1) which now reads

[TABLE]

This implies also that $\tilde{F}_{\textrm{c},n}^{(p)}=F_{(\pm)}^{(p)}=0$, and therefore the time evolution operator (III.3) that here reduces to

[TABLE]

Note that, for each fixed oscillator $p$, this coincides exactly with the expression obtained in the literature Bruschi and Xuereb (2018).

VI.1 Linearised regime

Our techniques allow, in principle, to fully grasp the nonlinear character of the system. Previously, such analytical understanding was not available and a standard approach has been to linearise the standard optomechanical Hamiltonian, i.e., the Hamiltonian (26). This procedure is standard and we refer to the literature for an extensive justification Aspelmeyer, Kippenberg, and Marquardt (2014).

In brief, we consider the interaction term $\hbar\,g_{np}\hat{a}^{\dagger}_{n}\hat{a}_{n}\hat{B}^{(+)}_{p}$ and replacing $\hat{a}_{n}\rightarrow\alpha_{n}+\delta\hat{a}_{n}$, where $\alpha_{n}\gg 1$ represents a large mean value for the cavity operators and the “new” operator $\delta\hat{a}_{n}$ represents small excitations around the large classical mean $\alpha_{n}$, which we assume real for simplicity. The operator $\delta\hat{a}_{n}$ has vanishing first moment, i.e., $\langle\delta\hat{a}_{n}\rangle=0$. In addition, the linearisation procedure requires us to start in an initial state where $\mu_{n}=0$ for all $n$, that is., all cavity modes are initially in the vacuum. The final step is to keep only terms proportional to $\alpha_{n}$ and $\alpha_{n}^{2}$. What is left is called the linearised Hamiltonian $\hat{H}_{\textrm{lin}}$, which reads

[TABLE]

where it is evident that the Hamiltonian contains only terms that are at most quadratic in the creation and annihilation operators. In particular, the last term is the central contribution arising from the linearisation of the nonlinear part of the initial Hamiltonian (1). Here we also have $\hat{A}^{(+)}_{n}:=\delta\hat{a}_{n}^{\dagger}+\delta\hat{a}_{n}$

The time evolution induced by the linearised Hamiltonian (28) reads $\hat{U}_{\textrm{lin}}(t)=\hat{U}_{0}(t)\,\overset{\leftarrow}{\mathcal{T}}\,\exp[-i/\hbar\,\int_{0}^{t}\,dt^{\prime}\,\hat{H}_{1}]$, where we have introduced $\hat{U}_{0}(t)=\exp[-i(\sum_{n}\omega_{\textrm{c},n}\delta\hat{a}_{n}^{\dagger}\delta\hat{a}_{n}+\sum_{p}\omega_{\textrm{m},p}\,\hat{b}_{p}^{\dagger}\hat{b}_{p})\,t]$ and

[TABLE]

Decoupling of the time evolution $\hat{U}_{\textrm{lin}}$ induced by this Hamiltonian can in principle be done analytically Bruschi, Lee, and Fuentes (2013). It requires $(N+M)\,(2\,(N+M)+1)$ terms for the quadratic part and $2\,M$ for the order one part (i.e., the displacement part of the hamiltonian) in the decomposition (4). Here $N$ and $M$ indicate the number of cavity modes and mechanical modes respectively.

We do not want to proceed to obtain such solution, since it would require extremely cumbersome algebra which is unnecessary for our purposes here. Instead, in the following, we will use the full result and the linearised one to discuss the usefulness of the work done here.

VI.2 The weak coupling regime

We have obtained exact expressions for all quantities of interest in this work. We can now look at a regime of weak coupling between light and matter, that is, to assume that the couplings $g_{np}$ in the Hamiltonian (26) are small, i.e., $g_{np}(t)=\epsilon\,\tilde{g}_{np}(t)$, where $\epsilon\ll 1$ and $\tilde{g}_{np}(t)$ is some time-dependent function that remains finite at all times, such that $g_{np}^{(+)}(t)\ll 1$. We assume that the parameter $\epsilon$ is the same for all couplings for the sake of simplicity.

VI.2.1 The weak coupling regime: bipartite coherence

We now compute the bipartite coherence induced by the full nonlinear Hamiltonian and by the linearised Hamiltonian (VI.1).

For the full nonlinear setup we have

[TABLE]

while, in the linearised regime we find

[TABLE]

These two expressions are different. If coherence could be measured experimentally, we should find agreement with the expression obtained through the full nonlinear Hamiltonian.

VI.2.2 The weak coupling regime: mixedness of the mechanical resonator subsystem

We can now ask the following question: how much mixedness is induced if the nonlinear coupling is weak?

It is not difficult to show, using perturbation theory, that the linear entropy $S_{N}$ for the full nonlinear case and the linearized case reads

[TABLE]

This implies that, to the order we are interested in, the optical and mechanical subsystem remain separable.

We note here that there is coherence induced between the optical modes and the mechanical ones at first order in $\epsilon$ in this regime, as seen in (VI.2.1) and (VI.2.1). However, the mixedness induced between the states of the subsystems occurs only to second order in $\epsilon$. This is not an unknown feature of such properties of quantum states in perturbation theory Sabín et al. (2016).

VI.3 Applications to optomechanics: single mode & multi-resonator cavity

We specialise to standard optomechanical scenarios with one cavity mode $\hat{a}_{\tilde{k}}$ and an arbitrary number of mechanical resonators $\hat{b}_{p}$. This means that the set $\mathcal{I}=\{\tilde{k}\}$, which has one element.

The first quantities we can compute are the time evolution of the mode operators. Using the results obtained above we find

[TABLE]

It is then easy to check that, given our chosen initial state, we have $\langle\hat{a}_{\tilde{k}}^{\dagger}(t)\,\hat{a}_{\tilde{k}}(t)\rangle=\langle\hat{a}_{\tilde{k}}^{\dagger}(0)\,\hat{a}_{\tilde{k}}(0)\rangle$ and $\langle\hat{b}_{p}^{\dagger}(t)\,\hat{b}_{p}(t)\rangle=N_{\mathrm{i},p}+|F_{\tilde{k}}^{(p)}|^{2}\,(\mu_{\tilde{k}}^{2}+\mu_{\tilde{k}}^{4})$.

We can also compute the first-order bipartite quantum coherence (18) for the field mode and an oscillator which reads

[TABLE]

and for two oscillators, which gives us

[TABLE]

Above, we have defined the angle $\phi_{\tilde{k}}:=\frac{1}{2}\,\sum_{p}\,(F^{(p)}_{\{\tilde{k}\tilde{k}\}}+2\,F^{(p,+)}_{\tilde{k}}\,F^{(p,-)}_{\tilde{k}})$. Note that, for zero temperature we have $r_{p}=0$ for all $p$, and the expressions (34) and (35) simplify and reduce to

[TABLE]

and $g^{(1)}_{pp^{\prime}}(t)=1$. This means that, while the mode of light and any resonator mode are coherent with a strength that depends on the parameters of the problem, any pair of resonators is perfectly coherent at zero temperature. We also note that, in this regime, the light and resonator coherence decreases exponentially with $|\mu_{\tilde{k}}|$, unless $\phi_{\tilde{k}}=0$. In the limit $|\mu_{\tilde{k}}|\rightarrow\infty$ we have that $g^{(1)}_{\tilde{k}p^{\prime}}(t)=\exp[-\frac{1}{2}\,\sum_{p}\,|F_{\tilde{k}}^{(p)}|^{2}]\neq 0$ only at the times $t_{n}$ such that $\phi_{\tilde{k}}(t_{n})=0$ and vanishes for all other times. This implies that, to be able to verify the coherence between light and a single resonator it is necessary to reduce the number of photons in the coherent state as much as possible.

For finite temperature the issue becomes more delicate. It is clear from (34) and (35) that an increase in temperature (i.e., an increase in $r_{p}$ for all $p$) implies that it is more difficult to establish the desired coherence. Therefore, reducing the temperature to levels where the initial phononic population $N_{\mathrm{i},p}$ fore each resonator becomes small is paramount.

We can also compute the mixedness of the reduced state of the oscillators, which employs algebraic manipulations that can be found in D.1. Finally, we show that it reads

[TABLE]

In this formula we have introduced the modified Bessel functions $I_{n}(z)$. Notice that, when $F^{(p)}_{\tilde{k}}=0$ we recover immediately $S_{N}=S^{\textrm{in}}_{N}$ as expected.

VI.4 Modulated scenario for a single mode & multi-resonator cavity

We can now ask another question of potential physical interest, which can also be used to highlight the difference in the predictions obtained with the full time-evolution obtained here, and with the linearised approach. We consider a scenario where the couplings have a time dependence of the form $g_{\tilde{k}p}(t)=g_{\tilde{k}p}\,(1+\kappa\,\sin(\omega_{\mathrm{d}}\,t))$ or $g_{\tilde{k}p}(t)=g_{\tilde{k}p}\,(1+\kappa\,\cos(\omega_{\mathrm{d}}\,t))$ for some modulating frequency $\omega_{\mathrm{d}}$ and amplitude of the oscillation $\kappa$.

This form of the coupling can be used to model in a simple way an interaction strength between light and matter that is not constant, something that is typically assumed as the lowest order approximation in many experiments. However, no quantity of this kind can be fundamentally and ultimately constant. Therefore, the expression that we provide here is provides a reasonable toy model to study deviations from the ideal case, where the mean value $g_{\tilde{k}p}$ is the one measured in the laboratory. This form of coupling has already been proposed and used in for different studies Qvarfort et al. (2018).

Here we choose to look at two possible situations: $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}\pm\omega_{\textrm{m},\tilde{p}}$. Here $\tilde{k}$ and $\tilde{p}$ are a specific optical and mechanical mode respectively. In E we show that, depending on the choice, in the linearised scenario this corresponds to a two-mode squeezing operation (+ sign) or mode-mixing operation (- sign) after a sufficiently long time (i.e., an approximation that is also used in the rotating wave approximation). Here we do not discuss the details of how much time needs to pass for this approximation to become more and more accurate but we just note that such resonant regime naturally occurs across all of physics and we leave it to future work to study the details of potential implementations.

VI.5 Modulated scenario: linearised results

Let us start with the linearised case. When $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}+\omega_{\textrm{m},\tilde{p}}$ we find that the evolution of the operators $\delta\hat{a}_{\tilde{k}}$ and $\hat{b}_{\tilde{p}}$ in the linearised regime after sufficiently long time reads

[TABLE]

while, when $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}-\omega_{\textrm{m},\tilde{p}}$ we have, after sufficiently long time,

[TABLE]

These expressions allow us to compute the expectation values of the population operators in these regimes. Recalling that we defined $\hat{a}_{k}=\alpha_{k}+\delta\hat{a}_{k}$, for $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}+\omega_{\textrm{m},\tilde{p}}$ we have

[TABLE]

while for $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}-\omega_{\textrm{m},\tilde{p}}$ we have

[TABLE]

Notice that both the resonant and asymptotic (i.e., after a sufficiently long time) behaviors (VI.5) and (VI.5) must be compatible with the fact that we have linearised the system, that is, $\delta\hat{a}$ is an operator that creates excitations around the vacuum. Therefore, we believe the validity of these expressions for “large” times, which do not invalidate the linearisation procedure. The exact details are left for future work.

VI.6 Modulated scenario: full results

In the case of the results obtained in this work, we note immediately through the expressions (III.3) and (VI.3) that we do not have any resonant behavior when $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}\pm\omega_{\textrm{m},\tilde{p}}$. In fact, the expressions for the time evolution operator (VI.3) imply that any resonance would occur for $\omega_{\mathrm{d}}=\omega_{\textrm{m},\tilde{p}}$, which is in stark contrast with the results obtained above and displayed in (VI.5) and (VI.5).

In particular, we obtain

[TABLE]

If we drive the system $\omega_{\mathrm{d}}=\omega_{\textrm{m},\tilde{p}}$, after a sufficiently long time we would obtain

[TABLE]

VI.7 Modulated scenario: considerations

In any realistic implementation one has $\omega_{\textrm{c},\tilde{k}}\gg\omega_{\textrm{m},\tilde{p}}$, and therefore we see that the the combination $\omega_{\textrm{c},\tilde{k}}\pm\omega_{\textrm{m},\tilde{p}}$ would never take values very close to $\omega_{\textrm{m},\tilde{p}}$.

For these reasons, the resonant regimes can be distinguished through experimental efforts, by driving the system at completely different drive frequencies. Our results show that, by driving the system with a coupling of the form $g\propto(1+\kappa\,\sin(\omega_{\mathrm{d}}\,t))$, where the drive frequency matches that of one of the mechanical resonators, we should observe a time dependent population of phonons predicted by (VI.6). By driving the system at the frequency $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}-\omega_{\textrm{m},\tilde{p}}$, we should not observe any resonant behavior, contrary to what predicted by the linearised regime in (VI.5).

VII Conclusion

In this work we found an analytic expression for the time evolution of an arbitrary number of coupled bosonic modes interacting through a time-dependent optomechanical-like Hamiltonian. Despite of the system having potentially an arbitrarily large amount of constituents, and despite of the nonlinearity present in the interaction driven by a time-dependent coupling, we were able to decouple the time-evolution operator using tools developed for this purpose Bruschi, Lee, and Fuentes (2013); Bruschi and Xuereb (2018). This result does not rely on any approximation and is therefore general: the only assumption made is that the Hamiltonian has the form (1), which occurs in many systems, such as optomechanical cavities Aspelmeyer, Kippenberg, and Marquardt (2014).

We were able to compute the time evolution of the expectation values of meaningful operators, such as the average photonic and phononic excitation, the first-order quantum bipartite coherence, and the mixedness of the reduced state of the mechanical oscillators. Our results allow us to study the coherence induced between the subsystems due to the nonlinear interaction of light and matter. Furthermore, they allow us to clearly quantify, using the linear entropy, the increase of mixedness of the subsystem of the resonators when the nonlinearity is switched on, as a function of time, the initial photonic population, and temperature inside the cavity. In addition, we compared some of the predictions with those arising from the linearised version of our system. Linearisation is typically used to describe optomechanical systems, and we have found that there are qualitative differences that arise when quantifying, for example, the expectation value of the cavity mode population or of the phononic population of the mechanical excitations. We gave a concrete example where the light matter coupling can be modulated periodically, and we were able to show that the linearised model and the full nonlinear one treated here, give different results.

Given the lack of a systematic understanding of the nonlinear nature of the interaction, and the nonlinearity induced by the coupling, this analytical insight can help shedding light onto intrinsic nonlinear aspects of quantum (opto) mechanical systems. In the end, these insights can also help in the quest of demonstrating in the laboratory the quantum nature of “macroscopic objects”, such as the mechanical resonators. Finally, the decoupling obtained here can be applied to many situations of theoretical and practical interest. We leave it to future work to pursue such new directions.

Acknowledgments

We acknowledge Alessio Serafini and André Xuereb for useful comments and discussions.

Appendix A Decoupling of the Hamiltonian

Here we show how to decouple the time-evolution operator $\hat{U}$ induced by the Hamiltonian (6), reprented here

[TABLE]

The procedure to be followed has been developed in the literature Bruschi, Lee, and Fuentes (2013). All details can be found there. The time-evolution operator is defined by $\hat{U}:=\overset{\leftarrow}{\mathcal{T}}\,\exp[-i\,\int_{0}^{t}\,dt^{\prime}\,\hat{H}(t^{\prime})]$.

As discussed in the main text, we make the decoupling ansatz

[TABLE]

where we have defined the Lie-algebra basis operators

[TABLE]

We now take the time derivative on both sides of the expression (45) and then multiply on the right by $\hat{U}^{\dagger}(t)$ to obtain

[TABLE]

We then use similarity relations of the form

[TABLE]

to obtain an explicit expression for (A) and, equating coefficients on both sides we obtain the differential equations for the time dependent $F$-functions in terms of the coefficients of the Hamiltonian, which read

[TABLE]

Note that the similarity relations (A) are not exhaustive nor complete. They are presented to give an intuitive understanding of the type of relations needed to simplify the expression (A). A complete list can be found in Bruschi and Xuereb (2018).

Appendix B Some useful expressions

In this section we provide some useful expression used in this work.

Let us start by the first, defined as $I_{\alpha}:=\langle\mu|\exp[-i\,\alpha\,\hat{a}^{\dagger}\,\hat{a}]\,|\mu\rangle$, where $\hat{a}|\mu\rangle=\alpha\,|\mu\rangle$. We have

[TABLE]

Notice that the expression $I^{(n)}_{\alpha}:=\langle\mu|\,(\hat{a}^{\dagger}\,\hat{a})^{n}\,\exp[-i\,\alpha\,\hat{a}^{\dagger}\,\hat{a}]\,|\mu\rangle$ can be computed as $\tilde{I}_{\alpha}=i^{n}\,\frac{d^{n}}{d\alpha^{n}}I_{\alpha}$. Therefore

[TABLE]

We continue by computing the expression $J_{\alpha}:=\sum_{n=0}^{+\infty}T^{2n}/C^{2}\langle n|\exp[\alpha\,\hat{a}^{\dagger}-\alpha^{*}\,\hat{a}]\,|n\rangle$, where $T:=\tanh r$ and $C:=\cosh r$. We have

[TABLE]

which gives us the final, simple and compact result

[TABLE]

In the above computations we have introduced the Confluent hypergeometric function ${}_{1}F_{1}(a,b;z)$ and the Laguerre polynomial $L_{n}(z)$. We have used the fundamental property ${}_{1}F_{1}(a,b;z)=\exp[z]\,{}_{1}F_{1}(b-a,b;-z)$ and the relation $L_{n}(z)={}_{1}F_{1}(-n,1;z)$. We have used the generating function for the Laguerre polynomial to be able to go to the last line of (B) from the second-to-last one.

We are also interested in $\tilde{J}_{\alpha}:=\sum_{n=0}^{+\infty}T^{2n}/C^{2}\langle n|\exp[\alpha\,\hat{a}^{\dagger}-\alpha^{*}\,\hat{a}]\,\hat{a}\,|n\rangle$. We can proceed as above and find

[TABLE]

which gives us the final, simple and compact result

[TABLE]

Note that here we have introduced the generalised Laguerre polynomials $L^{(q)}_{a}(z)$. Here as well, we have used the generating function for the generalised Laguerre polynomial to go to the last line of (B) from the second-to-last one.

Finally, we want to compute the following $\tilde{L}_{\alpha}:=\sum_{l,l^{\prime}=0}^{+\infty}T^{2(l+l^{\prime})}/C^{4}\,\left|\langle l|\exp[\alpha\,\hat{a}^{\dagger}-\alpha^{*}\,\hat{a}]|l^{\prime}\rangle\right|^{2}$. We once more proceed using similar techniques as above, skip some passages and find

[TABLE]

where ${}_{1}F_{1}(a,b;x)$ is the Confluent hypergeometric function that satisfies ${}_{1}F_{1}(a,b;x)=\exp[x]\,{}_{1}F_{1}(b-a,b;-x)$, and the functions $L_{a}(x)$ and $L^{(q)}_{a}(x)$ are the Laguerre and generalised Laguerre polynomials respectively.

We now use the Hardy-Hille formula for Laguerre polynomials and obtain

[TABLE]

Finally, we use the fundamental Jacobi-Anger expansion, which can be recast as the identity $\exp[z\,\cos\theta]=I_{0}(z)+2\,\sum_{n=1}^{\infty}I_{n}(z)\,\cos(n\,\theta)$ for modified Bessel functions, and note that in our case we have $z=2\frac{T^{2}}{1-T^{4}}\,\left|\alpha\right|^{2}$ and $\theta=0$. Therefore

[TABLE]

Appendix C Reduced state of the resonators

We want to compute the final reduced state $\hat{\rho}_{\textrm{m}}$ of the mechanical resonators. The reduced state is defined by $\hat{\rho}_{\textrm{m}}:=\textrm{Tr}_{\textrm{Phot}}(\hat{\rho}_{\textrm{NL}}(t))$. We assume that the initial state is $\hat{\rho}_{0}=\hat{\rho}_{\mathrm{c}}(0)\otimes\hat{\rho}_{\mathrm{m}}(0)$. We have

[TABLE]

where we have introduced $\sum_{\overset{\{n_{k}\}}{k\in\mathcal{I}}}:=\sum_{n_{1},n_{2},...,n_{N}}$ for $N$ modes that belong to the set of all possible combinations of excitations $\mathcal{I}$, while $\sum_{k\in\mathcal{I}}\,P_{k}:=P_{1}+P_{2}+...+P_{N}$ for any $k$-dependent quantities $P_{k}$. We have also introduced

[TABLE]

Finally, we have also introduced $F^{(p)}:=F^{(p)}_{+}+i\,F^{(p)}_{-}$ and $F^{(p)}_{k}:=F^{(p,+)}_{k}+i\,F^{(p,-)}_{k}$. The phases $e^{i\theta_{p}}$ are irrelevant since they cancel out in the expression (15). Note that we have $\textrm{Tr}(\hat{\rho}_{\textrm{m}}(t))=\sum_{\overset{\{n_{k}\}}{k\in\mathcal{I}}}p_{\{n_{k}\}}=1$ as expected.

Appendix D Computing mixedness for initial coherent/thermal states of field modes/resonators

We want to compute an expression for $\text{Tr}_{\{m_{k}\},\{n_{k}\}}$ which is defined as

[TABLE]

where

[TABLE]

and $\Delta^{(p)}_{\{n_{k},m_{k}\}}:=\sum_{k\in\mathcal{I}}\left(n_{k}-m_{k}\right)F^{(p)}_{k}$ for simplicity of presentation. Again, the phases $e^{i\theta_{p}^{\prime}}$ are irrelevant.

Our goal can be easily reached using the useful result (58). We first note that we have to consider many $\tilde{L}_{\Delta^{(p)}_{\{n_{k},m_{k}\}}}$ and $r_{p}$. We then note the important expression

[TABLE]

This easily gives us

[TABLE]

D.1 Computing mixedness for initial coherent/thermal states of field modes/resonators

Here we need to compute the mixedness (25). This requires us to compute the simpler contribution

[TABLE]

The calculations follow here

[TABLE]

Using the Jacobi-Anger expansion for the modified Bessel functions we obtain

[TABLE]

Appendix E Time evolution in the periodically driven linearised regime

Here we compute the time evolution induced by the linearized Hamiltonian $\hat{H}_{\textrm{lin}}$ when only one cavity mode $\tilde{k}$ is present. It reads

[TABLE]

This Hamiltonian induced the time evolution through the operator

[TABLE]

where we have defined

[TABLE]

The Hamiltonians in (E) read

[TABLE]

Here we will assume that the couplings $g_{\tilde{k}p}(t)$ are modulated periodically. In this work, the choice of modulation is not restricted by any physical constraint, therefore we can choose $g_{\tilde{k}p}(t)=g_{\tilde{k}p}\,(1+\kappa\,\sin(\omega_{\mathrm{d}}\,t))$ or $g_{\tilde{k}p}(t)=g_{\tilde{k}p}\,(1+\kappa\,\cos(\omega_{\mathrm{d}}\,t))$ as simple modulations. It is easy to see that the choice of modulation will “select” a particular term in the expression for $\hat{H}_{1}(t)$. What we mean by this is that

[TABLE]

which have particular, or resonant, behavior when $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}\pm\omega_{\textrm{m},p}$. Assuming that all frequencies are positive and that $\omega_{\textrm{c},\tilde{k}}>\omega_{\textrm{m},p}$, we then have

[TABLE]

when $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}\pm\omega_{\textrm{m},p}$.

In the following we will consider both cases separately.

E.1 Time evolution in the periodically driven linearised regime: mode-mixing drive

Here we consider the case where $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}-\omega_{\textrm{m},\tilde{p}}$ for a specific oscillator $\tilde{p}$ and a drive of the form $g_{\tilde{k}p}(t)=g_{\tilde{k}p}\,(1+\kappa\,\sin(\omega_{\mathrm{d}}\,t))$. This means that, after long enough time, the effective Hamiltonian contributions will be coming from the term

[TABLE]

This implies that \hat{U}_{1}(t)\sim\exp\bigr{[}\alpha_{\tilde{k}}\,\kappa\,g_{\tilde{k}\tilde{p}}\,\bigl{(}\delta\hat{a}^{\dagger}_{\tilde{k}}\hat{b}_{\tilde{p}}-\delta\hat{a}_{\tilde{k}}\hat{b}^{\dagger}_{\tilde{p}}\bigr{)}\bigl{]} after enough time.

This also means that it is easy to compute the time evolution through the linearised Hamiltonian in this regime, which after enough time reads

[TABLE]

Notice that, as expected for a resonant drive scenario, the mode-mixing angle $\theta_{\textrm{m},\tilde{p}}$ between the cavity mode and resonator, reads $\theta_{\textrm{m},\tilde{p}}=\alpha_{\tilde{k}}\,\kappa\,g_{\tilde{k}\tilde{p}}\,t$ and increases with time.

E.2 Time evolution in the periodically driven linearised regime: squeezing drive

Here we consider the case where $\omega_{\mathrm{d}}=\omega_{\textrm{c},\tilde{k}}+\omega_{\textrm{m},\tilde{p}}$ for a specific oscillator $\tilde{p}$ and a drive of the form $g_{\tilde{k}p}(t)=g_{\tilde{k}p}\,(1+\kappa\,\sin(\omega_{\mathrm{d}}\,t))$. This means that, after long enough time, the effective Hamiltonian contributions will be coming from the term

[TABLE]

This implies that \hat{U}_{1}(t)\sim\exp\bigl{[}\alpha_{\tilde{k}}\,\kappa\,g_{\tilde{k}\tilde{p}}\,\bigl{(}\delta\hat{a}^{\dagger}_{\tilde{k}}\hat{b}^{\dagger}_{\tilde{p}}-\delta\hat{a}_{\tilde{k}}\hat{b}_{\tilde{p}}\bigr{)}\bigr{]} after enough time.

This also means that it is easy to compute the time evolution through the linearised Hamiltonian in this regime, which after enough time reads

[TABLE]

Notice that, as expected for a resonant drive scenario, the squeezing $r_{\textrm{m},\tilde{p}}$ between the cavity mode and resonator, reads $r_{\textrm{m},\tilde{p}}=\alpha_{\tilde{k}}\,\kappa\,g_{\tilde{k}\tilde{p}}\,t$ and increases with time.