Time evolution of two harmonic oscillators with cross-Kerr interactions

TL;DR

This paper provides an exact analytical solution for the time evolution of two coupled quantum harmonic oscillators with cross-Kerr interactions, enabling detailed analysis of their nonlinear dynamics without approximations.

Contribution

It introduces a method to decouple the evolution operator and derive explicit solutions for a nonlinear two-oscillator system, advancing understanding of full nonlinearity effects.

Findings

Exact analytical expressions for system evolution

Ability to analyze nonlinear effects without approximations

Potential applications in optomechanics

Abstract

We study the time evolution of two coupled quantum harmonic oscillators interacting through nonlinear optomechanical-like Hamiltonians that include cross-Kerr interactions. We employ techniques developed to decouple the time-evolution operator and obtain the analytical solution for the time evolution of the system. We apply these results to obtain explicit expressions of a few quantities of interest. Our results do not require approximations and therefore allow to study the nature and implications of the full nonlinearity of the system. Potential applications and extensions to optomechanics are also discussed.

Equations

$${\hat{H}}_{\text{I}}(t)=$$

$$+\hslash g_2^{(+)}(t){\hat{a}}^{†}\hat{a}{\hat{B}}_{+}^{(2)}+\hslash g_2^{(-)}(t){\hat{a}}^{†}\hat{a}{\hat{B}}_{-}^{(2)}$$

$$+2\hslash g_2^{\prime }(t){\hat{a}}^{†}\hat{a}\hat{b}{}^{†}\hat{b},$$

$${\hat{H}}_{\text{I}}(\tau )=$$

$$+{\tilde{g}}_2^{(+)}(\tau ){\hat{a}}^{†}\hat{a}{\hat{B}}_{+}^{(2)}+{\tilde{g}}_2^{(-)}(\tau ){\hat{a}}^{†}\hat{a}{\hat{B}}_{-}^{(2)}.$$

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

$$\hat{U}(t)=\prod _n{\hat{U}}_n(t),$$

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

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

$${\hat{H}}_{\mathrm{sq}}(\tau )$$

$$+a^{†}\hat{a}\left[{\tilde{g}}_2^{(+)}(\tau )\sin (2{\hat{\theta }}_2)+{\tilde{g}}_2^{(-)}(\tau )\cos (2{\hat{\theta }}_2)\right]{\hat{B}}_{-}^{(2)}$$

$${\hat{H}}_{\mathrm{OM}}(\tau )$$

$$+a^{†}\hat{a}\left[{\tilde{g}}_1^{(+)}(\tau )\sin ({\hat{\theta }}_2)+{\tilde{g}}_1^{(-)}(\tau )\cos ({\hat{\theta }}_2)\right]{\hat{B}}_{-}.$$

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

$${\hat{U}}_{\mathrm{OM}}:=$$

$${\hat{U}}_{\mathrm{sq}}^{†}\hat{b}{\hat{U}}_{\mathrm{sq}}=$$

$${\hat{U}}_{\mathrm{sq}}^{†}{\hat{H}}_{\mathrm{OM}}(\tau ){\hat{U}}_{\mathrm{sq}}=$$

$$-i{\hat{a}}^{†}\hat{a}\hat{\Im }\left[{\tilde{g}}_1(\tau ){\hat{E}}_2({\hat{\alpha }}^{†}-{\hat{\beta }}^{†})\right]{\hat{B}}_{-},$$

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

$$\times e^{-i{\int }_0^{\tau }d{\tau }^{\prime }{\hat{a}}^{†}\hat{a}\Re \left[{\tilde{g}}_1({\tau }^{\prime }){\hat{E}}_2({\hat{\alpha }}^{†}+{\hat{\beta }}^{†})\right]{\hat{B}}_{+}}$$

$$\times e^{i{\int }_0^{\tau }d{\tau }^{\prime }{\hat{a}}^{†}\hat{a}\Im \left[{\tilde{g}}_1({\tau }^{\prime }){\hat{E}}_2({\hat{\alpha }}^{†}-{\hat{\beta }}^{†})\right]{\hat{B}}_{-}},$$

$${\hat{F}}^{(2)}=$$

$$\times {\int }_0^{{\tau }^{\prime }}d{\tau }^{\prime \prime }\Re \left[{\tilde{g}}_1({\tau }^{\prime \prime }){\hat{E}}_2({\hat{\alpha }}^{†}+{\hat{\beta }}^{†})\right].$$

$$\hat{\alpha }=$$

$$\hat{\beta }=$$

$${\tilde{\chi }}_2{\ddot{\hat{p}}}_{11}+2i\frac{d}{d\tau }({\varphi }_2+{\hat{\theta }}_2){\dot{\hat{p}}}_{11}-{\tilde{\chi }}_2(a^{†}\hat{a}{)}^2{\hat{p}}_{11}=$$

$$\hat{b}(\tau )=$$

$$N_{\text{b}}(\tau )=$$

$$+|{\beta }_n{|}^2-2n^2\Re \left({\alpha }_n{\beta }_n^{\ast }{I}_n^2\right)].$$

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

$$\times e^{-i{\int }_0^{\tau }d{\tau }^{\prime }{\hat{a}}^{†}\hat{a}\left[{\tilde{g}}_1^{(+)}(\tau )\cos ({\hat{\theta }}_2)-{\tilde{g}}_1^{(-)}(\tau )\sin ({\hat{\theta }}_2)\right]{\hat{B}}_{+}}$$

$$\times e^{i{\int }_0^{\tau }d{\tau }^{\prime }{\hat{a}}^{†}\hat{a}\left[{\tilde{g}}_1^{(+)}(\tau )\sin ({\hat{\theta }}_2)+{\tilde{g}}_1^{(-)}(\tau )\cos ({\hat{\theta }}_2)\right]{\hat{B}}_{-}}.$$

$${\hat{F}}^{(2)}=$$

$$\times {\int }_0^{{\tau }^{\prime }}d{\tau }^{\prime \prime }\left[{\tilde{g}}_1^{(+)}({\tau }^{\prime \prime })\cos ({\hat{\theta }}_2)-{\tilde{g}}_1^{(-)}({\tau }^{\prime \prime })\sin ({\hat{\theta }}_2)\right].$$

$$\hat{b}(\tau )=$$

$$\Delta N_{\text{b}}(\tau )=$$

$$S_N({\hat{\rho }}_{\text{m}}(\tau ))=$$

$${\tilde{\chi }}_2(\tau ){\ddot{p}}_{11}-i\left(1+\frac{d}{d\tau }{\varphi }_2\right){\dot{p}}_{11}-{\tilde{\chi }}_2(\tau )(a^{†}\hat{a}{)}^2{\hat{p}}_{11}=$$

Full text

I Introduction

The understanding of the dynamics of quantum systems is paramount to deepen our understanding of the the laws of Nature. Given an arbitrary Hamiltonian, it is in general impossible to obtain a full and analytical expression for the time evolution of the system it describes. Typically, approximations or ad-hoc methods are required, leaving open the question of the existence of a general framework or set of tools to address the problem systematically.

In recent years, specific techniques have been developed to tackle the problem for Gaussian states of bosonic systems interacting through quadratic Hamiltonians Puri (2001); Bruschi et al. (2013). This has led to the development of a mathematical framework aimed at obtaining the full time evolution of the system using the power of the covariance matrix formalism Adesso et al. (2014). Interestingly, these techniques do not apply to this class of states and transformations only, but it can be generalised to arbitrary states and arbitrary Hamiltonians. This novel approach has led to very interesting preliminary results, where complete and analytical solutions of the time evolution of the system were found for optomechanical-like Hamiltonians of arbitrary numbers of modes of light interacting with arbitrary numbers of mechanical resonators Bruschi and Xuereb (2018); Bruschi (2018). This has allowed, for example, for a quantitative analysis of the nonlinear character of an optomechanical system Qvarfort et al. (2018).

In this work we take another step forward and study the time evolution of a system of two bosonic modes interaction via an optomechanical-like Hamiltonian with cross Kerr terms. Such Hamiltonian does not model a realistic optomechanical systems since it does not take into account external drive, loss or dissipation. However, this Hamiltonian models ideal quantum systems interacting through nonlinear terms that can potentially be designed in realistic implementations where loss and decoherence can be ignored at first. We are able to find an analytical solution to the time evolution which we use to compute the expectation value of relevant quantities as a function of time, such as the number of ‘phonons’ in the system, as well as the mixedness induced by the interaction in the state of the resonator’s subsystem. We specialise our results to a few scenarios of interest, which model light-matter interaction that can be exploited for sensing and gravimentry Rätzel et al. (2018), as well as tests of fundamental physics Schmöle et al. (2016).

Ultimately, our work can aid in deepening our understanding of the mathematical framework of quantum mechanics, one of our two pillars of modern science.

The paper is structured as followed: In Section II we introduce the general framework necessary for our work. In Section III we decouple the time evolution of the system for the general case. In Section IV we specialise to scenarios of more practical interest. Finally, we briefly give concluding remarks in Section V.

II General framework

We start by introducing the Hamiltonian and the tools necessary to this work. More details about the computations are left to the appendices.

II.1 Hamiltonian

In this work we consider the general Hamiltonian $\hat{H}=\hat{H}_{0}(t)+\hat{H}_{\textrm{I}}(t)$ where $\hat{H}_{0}:=\hbar\,\omega_{\mathrm{c}}(t)\hat{a}^{\dagger}\hat{a}+\hbar\,\omega_{\mathrm{m}}\,\hat{b}^{\dagger}\hat{b}$ is the potentially time-dependent free Hamiltonian and

[TABLE]

we have defined $\hat{B}_{+}:=\hat{b}^{\dagger}+\hat{b}$, $\hat{B}_{-}:=i(\hat{b}^{\dagger}-\hat{b})$, $\hat{B}^{(2)}_{+}:=\hat{b}^{\dagger 2}+\hat{b}^{2}$ and $\hat{B}^{(2)}_{-}:=i(\hat{b}^{\dagger 2}-\hat{b}^{2})$ for notational convenience. These operators correspond to the quadrature operators, and the square of the quadrature operators, respectively. We choose to retain our notation because it is more natural to the techniques used in this work. An advantage of the expression (II.1) is that it allows us to consider each cross-Kerr term independently from each other.

The Hamiltonian (II.1) is a formal extension of well-known Hamiltonians that can model light and matter interaction, e.g., optomechanical Hamiltonians. The main difference is that we consider an “ideal” case here, namely, we do not include an external (laser) drive and we assume that the system is lossless and noiseless. In this sense, the system considered here does not model realistic optomechanical systems. Nevertheless, it provides a platform for studying cross-Kerr interactions, which can become of practical interest with the advance of technological control.

II.2 Dimensionless dynamics

To understand which are the relevant dimensionless parameters that govern the dynamics of the system we introduce dimensionless quantities and rescale the Hamiltonian (II.1). This is achieved rescaling all frequencies and time by $\omega_{\mathrm{m}}$. Therefore, the laboratory time $t$ will become $\tau=\omega_{\mathrm{m}}\,t$, where $\tau$ is the new dimensionless time. The couplings in the Hamiltonian are subsequently relabelled as follows: $\tilde{g}^{(\pm)}_{1}(\tau)=g^{(\pm)}_{1}(\omega_{\mathrm{m}}\,t)/\omega_{\mathrm{m}}$, $\tilde{g}^{(\pm)}_{2}(\tau)=g^{(\pm)}_{2}(t\omega_{\mathrm{m}})/\omega_{\mathrm{m}}$ and $\tilde{g}^{\prime}_{2}(\tau)=g^{\prime}_{2}(t\omega_{\mathrm{m}})/\omega_{\mathrm{m}}$. We also rescale the optical frequency to $\tilde{\omega}_{\mathrm{c}}(\tau):=\omega_{\mathrm{c}}(\omega_{\mathrm{m}}\,t)/\omega_{\mathrm{m}}$ and the Hamiltonian by $\hbar$. This implies that we use the rescaled Hamiltonian $\hat{H}(\tau)=\hat{H}_{0}(\tau)+\hat{H}_{\textrm{I}}(\tau)$, where we have $\hat{H}_{0}(\tau):=\tilde{\omega}_{\mathrm{c}}(\tau)\,\hat{a}^{\dagger}\hat{a}+\hat{b}{}^{\dagger}\hat{b}$ and

[TABLE]

II.3 Solving the dynamics

Here we outline the tools needed to solve the dynamics generated by (II.2). We refer the reader to Appendices A and B for detailed calculations.

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

[TABLE]

where $\overset{\leftarrow}{\mathcal{T}}$ is the time ordering operator Bruschi et al. (2013).111An alternative approach was attempted in Brown et al. (2013) This expression simplifies to $U(t)=\exp[-\frac{i}{\hbar}\,\hat{H}\,t]$ when the Hamiltonian $\hat{H}$ is time independent. However, we are interested in Hamiltonians with time dependent parameters.

Any Hamiltonian can be cast in the form $\hat{H}(t)=\sum_{n}\hbar\,g_{n}(t)\,\hat{G}_{n}$, where the $\hat{G}_{n}$ are time-independent Hermitian operators, the $g_{n}(t)$ are time-dependent real functions, and the choice of $\hat{G}_{n}$ is not unique.

In general, one aims to recast (3) in the form

[TABLE]

where we have defined $\hat{U}_{n}:=\exp[-i\,F_{n}(t)\,\hat{G}_{n}]$ and the real functions $F_{n}(t)$ are in general time dependent. If such an expression exists, we say that the time-evolution operator has been decoupled. The functions $F_{n}(t)$ can be found using the techniques developed in the literature Bruschi et al. (2013) and are determined solely by the parameters of the Hamiltonian. The order of the operators in (4) is not unique. A different order will change the form of the functions $F_{n}(t)$ but not the expectation value of measurable quantities. A more detailed outline of these techniques can be found in Appendix A.

II.4 Initial state

The last ingredient in this work is the choice of the initial state of the system. We assume that the optical mode is initially in a coherent state $\ket{\mu_{\mathrm{c}}}$, while the mechanical is in a thermal state $\hat{\rho}_{\mathrm{m}}(T)$. These states are defined by $\hat{a}\ket{\mu_{\mathrm{c}}}=\mu_{\mathrm{c}}\ket{\mu_{\mathrm{c}}}$ and $\hat{\rho}_{\mathrm{m}}(T)=\sum_{n}\frac{\tanh^{2n}r}{\cosh^{2}r}\ket{n}\bra{n}$, and $\tanh r:=\exp[-\frac{\hbar\,\omega_{\mathrm{m}}}{2\,k_{\mathrm{B}}\,T}]$.

Therefore, the initial state $|\psi(0)\rangle$ reads

[TABLE]

where $N_{\text{a}}(0)=|\mu|^{2}$ and $N_{\text{b}}(0)=\sinh^{2}r$ are the initial number of excitation of modes $\hat{a}$ and $\hat{b}$ respectively. This state is a good approximation for the initial state of realistic systems (i.e., optomechanical systems).

III Decoupling of a non-linear time-dependent optomechanical cross Kerr Hamiltonian

The first aim of this work is to obtains an analytical solution to the decoupled time-evolution operator (4) given our Hamiltonian (II.2). Techniques exist already that can be used to tackle to decoupling of time evolution operators in the way that is suitable for us Bruschi et al. (2013). These are the techniques we will employ here.

III.1 Decoupling algebra of the nonlinear Hamiltonian

Decoupling of the Hamiltonian (II.2) can be done using a choice of the Hermitian operators $\hat{G}_{n}$, see Appendix A.

The first step is to write

[TABLE]

where we have introduced $\hat{\theta}_{2}=\hat{\theta}_{2}(\hat{a}^{\dagger}\,\hat{a}):=\tau+2\,\int_{0}^{\tau}\,d\tau^{\prime}\,\tilde{g}_{2}^{\prime}(\tau^{\prime})\,\hat{a}^{\dagger}\hat{a}$ for notational convenience. Notice that $\hat{\theta}_{2}$ is a time-dependent operator-function of the Hermitian operator $\hat{a}^{\dagger}\,\hat{a}$. We can think of $\hat{U}_{\textrm{I}}(\tau)$ as the time evolution operator of the interaction picture.

We now need to consider the operator \hat{U}_{\textrm{I}}(\tau):=\overset{\leftarrow}{\mathcal{T}}\exp\bigl{[}-\frac{i}{\hbar}\int_{0}^{\tau}\hat{H}_{\textrm{I}}(\tau)\bigr{]} induced by the Hamiltonian $\hat{H}_{\textrm{I}}(\tau)=\hat{H}_{\mathrm{sq}}(\tau)+\hat{H}_{\mathrm{OM}}(\tau)$, where

[TABLE]

We then split the evolution as $\hat{U}_{\textrm{I}}(\tau)=\hat{U}_{\mathrm{sq}}(\tau)\,\hat{U}_{\mathrm{OM}}(\tau)$, where we have defined

[TABLE]

We need to supplement this ansatz with a second one, namely with the action of the two-mode squeezing-like operator $\hat{U}_{\mathrm{sq}}$ on the operators $\hat{b}$ and $\hat{b}^{\dagger}$. We have

[TABLE]

where $\hat{\alpha}=\hat{\alpha}(\hat{a}^{\dagger}\hat{a})$ and $\hat{\beta}=\hat{\beta}(\hat{a}^{\dagger}\hat{a})$ are the Bogoliubov coefficients. These coefficients are time-dependent functions of the operator $\hat{a}^{\dagger}\hat{a}$. They satisfy the Bogoliubov identities $\hat{\alpha}\,\hat{\alpha}^{\dagger}-\hat{\beta}\,\hat{\beta}^{\dagger}=\mathds{1}$ and $\hat{\alpha}\,\hat{\beta}^{\dagger}-\hat{\beta}\,\hat{\alpha}^{\dagger}=0$.

This allows us to obtain

[TABLE]

where we have defined $\hat{E}_{2}=\hat{E}_{2}(\hat{a}^{\dagger}\hat{a}):=\cos(\hat{\theta}_{2})+i\,\sin(\hat{\theta}_{2})$, the coupling $\tilde{g}_{1}(\tau):=\tilde{g}^{(+)}_{1}(\tau)+i\,\tilde{g}_{1}^{(-)}(\tau)$ and $\Re(\hat{O}):=\frac{1}{2}[\hat{O}+\hat{O}^{\dagger}]$ for convenience of notation.

Given all of the above, we can finally obtain

[TABLE]

which is the main expression for the decoupled time-evolution operator in this work.

To complete our main result (III.1) we require the expression of $\hat{F}^{(2)}$ and of the Bogoliubov coefficients $\hat{\alpha}$ and $\hat{\beta}$. The expression for $\hat{F}^{(2)}$ is readily found as

[TABLE]

This leaves us with the task of computing the Bogoliubov coefficients $\hat{\alpha}$ and $\hat{\beta}$ in order to obtain a fully analytical understanding of our system.

III.2 The action of the single-mode squeezing operator

We have noted that, in general, the action of $\hat{U}_{\mathrm{sq}}$ on the operators $\hat{b}$ has the form (9). Ideally, we would like to have an analytical expression for the functional form of $\hat{\alpha}$ and $\hat{\beta}$ in terms of $\hat{a}^{\dagger}\hat{a}$, and of the couplings of the system. Although this is not possible in general Moore and Bruschi (2016), we proceed to construct two uncoupled differential equations that relate the derivatives of $\hat{\alpha}$ and $\hat{\beta}$ with the couplings $\tilde{g}^{\pm}_{2}(\tau)$ and $\tilde{g}_{2}^{\prime}(\tau)$.

The expression for the single mode squeezing is given in (9). In Appendix B we show that

[TABLE]

where we have introduced the modulus $\tilde{\chi}(\tau):=\sqrt{\tilde{g}^{(+)2}_{2}(\tau)+\tilde{g}^{(-)2}_{2}(\tau)}$, the angle $\phi_{2}$ through $\tan(2\,\phi_{2})=\tilde{g}^{(-)}_{2}/\tilde{g}^{(+)}_{2}$ and the operator-function $\hat{p}_{11}=\hat{p}_{11}\left(\int_{0}^{\tau}\,d\tau^{\prime}\,\tilde{\chi}(\tau^{\prime})\right)$, which satisfies the second-order differential equation

[TABLE]

The derivative here is with respect to $y(\tau):=\int_{0}^{\tau}\,d\tau^{\prime}\tilde{\chi}_{2}(\tau^{\prime})$. These differential equations have to be supplemented by the initial conditions $\hat{p}_{11}(0)=\mathds{1}$ and $\dot{\hat{p}}_{11}(0)=0$.

The evolution of the initial state $\hat{\rho}(0)$ is obtained by the usual Heisenberg equation $\hat{\rho}(\tau)=\hat{U}(\tau)\,\hat{\rho}(0)\,\hat{U}^{\dagger}(\tau)$. The full expression is not illuminating and we do not print it here. However, we can ask what is the expression for the time evolution of the mode operators $\hat{a}$ and $\hat{b}$. These expressions allow us, in principle, to compute the expectation value of most quantities of interest.

We define $\hat{a}(\tau):=\hat{U}^{\dagger}(\tau)\,\hat{a}\,\hat{U}(\tau)$ and $\hat{b}(\tau):=\hat{U}^{\dagger}(\tau)\,\hat{b}\,\hat{U}(\tau)$ and, using the expression (III.1), we find

[TABLE]

Here $\hat{I}=\hat{I}(\hat{a}^{\dagger}\,\hat{a}):=\int_{0}^{\tau}\,d\tau^{\prime}\left[\tilde{g}_{1}(\tau^{\prime})\,\hat{E}_{2}\,\hat{\alpha}^{\dagger}+\tilde{g}_{1}^{*}(\tau^{\prime})\,\hat{E}^{\dagger}_{2}\,\hat{\beta}\right]$ for notational convenience. The expression for $\hat{a}(\tau)$ depends on the explicit functional form of $\hat{F}^{(2)}$, \Im\bigl{[}\tilde{g}_{1}(\tau^{\prime})\,\hat{E}_{2}\,(\hat{\alpha}^{\dagger}-\hat{\beta}^{\dagger})\bigr{]} and \Re\bigl{[}\tilde{g}_{1}(\tau^{\prime\prime})\,\hat{E}_{2}\,(\hat{\alpha}^{\dagger}+\hat{\beta}^{\dagger})\bigr{]} in terms of $a^{\dagger}\hat{a}$. This can be computed once the scenario of interest has been determined.

The number $N_{\text{b}}(\tau):=\langle\hat{U}^{\dagger}(\tau)\,\hat{b}^{\dagger}\,\hat{b}\,\hat{U}(\tau)\rangle$ of ‘phonons’ at any time $\tau$ can be computed using (15), and it reads

[TABLE]

In the expression of (III.2) the subscript $n$ means within each operator we need to replace $n\rightarrow\hat{a}^{\dagger}\,\hat{a}$.

IV Applications: time evolution of specific systems with cross-Kerr terms

We are now in the position to study the time evolution induced by the Hamiltonian (III.1) within some specific cross-Kerr scenarios. We will be able to obtain some analytical expressions for meaningful quantities that encode the full nonlinear character of the system.

IV.1 Cross-Kerr without squeezing

Here we start with the scenariowhere $\tilde{g}^{(\pm)}_{2}(\tau)=\tilde{\chi}_{2}(\tau)=0$. This implies that the only cross Kerr term that does not vanish is the term $\tilde{g}_{2}^{\prime}(\tau)\,\hat{a}^{\dagger}\hat{a}\,\hat{b}{}^{\dagger}\hat{b}$ in (II.2). Since $\tilde{\chi}_{2}(\tau)=0$, we can immediately see that $\hat{p}_{11}=\mathds{1}$ and therefore $\hat{\alpha}=\mathds{1}$ and $\hat{\beta}=0$. Some algebra allows us then to find the expression for (III.1) in this case, which reads

[TABLE]

where $\hat{F}^{(2)}$ reads

[TABLE]

Notice that, for $\tilde{g}_{2}^{\prime}(\tau)=\tilde{g}^{(-)}_{1}=0$, we recover the results found in the literature, as expected Bruschi and Xuereb (2018).

The expressions (IV.1) allow us to find

[TABLE]

where we have omitted the expression for $\hat{a}(\tau)$ which can be computed but is not illuminating.

The change $\Delta N_{\text{b}}(\tau):=N_{\text{b}}(\tau)-N_{\text{b}}(0)$ in the number of phonons, given our initial state (5), reads

[TABLE]

where $I_{n}=\int_{0}^{\tau}\,d\tau^{\prime}\tilde{g}_{1}(\tau^{\prime})\,\exp[i\,(\tau+2\,n\,\int_{0}^{\tau}\,d\tau^{\prime}\,\tilde{g}_{2}^{\prime}(\tau^{\prime}))]$ for this case.

Finally, we can also compute the mixedness of the reduced state. The calculations can be found in Appendix D. Lengthy algebra, and the use of a similar approach developed in the literature Bruschi (2018), allow us to find

[TABLE]

where we have defined \Delta_{nn^{\prime}}:=\int_{0}^{\tau}\,d\tau^{\prime}\,\tilde{g}_{1}(\tau^{\prime})\bigl{(}n\,e^{i\,\theta_{2,n}(\tau^{\prime})}-n^{\prime}\,e^{\theta_{2,n^{\prime}}(\tau^{\prime})}\bigr{)} and $\theta_{2,n}(\tau):=\tau+2\,i\,n\,\int_{0}^{\tau}\,d\tau^{\prime\prime}\,\tilde{g}_{2}^{\prime}(\tau^{\prime\prime})$ for convenience of presentation. It is easy to find the result for zero temperature: it is sufficient to set $r_{T}=0$ in (21). Notice that when $\Delta_{nn^{\prime}}=0$, i.e., $\tilde{g}_{1}(\tau)=0$, we obtain $S_{N}(\hat{\rho}_{\textrm{m}}(\tau))=S_{N}(\hat{\rho}_{\textrm{m}}(0))$ as expected.

IV.2 Cross-Kerr with squeezing and without diagonal term

Here we consider the simpler scenario when $\tilde{g}_{2}^{\prime}(\tau)=0$. This implies that the cross Kerr term that vanishes is $\tilde{g}_{2}^{\prime}(\tau)\,\hat{a}^{\dagger}\hat{a}\,\hat{b}{}^{\dagger}\hat{b}$ in (III.1). Since $\tilde{g}_{2}^{\prime}(\tau)=0$, the main differential equation (14) reads

[TABLE]

This case, requires numerical integration of the differential equation above, since there is no analytical solution for a general form of the coupling $\tilde{\chi}_{2}(\tau)$. A numerical approach can then enable the used of (III.1).

IV.3 Cross-Kerr with constant squeezing and constant diagonal term

In this final case we consider $\tilde{g}_{2}^{\prime}(\tau)=\tilde{g}_{2}^{\prime}$ and $\tilde{g}_{2}(\tau)=\tilde{g}_{2}$ constant. This implies that also $\tilde{\chi}_{2}$ and $\phi_{2}$ are constant. In turn, this means that (14) reads

[TABLE]

given that $\hat{\theta}_{2}=(1+2\,\tilde{g}_{2}^{\prime}\,\hat{a}^{\dagger}\hat{a})\,\tau$ in this scenario. This allows us to find the solution

[TABLE]

where we have introduced $\hat{K}:=1+2\,\tilde{g}_{2}^{\prime}\,\hat{a}^{\dagger}\hat{a}$ and $\hat{\Lambda}:=\sqrt{\hat{K}^{2}-4\,\tilde{\chi}_{2}^{2}\,(a^{\dagger}\hat{a})^{2}}$ for convenience of presentation. Note that $\hat{\Lambda}^{-1}\,\sin(\hat{\Lambda}\,\tau)$ is a well defined operator, since it is the short hand notation for the expression \hat{\Lambda}^{-1}\,\sin\bigl{(}\hat{\Lambda}\,\tau\bigr{)}=\bigl{[}1-(\hat{\Lambda}\,\tau)^{2}/3!+(\hat{\Lambda}\,\tau)^{4}/5!-\ldots\bigr{]}\,\tau. Also note that $\hat{\theta}_{2}=\hat{K}\,\tau$ in our present case. Finally, note that if we were to set $\tilde{\chi}_{2}=0$, we would have $\hat{\Lambda}=\hat{K}$ and we can immediately see that $\hat{p}_{11}(\tau)=\mathds{1}$, as expected from our results above.

Recalling that $\hat{p}_{11}(\tau)$ in (24) gives us one of the Bogoliubov coefficients directly, i.e., $\hat{\alpha}=\hat{p}_{11}(\tau)$, we can easily find the other coefficient by employing (III.2), which gives

[TABLE]

It is immediate to check that $\hat{\alpha}\,\hat{\alpha}^{\dagger}-\hat{\beta}\,\hat{\beta}^{\dagger}=\mathds{1}$.

This allows us to obtain an analytical expression for the time evolution operator (III.1). We first obtain $\hat{F}^{(2)}$, which reads (III.1), together with an expression for $\Re\left[\tilde{g}_{1}(\tau)\,\hat{E}_{2}\,(\hat{\alpha}^{\dagger}+\hat{\beta}^{\dagger})\right]$ and $\Im\left[\tilde{g}_{1}(\tau)\,\hat{E}_{2}\,(\hat{\alpha}^{\dagger}-\hat{\beta}^{\dagger})\right]$, whose full expressions can be found in (B.3).

We do not need an analytical expression for the decoupled form of the operator $\hat{U}_{\mathrm{sq}}$ for this case, since we have found the Bogoliubov coefficients (IV.3).

This means that our time evolution operator (III.1) has the explicit expression (B.3), which can be specialised to any desired functional expression of the $\tilde{g}_{1}$ drive.

V Conclusions

We have studied the time evolution induced by a class of Hamiltonians of two interacting quantum harmonic oscillators which include cross-Kerr interactions. We employed tools developed to obtain an analytical expression for the time-evolution operator, which in turn allowed us to compute explicitly quantities of interest. These include the number expectation value of the mechanical resonator and the mixedness of the state of the resonator. Our results are free from approximations and therefore encode the full nonlinear character of the interaction.

These techniques, and the control gained by employing them, can have many applications. For example, similar decoupling techniques have been recently applied the the analysis of correlations within tripartite coupled bosonic system interacting with quadratic Hamiltonians Bruschi et al. (2017), which prompted and facilitated a successful experiment investigating the resulting entanglement between the bosonic modes Lähteenmäki et al. (2016). Among other possible applications of our tools there is quantum control Dong and Petersen (2010), understanding of the interplay between the linear and nonlinear character of quantum mechanical systems Qvarfort et al. (2018) and extending current studies of hidden quantum correlations in existing electromechanical measurements Ockeloen-Korppi et al. (2018).

Finally, it is important to note that our system is ideal and closed, in the sense that we have assumed it to be isolated from the environment, without losses and without decoherence. We leave it to future work to include these important aspects in order to achieve a more realistic and concrete analysis of optomechanical physics driven by cross-Kerr interactions.

Acknowledgments

We would like to thank Anja Metelmann, Sorin Paraoanu, Alessio Serafini and Myung-Joong Hwang for useful comments. D. E. B. would like to thank the University of Vienna and the Erwin Schrödinger Institute (ESI) for hospitality.

Appendix A Decoupling of techniques for time evolution

In this appendix, we outline the general decoupling techniques that we shall be using throughout this work to find a decoupled time-evolution operator generated by the Hamiltonian in (II.2).

A.1 Decoupling for arbitrary Hamiltonians

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

[TABLE]

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 Bruschi et al. (2013) that it is always possible to obtain the decoupling

[TABLE]

where we have defined $\hat{U}_{n}:=\exp[-i\,F_{n}(t)\,\hat{G}_{n}]$ and the real, time-dependent functions $F_{n}(t)$.

The functions $F_{n}(t)$ are uniquely determined by the coupled, nonlinear, first order differential equations

[TABLE]

This is the general method we has been employed in this work, and in previous related one .

Here, we find it convenient to consider the closed finite $9$-dimensional Lie algebra generated by the following set of Hermitian basis operators

[TABLE]

which are the maximal extension of those that generate the Hamiltonian (II.2). Note that there are an infinite of these operators, i.e., there are operators for all $N\in\mathbb{N}$.

A.2 Decoupling for quadratic Hamiltonians

If the Hamiltonian is quadratic in the mode operators the techniques described in the previous subsection have a more powerful representation. Here we proceed to describe these techiques.

A.3 Continuous variables and Covariance Matrix formalism

In quantum mechanics, the initial state $\hat{\rho}_{\textrm{I}}$ of a system of $N$ bosonic modes with operators $\{\hat{a}_{n},\hat{a}^{{\dagger}}_{n}\}$ evolves to a final state $\hat{\rho}_{f}$ through the standard Heisenberg equation $\hat{\rho}_{f}=\hat{U}^{{\dagger}}\,\hat{\rho}_{\textrm{I}}\,\hat{U}$, where $U$ implements the transformation of interest, such as time evolution. If the state $\hat{\rho}$ is Gaussian and the Hamiltonian $H$ is quadratic in the operators, it is convenient to introduce the vector $\hat{\mathbb{X}}=(\hat{a}_{1},\ldots,\hat{a}_{N},\hat{a}^{{\dagger}}_{1},\ldots,\hat{a}^{{\dagger}}_{N})^{Tp}$, the vector of first moments $d:=\langle\hat{\mathbb{X}}\rangle$ and the covariance matrix $\boldsymbol{\sigma}$ defined by $\sigma_{nm}:=\langle\{\hat{X}_{n},\hat{X}^{{\dagger}}_{m}\}\rangle-2\langle\hat{X}_{n}\rangle\langle\hat{X}_{m}^{{\dagger}}\rangle$, where $\{\cdot,\cdot\}$ stands for anticommutator and all expectation values of an operator $\hat{\mathcal{A}}$ are defined by $\langle\hat{\mathcal{A}}\rangle:=\text{Tr}(\hat{\mathcal{A}}\,\hat{\rho})$. In this language, the canonical commutation relations read $[\hat{X}_{n},\hat{X}_{m}^{{\dagger}}]=i\,\Omega_{nm}$, where the $2N\times 2N$ matrix $\boldsymbol{\Omega}$ is known as the symplectic form Adesso et al. (2014). We then notice that, while arbitrary states of bosonic modes are, in general, characterised by an infinite amount of degrees of freedom, a Gaussian state is uniquely determined by its first and second moments, $d_{n}$ and $\sigma_{nm}$ respectively Adesso et al. (2014). Furthermore, quadratic (i.e., linear) unitary transformations, such as Bogoliubov transformations, preserve the Gaussian character of the Gaussian state and can always be represented by a $2N\times 2N$ symplectic matrix $\boldsymbol{S}$ that preserves the symplectic form, i.e., $\boldsymbol{S}^{{\dagger}}\,\boldsymbol{\Omega}\,\boldsymbol{S}=\boldsymbol{\Omega}$.

The Heisenberg equation can be translated in this language to the simple equation $\boldsymbol{\sigma}_{f}=\boldsymbol{S}^{{\dagger}}\,\boldsymbol{\sigma}_{\textrm{I}}\,\boldsymbol{S}$ for the second moments, and $\boldsymbol{r}_{f}=\boldsymbol{S}\,\boldsymbol{r}_{\textrm{I}}$ for the first moments, which shifts the problem of usually untreatable operator algebra to simple $2N\times 2N$ matrix multiplication. Finally, Williamson’s theorem guarantees that any $2N\times 2N$ hermitian matrix, such as the covariance matrix $\boldsymbol{\sigma}$, can be decomposed as $\boldsymbol{\sigma}=\boldsymbol{S}^{{\dagger}}\,\boldsymbol{\nu}_{\oplus}\,\boldsymbol{S}$, where $\boldsymbol{S}$ is an appropriate symplectic matrix, the diagonal matrix $\boldsymbol{\nu}_{\oplus}=\textrm{diag}(\nu_{1},\dots,\nu_{N},\nu_{1},\dots,\nu_{N})$ is known as the Williamson form of the state and $\nu_{n}:=\coth(\frac{2\,\hbar\,\omega_{n}}{k_{B}\,T})\geq 1$ are the symplectic eigenvalues of the state Williamson (1936).

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

The time evolution operator has a symplectic representation $\boldsymbol{S}(t)=\overset{\leftarrow}{\mathcal{T}}\exp[\boldsymbol{\Omega}\,\int_{0}^{\tau}\boldsymbol{H}]$, where $\hat{H}=\frac{1}{2}\,\mathbb{X}^{{\dagger}}\,\boldsymbol{G}_{n}\,\mathbb{X}$, and the decoupled ansatz (27) has the form

[TABLE]

where we have introduced $\boldsymbol{S}_{n}:=\exp[F_{n}(t)\,\boldsymbol{\Omega}\,\boldsymbol{S}_{n}]$, the matrices $\boldsymbol{S}_{n}$ through $\hat{G}_{n}=\frac{1}{2}\,\mathbb{X}^{{\dagger}}\,\boldsymbol{G}_{n}\,\mathbb{X}$, and the real time-dependent functions $F_{n}(t)$ that are the same as they would be obtained with the technique above.

The real, time dependent functions $F_{n}(\tau)$ can be obtained by solving the following system of coupled nonlinear first order differential equations

[TABLE]

This is the matrix version of the operator differential equations (28) for quadratic Hamiltonians, which reduces the problem of operator algebra to matrix multiplication.

Appendix B Decoupling of the squeezing Hamiltonian

Here we employ the techniques of Appendix (A) to decouple the operator $\hat{U}_{\mathrm{sq}}=\overset{\leftarrow}{\mathcal{T}}\exp\left[-\frac{i}{\hbar}\int_{0}^{\tau}\,d\tau^{\prime}\,\hat{H}_{\mathrm{sq}}(\tau^{\prime})\right]$ induced by the squeezing Hamiltonian $\hat{H}_{\mathrm{sq}}(\tau)$, which we reprint here

[TABLE]

B.1 Preliminaries

We start by defining $\mathbb{X}:=(\hat{b},\hat{b}^{\dagger})^{\mathrm{T}}$ an noting that we can write

[TABLE]

where the $2\times 2$ symplectic matrix $\boldsymbol{S}_{\mathrm{sq}}(\tau)$ is the symplectic representation of $\hat{H}_{\mathrm{sq}}(\tau)$ and satisfies $\boldsymbol{S}_{\mathrm{sq}}^{\dagger}(\tau)\,\boldsymbol{\Omega}\,\boldsymbol{S}_{\mathrm{sq}}(\tau)=\boldsymbol{\Omega}$. Here $\boldsymbol{\Omega}=\text{diag}(-i,i)$ is the symplectic form. The matrix $\boldsymbol{S}_{\mathrm{sq}}(\tau)$ therefore has the expression $\boldsymbol{S}_{\mathrm{sq}}(\tau)=\overset{\leftarrow}{\mathcal{T}}\,\exp[\boldsymbol{\Omega}\,\int_{0}^{\tau}\,d\tau^{\prime}\,\boldsymbol{H}_{\mathrm{sq}}(\tau^{\prime})]$. In this case we have

[TABLE]

where we have defined $\tilde{g}_{2}(\tau)=\tilde{\chi}_{2}(\tau)\,\exp[2\,i\,\phi_{2}]:=\tilde{g}^{(+)}_{2}(\tau)+i\,\tilde{g}^{(-)}_{2}(\tau)$, and therefore we have $\tilde{\rho}(\tau)=\sqrt{\tilde{g}^{(+)2}_{2}(\tau)+\tilde{g}^{(-)2}_{2}(\tau)}\geq 0$ and $\tan(2\,\phi_{2})=\tilde{g}^{(-)}_{2}/\tilde{g}^{(+)}_{2}$.

Therefore, through simple algebra we obtain

[TABLE]

B.2 Solving the matrix time-ordered exponential

Here we wish to find a formal expression for (35). We start by noticing that, if we wrote down the time ordered exponential we would be able to write

[TABLE]

where the matrix $\boldsymbol{K}$ is defined as

[TABLE]

and the diagonal matrix $\boldsymbol{P}$ is the object that we need to compute now. Note that it is straightforward to check that $\boldsymbol{P}$ is diagonal.

We use the fact that

[TABLE]

to find the equation

[TABLE]

Since $\boldsymbol{K}$ is invertible (for all cases except when $\tilde{\chi}_{2}(\tau)$=0, which implies $\boldsymbol{P}=\mathds{1}$), we can manipulate this equation and obtain, after some algebra,

[TABLE]

We can now solve the four differential equations contained in (40), two of which are trivial and read $P_{12}=P_{21}=0$, which read

[TABLE]

The differential equations (B.2) do not admit an explicit solution in general. However, they can be integrated numerically when an explicit form of $\tilde{g}_{2}(\tau)$ and $\tilde{g}^{\prime}(\tau)$ are given.

These differential equations have to be supplemented by initial conditions. We note that, since the left hand side of (36) is the identity matrix for $\tau=0$, we have that $\boldsymbol{P}(0)=0$ which implies $\hat{p}_{11}(0)=P_{22}(0)=1$. In addition, taking the time derivative of both sides of (36) and setting $t=0$ it is easy to check that this implies $\dot{P}_{11}(0)=\dot{P}_{22}(0)=0$.

Given that there differential equations are valid only when $\tilde{g}_{2}(\tau)\neq 0$, it is convenient to introduce the functions $\hat{p}_{11}(y)$ and $\hat{p}_{22}(y)$ defined as

[TABLE]

where we have introduced $y(\tau):=2\,\int_{0}^{\tau}\,d\tau^{\prime}\tilde{\chi}_{2}(\tau)$. In terms of these functions we have that the differential equations (B.2) read

[TABLE]

where the derivative stands for derivative with respect to the variable $y$.

This has allowed us to find

[TABLE]

Given that we know that, in general, one has

[TABLE]

this immediately allows us to identify the Bogoliubov coefficients as

[TABLE]

with the auxiliary consistency condition $\hat{p}_{22}=\hat{p}^{\dagger}_{11}$. This condition follows also from (B.2).

B.3 Useful expressions

Here we present a list of useful expressions that are too cumbersome to appear in the main text, but include key steps for the obtainment of the final results.

We start by presenting the expressions for $\Re\left[\tilde{g}_{1}(\tau)\,\hat{E}_{2}\,(\hat{\alpha}^{\dagger}+\hat{\beta}^{\dagger})\right]$ and $\Im\left[\tilde{g}_{1}(\tau)\,\hat{E}_{2}\,(\hat{\alpha}^{\dagger}-\hat{\beta}^{\dagger})\right]$ for the case of constant couplings $\tilde{\chi}_{2}$ and $\tilde{g}^{\prime}_{2}$, which read

[TABLE]

This allows us to find

[TABLE]

Appendix C Properties of the final state

The final state $\hat{\rho}(\tau)$ can be partially computed analytically, given the expressions for our initial state (5) and the time-evolution operator (III.1), and reads

[TABLE]

which can be made slightly more explicit as

[TABLE]

In the expression (C), the subscript $n$ and $m$ mean that the respective quantities are evaluated for $\hat{a}^{\dagger}\hat{a}\rightarrow n,m$.

C.1 Reduced final state of the mechanical oscillator

We can compute the reduced the final reduced state $\hat{\rho}_{\textrm{m}}(\tau)$ of the mechanical resonator. We can be obtain it as $\hat{\rho}_{\textrm{m}}(\tau):=\textrm{Tr}_{\textrm{c}}(\hat{\rho}(\tau))$, which we can compute using (C). We find

[TABLE]

Note that here, again, the subscript $n$ means that we need to replace $n\rightarrow\hat{a}^{\dagger}\,\hat{a}$.

Appendix D Mixedness in the final reduced state of the mechanical oscillator

Here we compute the mixedness of the reduced state of the mechanical oscillator. Given that any pure state $\hat{\rho}$ has the property that $\hat{\rho}^{2}=\hat{\rho}$, it follows that the mixedness of a state can be quantified by the linear entropy $S_{N}$ defined as $S_{N}(\hat{\rho}):=1-\textrm{Tr}(\hat{\rho}^{2})$.

In general, although the expression (C.1) gives us some insight on the reduced state of the mechanical mode, it is difficult to obtain the mixedness for arbitrary parameters of the Hamiltonian. We therefore proceed to specialise to interesting, yet general enough, scenarios below.

D.1 Mixedness of the reduced state of the mechanical oscillator: no cross Kerr squeezing

In the case where $\tilde{g}_{2}=0$ we have $\hat{U}_{\mathrm{sq},n}=\hat{\alpha}=\mathds{1}$ and $\hat{\beta}=0$. Therefore, with some algebra we can show that the expression (C.1) simplifies to

[TABLE]

Therefore, we can easily show that the linear entropy $S_{N}(\hat{\rho}_{\textrm{m}}(\tau))$ associated to this state reads

[TABLE]

where we have introduced

[TABLE]

and defined \Delta_{nn^{\prime}}:=\int_{0}^{\tau}\,d\tau^{\prime}\,\tilde{g}_{1}(\tau^{\prime})\bigl{(}n\,e^{i\,\tau^{\prime}+2\,i\,n\,\int_{0}^{\tau^{\prime}}\,d\tau^{\prime\prime}\,\tilde{g}_{2}^{\prime}(\tau^{\prime\prime})}-n^{\prime}\,e^{i\,\tau^{\prime}+2\,i\,n^{\prime}\,\int_{0}^{\tau^{\prime}}\,d\tau^{\prime\prime}\,\tilde{g}_{2}^{\prime}(\tau^{\prime\prime})}\bigr{)}. Note that the exact form of $e^{i\,\psi_{nn^{\prime}}}$ does not matter since it cancels out in (53).

It is also easy to check that $S_{N}(\hat{\rho}_{\textrm{m}}(\tau))$ can be also written as

[TABLE]

Now, we can follow the same procedure introduced in Bruschi (2018) and compute (D.1) explicitly. We find

[TABLE]

We note that, as expected, the result (56) reduces to the existing one in the literature for $\tilde{g}_{2}^{\prime}(\tau)=0$, see Bruschi (2018).