Higher-order interactions in quantum optomechanics: Analytical solution of nonlinearity

TL;DR

This paper introduces a perturbation method for solving nonlinear quantum Langevin equations involving quadratic interactions, providing analytical spectral densities and correlation functions, with applications to quantum systems like lasers and amplifiers.

Contribution

It proposes first and second order truncation schemes using higher-order operators for analytical solutions in nonlinear quantum optomechanics.

Findings

Derived spectral densities of higher-order operators.

Found an analytical expression for the second-order correlation function.

Predicted cavity photon occupation at laser threshold as √6−2.

Abstract

A method is described to solve the nonlinear Langevin equations arising from quadratic interactions in quantum mechanics. While, the zeroth order linearization approximation to the operators is normally used, here first and second order truncation perturbation schemes are proposed. These schemes employ higher-order system operators, and then approximate number operators with their corresponding mean boson numbers, only where needed. Spectral densities of higher-order operators are derived, and an expression for the second-order correlation function at zero time-delay has been found, which reveals that the cavity photon occupation of an ideal laser at threshold reaches $\sqrt{6}-2$, in good agreement with extensive numerical calculations. As further applications, analysis of the quantum anharmonic oscillator, calculation of $Q-$functions, analysis of quantum limited amplifiers, and nondemoliton measurements.

Figures

Equations

$$\mathbb{H}=\hslash \gamma (\hat{b}\pm {\hat{b}}^{†}{)}^2(\hat{a}\pm {\hat{a}}^{†}{)}^2,$$

$$\mathbb{H}=\hslash \gamma {\hat{a}}^{†}\hat{a}(\hat{b}\pm {\hat{b}}^{†}{)}^2,$$

$$\hat{c}$$

$$\hat{d}$$

$$[\hat{c},{\hat{c}}^{†}]$$

$$[\hat{c},\hat{n}]$$

$$[{\hat{c}}^{†},\hat{n}]$$

$$[\hat{c},{\hat{a}}^{†}]$$

$$[\hat{d},{\hat{d}}^{†}]$$

$$[\hat{d},\hat{m}]$$

$$[{\hat{d}}^{†},\hat{m}]$$

$$[\hat{d},{\hat{b}}^{†}]$$

$$\frac{d}{dt}\{A\}=[\text{M}]\{A\}-\sqrt{[\Gamma ]}\{A_{\mathrm{in}}\},$$

$$\{A_{\mathrm{out}}\}=\{A_{\mathrm{in}}\}+\sqrt{[\Gamma ]}\{A\}.$$

$$\{A_{\mathrm{out}}\}=[\text{S}]\{A_{\mathrm{in}}\}.$$

$$[\text{S}(w)]=[\text{I}]-\sqrt{[\Gamma ]}{\left(iw[\text{I}]+[\text{M}]\right)}^{-1}\sqrt{[\Gamma ]}.$$

$$\dot{\hat{z}}=\frac{d}{dt}\hat{z}=-\frac{i}{\hslash }[\hat{z},\mathbb{H}]-[\hat{z},{\hat{x}}^{†}](\frac{1}{2}\Gamma \hat{x}+\sqrt{\Gamma }{\hat{z}}_{\mathrm{in}})+(\frac{1}{2}\Gamma {\hat{x}}^{†}+\sqrt{\Gamma }{\hat{z}}_{\mathrm{in}}^{†})[\hat{z},\hat{x}],$$

$$\mathbb{H}=2\hslash \gamma \hat{n}(\hat{d}+{\hat{d}}^{†}+\hat{m}).$$

$$\dot{\hat{a}}$$

$$\dot{\hat{d}}$$

$${\dot{\hat{d}}}^{†}$$

$$\dot{\hat{m}}$$

$$\dot{\hat{a}}$$

$$\dot{\hat{d}}$$

$${\dot{\hat{d}}}^{†}$$

$$\dot{\hat{m}}$$

$$\text{raggedleft}\frac{d}{dt}\left\{\begin{array}{c}\hat{a}\\ \hat{d}\\ {\hat{d}}^{†}\\ \hat{m}\end{array}\right\}\text{@add@raggedleft}$$

$$\mathbb{H}=4\hslash \gamma (\hat{d}+{\hat{d}}^{†}+\hat{m})(\hat{c}+{\hat{c}}^{†}+\hat{n}),$$

$$\text{raggedleft}\dot{\hat{c}}\text{@add@raggedleft}$$

$${\dot{\hat{c}}}^{†}$$

$$\dot{\hat{n}}$$

$$\dot{\hat{d}}$$

$${\dot{\hat{d}}}^{†}$$

$$\dot{\hat{m}}$$

$$\dot{\hat{c}}$$

$${\dot{\hat{c}}}^{†}$$

$$\dot{\hat{n}}$$

$$\dot{\hat{d}}$$

$${\dot{\hat{d}}}^{†}$$

$$\dot{\hat{m}}$$

$$\left[\begin{array}{cccccc}-i8\gamma \bar{m}-\frac{2\bar{n}+1}{4}{\Gamma }_1& 0& -i4\gamma \bar{m}& -i2\gamma (2\bar{n}+1)& -i2\gamma (2\bar{n}+1)& -i2\gamma (2\bar{n}+1)\\ 0& i8\gamma \bar{m}-\frac{2\bar{n}+1}{4}{\Gamma }_1& i4\gamma \bar{m}& i2\gamma (2\bar{n}+1)& i2\gamma (2\bar{n}+1)& i2\gamma (2\bar{n}+1)\\ i8\gamma \bar{m}& -i8\gamma \bar{m}& 0& 0& 0& 0\\ -i2\gamma (2\bar{m}+1)& -i2\gamma (2\bar{m}+1)& -i2\gamma (2\bar{m}+1)& -i8\gamma \bar{n}-\frac{2\bar{m}+1}{4}{\Gamma }_2& 0& -4i\gamma \bar{n}\\ i2\gamma (2\bar{m}+1)& i2\gamma (2\bar{m}+1)& i2\gamma (2\bar{m}+1)& 0& i8\gamma \bar{n}-\frac{2\bar{m}+1}{4}{\Gamma }_2& i4\gamma \bar{n}\\ 0& 0& 0& i8\gamma \bar{n}& -i8\gamma \bar{n}& 0\end{array}\right]$$

$$\times \left\{\begin{array}{c}\hat{c}\\ {\hat{c}}^{†}\\ \hat{n}\\ \hat{d}\\ {\hat{d}}^{†}\\ \hat{m}\end{array}\right\}-\left\{\begin{array}{c}\sqrt{{\Delta }_1}{\hat{c}}_{\text{in}}\\ \sqrt{{\Delta }_1}{\hat{c}}_{\text{in}}^{†}\\ 0\\ \sqrt{{\Delta }_2}{\hat{d}}_{\text{in}}\\ \sqrt{{\Delta }_2}{\hat{d}}_{\text{in}}^{†}\\ 0\end{array}\right\}=\frac{d}{dt}\left\{\begin{array}{c}\hat{c}\\ {\hat{c}}^{†}\\ \hat{n}\\ \hat{d}\\ {\hat{d}}^{†}\\ \hat{m}\end{array}\right\},$$

Full text

I Introduction

In quantum optomechanics the standard interaction Hamiltonian is simply the product of photon number $\hat{n}=\hat{a}^{\dagger}\hat{a}$ and the position $x_{\textrm{zp}}(\hat{b}+\hat{b}^{\dagger})$ operators 1 ; 2 ; 2a ; 3 ; 4 ; 5 , where $x_{\textrm{zp}}$ is the zero-point motion, and $\hat{a}$ and $\hat{b}$ are respectively the photon and phonon annihilators. This type of interaction can successfully describe a vast range of phenomena, including optomechanical arrays 6 ; 7 ; 8 ; 9 ; 10 ; 11 ; 12 , squeezing of phonon states 14 ; 14a ; 14b , non-reciprocal optomechanics 15 ; 16 ; 17 ; 17a , Heisenberg’s limited measurements 17b , sensing 18 ; 19 ; 20 , engineered dissipation and states 21 ; 22 , and non-reciprocal acousto-optics 23 . In all these applications, the mathematical toolbox to estimate the measured spectrum is Langevin equations 24 ; 25 ; 26 ; 26a .

Usually, the analysis of quantum optomechanics is done within the linearized approximation of photon ladder operators, normally done as $\hat{a}\rightarrow\bar{a}+\delta\hat{a}$ with $|\bar{a}|^{2}=\bar{n}$ being the mean cavity photon number, while nonlinear terms in $\delta\hat{a}$ are ignored. But this suffers from limited accuracy wherever the basic optomechanical interaction $\mathbb{H}_{\rm OM}=\hbar g_{0}\hat{n}(\hat{b}+\hat{b}^{\dagger})$ is either vanishingly small or non-existent. In fact, the single-photon interaction rate $g_{0}$ can be identically made zero by appropriate design 27 ; 28 ; 29 ; 30 , when quadratic or even quartic effects are primarily pursued. This urges need for accurate knowledge of higher-order interaction terms.

Some other optomechanical phenomena such as four-wave mixing, also can be suitably understood by incorporation of higher-order interaction terms 31 . Recent experiments 32 ; 33 have already established the significance and prominent role of such type of nonlinear interactions. In fact, quadratic nonlinear optomechanics 33.0 ; 33.1 ; 33.2 ; 33.3 ; 33.4 ; 33.5 ; 33.6 ; 33.6.0 ; 33.6.1 ; 33.6.2 ; 33.6.3 ; 33.6.4 ; 33.6.5 ; 33.6.5a ; 33.6.6 ; 33.6.7 is now a well recognized subject of study even down to the single-photon level 33.7 , for which circuit analogues have been constructed 33.8.0 ; 33.8 and may be regarded as fairly convenient simulators 33.8a ; 33.8b ; 33.8c of much more complicated experimental optomechanical analogues. Dual formalisms of quadratic optomechanics are also found in ultracold atom traps 33.9 ; 33.9a as well as optical levitation 33.10 . Such types of nonlinear interactions also appear elsewhere in anharmonic quantum circuits 33.10.0 . Quadratic interactions are in particular important for energy and non-demolition measurements of mechanical states 1 ; 2 ; 3 ; 33.11 ; 33.12 ; 33.13 . While the simple linearization of operators could be still good enough to explain some of the observations, there remains a need for an exact and relatively simple mathematical treatment. Method of Langevin equations also normally fails, and other known methods such as expansion unto number states and master equation, require lots of computation while giving little insight to the problem.

Perturbative expansions and higher-order operators have been used by other researchers to study noise spectra of lasers 33.15 ; 33.16 ; 33.19 ; 33.20 . Also, the master equation approach 33a ; 33b can be used in combination with the quasi-probablity Wigner functions 33c ; 33d to yield integrable classical Langevin equations. Nevertheless, a method recently has been proposed 34 , which offers a truncation correlation scheme for solution of driven-dissipative multi-mode systems. While being general, it deals with the time evolution of expectation values instead of operators within the truncation accuracy, so the corresponding Langevin equations cannot be analytically integrated.

Alternatively, a first-order perturbation has been proposed to tackle the nonlinear quadratic optomechanics 34a . This method perturbatively expands the unknown parameters of classical Langevin equations for the nonlinear system, and proceeds to the truncation at first order. However, the expansion is accurate only where the ratio of photon loss rate to mechanical frequency $\kappa/\Omega$ is large. This condition is strongly violated for instance in superconductive electromechanical systems.

So far and to the best knowledge of author, no treatment of quadratic interactions using Langevin equations for operators has been reported. This paper presents a perturbative mathematical treatment within the first and second order approximations to the nonlinear system of Langevin equations, which ultimately result in an integrable system of quantum mechanical operators. The trick here is to introduce operators of higher dimensionality into the solution space of the problem. Having their commutators calculated, it would be possible to set up an extended system of Langevin equations which could be conveniently solved by truncation at the desirable order. To understand how it works, one may consider the infamous first order quadratic nonlinear Riccati differential equation 34b ; 34c , which is exactly integrable if appropriately transformed as a system of two coupled linear first order differential equations. Alternatively, Riccati equation could be exactly transformed into a linear second order differential equation, too. But this is not what we consider here, since it will result in a much more complicated second-order system of Langevin equations involving derivatives of noise terms.

The method introduced here is useful in other areas of quantum physics 33.9 ; 33.10 than optomechanics, where nonlinearities such as anharmonic or Kerr interactions are involved. We also describe how the $Q-$functions could be obtained for the anharmonic oscillator. Further applications of nonlinear stochastic differential equations 34.8 ; 34.9 ; 34.10 beyond stochastic optomechanics 33.6.5 ; 33.6.5a includes finance and stock-market analysis 34.11 , turbulence 34.12 ; 34.12a , hydrology and flood prediction 34.13 , and solar energy 34.14 . Also, the Fokker-Planck equation 33.20 ; 34.15 ; 34.16 ; 34.16a ; 34.16b is actually equivalent to the nonlinear Schrödinger equation with bosonic operator algebra, and its moments 34.17 translate into nonlinear Langevin equations. Similarly, this method can deal with side-band generation in optomechanics 35 , superconducting circuits 35a , as well as spontaneous emission in open systems 36 ; 36.0 . Applications in estimation of other parameters such as the second order correlation $g^{(2)}(0)$ 36a ; 36b ; 36c ; 36d , quantum limited amplifiers 36e ; 36f and quantum nondemolition measurements 36f ; 36g ; 36h ; 36i are demonstrated, and furthermore it is found that an unsqueezed ideal laser reaches $\sqrt{6}-2$ cavity photons at threshold.

II Theory

II.1 Hamiltonian

A nonlinear quadratic optomechanical interaction in the most general form 37 is here defined as

[TABLE]

where $\gamma$ is the interaction rate. Furthermore, bosonic photon $\hat{a}$ and phonon $\hat{b}$ ladder operators satisfy $[\hat{b},\hat{b}^{\dagger}]=[\hat{a},\hat{a}^{\dagger}]=1$ as well as $[\hat{b},\hat{a}]=[\hat{b},\hat{a}^{\dagger}]=0$. Meanwhile, quadratic interactions normally are 1 ; 2 ; 3

[TABLE]

which by defining the photon number operator $\hat{n}=\hat{a}^{\dagger}\hat{a}$ takes essentially the same algebraic form.

Direct expansion of (1) shows that it essentially brings in a different interaction type compared to (2). Doing so, we obtain $\mathbb{H}=\hbar\gamma(\hat{b}^{2}+\hat{b}^{\dagger 2}\pm 2\hat{m}\pm 1)^{2}(\hat{a}^{2}+\hat{a}^{\dagger 2}\pm 2\hat{n}\pm 1)$ where $\hat{m}=\hat{b}^{\dagger}\hat{b}$. Hence, (1) includes interactions of type $\hat{a}^{2}\hat{b}^{2}$, $\hat{a}^{2}\hat{b}^{\dagger 2}$, and so on, which are absent in (2). It should be noticed that the widely used standard optomechanical interaction $\mathbb{H}_{\rm OM}$ results in nonlinear and linear Langevin equations when expressed respectively in the terms of $\{\hat{a},\hat{b}\}$ and $\{\hat{n},\hat{x}\}$. Hence, this type of interaction is not addressed here. In addition to the above Hamiltonians (1,2), there exist still other types of nonlinear optomechanical interactions 14b ; 38 such as $\mathbb{H}=\hbar g(\hat{b}\pm\hat{b}^{\dagger})(\hat{a}^{2}\pm\hat{a}^{\dagger 2})$, which is also not considered explicitly here, but can be well treated using the scheme presented in this article.

II.2 Linear Perturbation

This approach is being mostly used by authors to solve the systems based on either (1) or (2). To this end, ladder field operators are replaced with their perturbations, while product terms beyond are neglected and truncated. Obviously, this will give rise to interactions of the type $\hbar(\hat{b}\pm\hat{b}^{\dagger})^{2}(q\delta\hat{a}+q^{\ast}\delta\hat{a}^{\dagger})$, where $q=2\gamma(\bar{a}\pm\bar{a}^{\ast})$ for (1) and $q=\gamma\bar{a}$ for (2) is some complex constant in general, and $\delta\hat{a}$ now represents the perturbation term around the steady state average $|\bar{a}|=\sqrt{\bar{n}}$. This technique is mostly being referred to as the linearization of operators, and directly leads to an integrable set of Langevin equations if also applied to the mechanical displacement as well.

II.3 Square Field Operators

Here, we define the square field operators 37

[TABLE]

for photons, which obviously satisfy $[\hat{c},\hat{a}]=[\hat{c},\hat{b}]=[\hat{d},\hat{a}]=[\hat{d},\hat{b}]=[\hat{c},\hat{d}]=0$. Now, it is not difficult to verify that these operators furthermore satisfy the commutation relationships

[TABLE]

Defining the phonon number operator as $\hat{m}=\hat{b}^{\dagger}\hat{b}$, in a similar manner we could write

[TABLE]

The set of commutator equations (4) and (5) enables us to treat the quadratic nonlinear interaction perturbatively to the desirable accuracy, as is described in the following.

II.4 Langevin Equations

The input/output formalism 24 ; 25 ; 26 ; 26a can be used to assign decay channels to each of the quantum variables of the system. This will result in the set of Langevin equations

[TABLE]

where $\{A\}$ is the system vector, $[\textbf{M}]$ is the coefficients matrix whose eigenvalues need to have negative or vanishing real parts to guarantee stability, and $[\Gamma]$ is a real-valued matrix which is diagonal if all noise terms corresponding to the members of $\{A\}$ are mutually independent. When $[\textbf{M}]$ is independent of $\{A\}$, (6) is linear and integrable and otherwise nonlinear and non-integrable. If $[\textbf{M}(t)]$ is a function of time, then (6) is said to be time-dependent. Furthermore, $\{A_{\rm in}\}$ represents the input fields to the system at the respective ports, and $\{A_{\rm out}\}$ is the output fields, which are related together as 4 ; 5 ; 6

[TABLE]

Here, $[\Gamma]$ is supposed to be diagonal for simplicity. From the scattering matrix formalism we also have

[TABLE]

Hence, taking $w$ as the angular frequency and performing a Fourier transform on (6), the scattering matrix is found by using (7) and (8) as

[TABLE]

Hence, $[\textbf{S}]$ is well-defined if $[\textbf{M}]$ is known. This can be obtained by using the Langevin equations

[TABLE]

where $\hat{x}$ is any system operator, which is here taken to be the same as $\hat{z}$ to comply with (8).

By setting either $\hat{z}=\hat{c}$ or $\hat{z}=\hat{d}$ the commutators in (10) by (4) or (5) always lead back to the same linear combination of these forms. Thus, the new set of Langevin equations is actually linear in terms of the square or higher-order operators, if perturbatively truncated at a finite order. So, instead of solving the nonlinear system in linearized $2\times 2$ space $\{A\}^{\rm T}=\{\hat{a},\hat{b}\}$, one may employ an expanded dimensional space with increased accuracy. There, truncation and sometimes mean field approximations are necessary to restrict the dimension, since commutators of new operators mostly lead to even higher-orders and are thus not closed under commutation. As examples, a $4\times 4$ space $\{A\}^{\rm T}=\{\hat{a},\hat{d},\hat{d}^{\dagger},\hat{m}\}$ truncated at the first-order, or a $6\times 6$ space $\{A\}^{\rm T}=\{\hat{c},\hat{c}^{\dagger},\hat{n},\hat{d},\hat{d}^{\dagger},\hat{m}\}$ truncated at the second-order could be used for (1,2). To illustrate the application of this method, we describe two examples in the next section. It could be extended to the accuracy of the second-order perturbation too, by defining appropriate cross product operator terms between photonic and phononic partitions.

III Examples

Here, we describe two examples from the nonlinear interactions of having type (1) or (2).

III.1 Standard Quadratic Interaction (2)

Analysis of such systems requires analysis in a 4-dimensional space, spanned by $\{A\}^{\rm T}=\{\hat{a},\hat{d},\hat{d}^{\dagger},\hat{m}\}$. Taking the plus sign here without loss of generality and after dropping a trivial non-interacting term $\mathbb{H}_{0}=\hbar\gamma\hat{n}$, the nonlinear interaction is

[TABLE]

This can be found by expansion of (2), plugging in (3) and $[\hat{b},\hat{b}^{\dagger}]=1$, and dropping a trivial term $\hbar\gamma\hat{n}$. Using (5), $[\hat{a},\hat{n}]=\hat{a}$ and $[\hat{a}^{\dagger},\hat{n}]=-\hat{a}^{\dagger}$ in the non-rotating frame of operators, and ignoring the self-energy Hamiltonian $\mathbb{H}_{\rm self}=\hbar(\omega+\gamma)\hat{n}+\hbar\Omega\hat{m}$ for the moment, Langevin equations become

[TABLE]

So far, the set of equations (12) is exact. However, integration of (12) is still not possible at this stage, and taking Fourier transformation must be done later when arriving at a linear operator system. We present a first-order and second-order perturbative method to deal with this difficulty.

It should be furthermore noticed that using a non-rotating frame with the self-energy Hamiltonian $\mathbb{H}_{\rm self}$ not ignored, would have resulted in identical equations, except with the addition of the trivial terms $-i\Delta\hat{a}$, $-i2\Omega\hat{d}$, and $+i2\Omega\hat{d}^{\dagger}$ respectively to the first three equations, where $\Delta=\omega+\gamma-\nu$ is the optical detuning with $\nu$ being the cavity optical resonance frequency, and $\omega$ and $\Omega$ are respectively the optical and mechanical frequencies. Also, the damping coefficient in high mechanical quality factor $Q_{\rm m}$ limit could be estimated as $\Gamma_{2}=2\Gamma_{\rm m}$, where $\Gamma_{\rm m}$ is the damping rate of the $\hat{b}$ phononic field. Here, it is preferable not to use the rotating frames since the coefficients matrix $[\textbf{M}]$ becomes time-dependent.

III.1.1 First-order Perturbation to (12)

Now, if the photon and phonon baths each have a mean boson number respectively as $\left<\hat{n}\right>=\bar{n}$ and $\left<\hat{m}\right>=\bar{m}$, we could immediately write down the linear system of equations in the non-rotating frame of operators and neglection of self-energies $\mathbb{H}_{\rm self}$ as

[TABLE]

which is now exactly integrable. Here, we use the linearization $2\hat{a}\hat{d}=(\bar{a}+\delta\hat{a})\hat{d}+\hat{a}(\bar{d}+\delta\hat{d})\rightarrow\bar{a}\hat{d}+\bar{d}\hat{a}$, where $\bar{d}=\frac{1}{2}\bar{a}^{2}$ and higher-order terms of the form $\delta\hat{a}\delta\hat{d}$ are dropped, and so on. But this cannot be applied to $\hat{n}\hat{m}=\hat{a}^{\dagger}\hat{a}\hat{m}$ since $\hat{n}$ and $\hat{a}^{\dagger}$ are absent from the basis. Furthermore, any linearization of this expansion would generate terms $\hat{a}\hat{m}$ and $\hat{a}^{\dagger}\hat{m}$ which are still nonlinear. Both of these issues can be resolved by a second-order perturbation as follows next. This results in the operator equations

[TABLE]

where $\sqrt{\Delta_{1}}=\sqrt{\Gamma_{1}}$ and $\sqrt{\Delta_{2}}=\left(\bar{m}+\frac{1}{2}\right)\sqrt{\Gamma_{2}}$. The set of equations (14) is linear and can be easily addressed by standard methods of stochastic Langevin equations used in optomechanics 1 ; 2 ; 3 ; 24 ; 25 and elsewhere. More specifically, one may employ analytical Fourier methods in frequency domain as an matrix algebraic problem to obtain spectra of variables, or integrate the system numerically by stochastic numerical methods in time domain to obtain time dependent behavior of expectation values.

All that remains is to find the average cavity boson numbers for photons $\bar{n}$ and phonons $\bar{m}$. In order to do this, one may first arbitrate $d/dt=0$ in (13) at steady state, and then use the equality of real parts in first equation to find the expression for $\bar{n}$. Doing this, results in $\bar{n}=4|\bar{a}_{\rm in}|^{2}/\Gamma_{1}$ where $|\bar{a}_{\rm in}|$ represents the amplitude of coherent laser input. Also, the initial cavity phonon occupation number at $t=0$ could be estimated simply as $\bar{m}=1/\left[\exp(\hbar\Omega/k_{\rm B}T)-1\right]$ 24 ; 25 , where $k_{\rm B}T$ is the thermal energy with $k_{\rm B}$ and $T$ being respectively the Boltzmann’s constant and absolute temperature. Detailed numerical examinations reveal that the system of equations (14) is generally very well stable with $\Re\{{\rm eig}[{\bf M}]\}<0$ at sufficiently low optical intensities.

III.2 Full Quadratic Interaction (1)

Analysis of a fully quadratic system requires analysis in a $6\times 6$ dimensional space, spanned by $\{A\}^{\rm T}=\{\hat{c},\hat{c}^{\dagger},\hat{n},\hat{d},\hat{d}^{\dagger},\hat{m}\}$. Taking both of the plus signs here, the Hamiltonian could be written as

[TABLE]

where a trivial non-interacting term $\mathbb{H}_{0}=2\hbar\gamma(1+\hat{n}+\hat{m}+\hat{d}+\hat{c}+\hat{d}^{\dagger}+\hat{c}^{\dagger})$ is dropped. The set of Langevin equations can be obtained in a similar manner, and in non-rotating frame of operators with neglection of self-energies $\mathbb{H}_{\rm self}=\hbar(\omega+2\gamma)\hat{n}+\hbar(\Omega+2\gamma)\hat{m}$ for the moment, results in

[TABLE]

Similar to (12), the damping rate for sufficiently high optical quality factors $Q$ could be estimated as $\Gamma_{1}=2\kappa$, where $\kappa$ is the damping rate of the $\hat{a}$ photonic field.

Quite clearly, should we have not ignored the self-energy Hamiltonian $\mathbb{H}_{\rm self}$, then addition of the diagonal terms $-i2\Delta\hat{c}$, $+i2\Delta\hat{c}^{\dagger}$ to the first two where $\Delta=\omega+2\gamma-\nu$ with $\nu$ being the optical cavity resonance frequency, and similarly $-i2\Omega\hat{d}$ and $+i2\Omega\hat{d}^{\dagger}$ to the fourth and fifth equations would have been necessary. These are not shown here only for the sake of convenience. Again, it is emphasized that transformation to the rotating frame of operators here would make the coefficients time-dependent in an oscillating manner, and it is far better to be avoided for these classes of nonlinear problems.

III.2.1 First-order Perturbation to (16)

In a similar manner to (13), we may assume photon and phonon baths each have a mean boson number respectively as $\left<\hat{n}\right>=\bar{n}$ and $\left<\hat{m}\right>=\bar{m}$, which gives

[TABLE]

We here need to assume the redefinition $\sqrt{\Delta_{1}}=(\bar{n}+\frac{1}{2})\sqrt{\Gamma_{1}}$. Now, without taking $\mathbb{H}_{\rm self}$ into account, this will lead to the linear system of matrix Langevin equations

[TABLE]

which is, of course, integrable now. The initial cavity boson numbers $\bar{n}$ and $\bar{m}$ can be set in the same manner which was done for the system of equations (14). Numerical tests reveal that (18) is conditionally stable if the optical intensity is kept below a certain limit on the red detuning, and is otherwise unstable.

III.3 Second order Perturbation to (12,16)

The set of Langevin equations (12,16) can be integrated with much more accuracy, if we first identify and sort out the cross terms as individual operators. For instance, (16) contains the cross operators $\hat{c}\hat{d}$, $\hat{c}\hat{d}^{\dagger}$, $\hat{c}\hat{m}$, $\hat{c}^{\dagger}\hat{d}$, $\hat{c}^{\dagger}\hat{d}^{\dagger}$, $\hat{c}^{\dagger}\hat{m}$, $\hat{n}\hat{d}$, $\hat{n}\hat{d}^{\dagger}$, as well as $\hat{n}\hat{m}$ which is self-adjoint. These constitute an extra set of nine cross operators to be included in the treatment. All these cross operators are formed by multiplication of photonic and phononic single operators, whose notation order, such as $\hat{c}\hat{d}=\hat{d}\hat{c}$ and so on, is obviously immaterial.

Now, one may proceed first to determine the commutators between these terms where relevant, which always result in linear combinations of the other existing terms. This will clearly enable a more accurate formulation of (16) but in a $6+9=15$ dimensional space, which is given by the array of operators $\{A\}^{\rm T}=\{\hat{c},\hat{c}^{\dagger},\hat{n},\hat{d},\hat{d}^{\dagger},\hat{m},\hat{c}\hat{d},\hat{c}\hat{d}^{\dagger},\hat{c}\hat{m},\hat{c}^{\dagger}\hat{d},\hat{c}^{\dagger}\hat{d}^{\dagger},\hat{c}^{\dagger}\hat{m},\hat{n}\hat{d},\hat{n}\hat{d}^{\dagger},\hat{n}\hat{m}\}$.

The independent non-trivial quadratic commutator equations among cross operators here are found after tedious but straightforward algebra as

[TABLE]

The rest of commutators among cross operators are either adjoints of the above, or have a common term which makes their evaluation possible using either (4) or (5). Commutators among cross operators and single operators can be always factored, such as $[\hat{c}\hat{d},\hat{n}]=[\hat{c},\hat{n}]\hat{d}$. Commutators among single operators are already known (4,5). It can be therefore seen that commutators (19) always lead to operators of higher orders yet, so that they do not terminate at any finite order of interest by merely expansion of operators basis. This fact puts the perturbative method put into work. There are, however, nonlinear systems such as semiconductor optical cavities 33.19 ; 33b in which higher-order operators yield an exact closed algebra and satisfy a closedness property within the original space by appropriate definition.

The set of ten commutators now can be perturbatively linearized as a second-order approximation, by replacing the number operators with their mean values, wherever needed to reduce the set of operators back to the available 15 dimensional space. This will give rise to the similar set of equations after some algebra

[TABLE]

where the reduction of triple operator products among single and cross operators as $4\hat{x}\hat{y}\hat{z}\rightarrow\bar{x}\hat{y}\hat{z}+\bar{x}\bar{y}\hat{z}+\bar{y}\bar{z}\hat{x}+\bar{z}\bar{x}\hat{y}$ is used where appropriate. For instance, the term $4\hat{n}\hat{m}^{2}$ is replaced as $\bar{m}\hat{m}\hat{n}+2\bar{m}\bar{n}\hat{m}+\bar{m}^{2}\hat{n}$ and so on. Also, similar to (13), products among single operators are reduced as $2\hat{x}\hat{y}\rightarrow\bar{x}\hat{y}+\bar{y}\hat{x}$. This is somewhat comparable to the mean field approach in cross Kerr optomechanics 38a .

There are two basic reasons why we have adopted this particular approach to the linearization and cuting off the diverging operators of higher orders. The first reason is that number operators vary slowly in time as opposed to their bosonic counterparts which oscillate rapidly in time, given the fact that the use of rotating frames is disallowed here. Secondly, number operators are both positive-definite and self-adjoint, and thus can be approximated by a positive real number. These properties makes the replacements $\hat{n}\rightarrow\bar{n}$ and $\hat{m}\rightarrow\bar{m}$ reasonable approximations, and the replacement with mean values needs only to be restricted to the number operators, to yield a closed algebra necessary for construction of Langevin equations. Hence, the correct application of replacements only to the triple operator products appearing in (19) will make sure that no operator having an order beyond than that of cross operators will appear in the formulation.

Anyhow, it can be seen now that all approximate commutators in (20) allow the set of operators $\{A\}^{\rm T}\cup\{\hat{1}\}=\{\hat{1},\hat{c},\hat{c}^{\dagger},\hat{n},\hat{d},\hat{d}^{\dagger},\hat{m},\hat{c}\hat{d},\hat{c}\hat{d}^{\dagger},\hat{c}\hat{m},\hat{c}^{\dagger}\hat{d},\hat{c}^{\dagger}\hat{d}^{\dagger},\hat{c}^{\dagger}\hat{m},\hat{n}\hat{d},\hat{n}\hat{d}^{\dagger},\hat{n}\hat{m}\}$ to take on linear combinations of its members among every pair of commutations possible, where $\hat{1}$ is the identity operator. Obviously, this approximate closedness property now makes the full construction of Langevin equations for the operators belonging to $\{A\}$ possible. It is noted that $\hat{1}$ is not an identity element for the commutation.

We can now define the set $\{S\}={\rm span}(\{A\}\cup\{\hat{1}\})$, which is spanned by all possible linear combinations of $\{\hat{1}\}$ and the members of $\{A\}$ together with the associative binary commutation operation $[]$ defined in (4,5,20). The ordered pair $(\{{S}\},[])$ is now a semigroup.

Having therefore these ten commutators (20) known, we may proceed now to composing the second-order approximation to the nonlinear Langevin equations (16), from which a much more accurate solution could be obtained. Here, the corresponding Langevin equations may be constructed at each step by setting both $\hat{z}$ and $\hat{x}$ in (10) equal to either of the 15 operators, while the noise input terms for cross operators is a simple product of related individual noise terms. The linear damping rates of higher-order operators is furthermore simply the sum of individual damping rates of corresponding single operators, which completes the needed parameter set of Langevin equations.

IV Further Considerations

IV.1 Optomechanical Interaction & Drive Terms

The method described in the above can be simultaneously used if other terms such as the standard optomechanical interaction $\mathbb{H}_{\textrm{OM}}$ is non-zero, or there exists a coherent pumping drive term which can be expressed as $\mathbb{H}_{\textrm{d}}=\sum_{k}F_{k}\hat{b}^{\dagger}+F_{k}^{*}\hat{b}$, where $F_{k}$ are time-dependent drive amplitudes. While $\mathbb{H}_{\textrm{d}}$ does not appear directly in the Langevin equations, treatment of $\mathbb{H}_{\textrm{OM}}$ requires inclusion of additional Langevin equations for $\hat{a}$ and $\hat{b}$ where appropriate, as well as few extra terms in the rest. This can be done in a pretty standard way, and is not repeated here for the sake of brevity 1 ; 2 ; 24 ; 25 ; 26 ; 26a .

IV.2 Multi-mode Fields

The analysis is also essentially unaltered if there are more than one mechanical mode to be considered 17 ; 39 ; 39a , and the method is still easily applicable with no fundamental change. Suppose that there are a total of $M$ mechanical modes with the corresponding bosonic operators $\hat{b}_{k}$ and $\hat{b}_{k}^{\dagger}$ where $k\in[1,M]$. Then, these modes are mutually independent in the sense that $[\hat{b}_{j},\hat{b}_{k}]=0$ and $[\hat{b}_{j},\hat{b}_{k}^{\dagger}]=\delta_{jk}$. The set of commutators (5) will be usable for all $M$ modes individually and as a result (19) and therefore (20) may be still used. The first and second order perturbations will respectively result in $3+3M=3(M+1)$ and $3+3M+9M=3(4M+1)$ equations. The redefined set of operators will be respectively now $\{A\}^{\rm T}=\{\hat{c},\hat{c}^{\dagger},\hat{n},\hat{d}_{k},\hat{d}_{k}^{\dagger},\hat{m}_{k};k\in[1,M]\}$ and $\{A\}^{\rm T}=\{\hat{c},\hat{c}^{\dagger},\hat{n},\hat{d}_{k},\hat{d}_{k}^{\dagger},\hat{m}_{k},\hat{c}\hat{d}_{k},\hat{c}\hat{d}_{k}^{\dagger},\hat{c}\hat{m}_{k},\hat{c}^{\dagger}\hat{d}_{k},\hat{c}^{\dagger}\hat{d}_{k}^{\dagger},\hat{c}^{\dagger}\hat{m}_{k},\hat{n}\hat{d}_{k},\hat{n}\hat{d}_{k}^{\dagger},\hat{n}\hat{m}_{k};k\in[1,M]\}$.

Similarly, in case of $N$ optical modes satisfying $[\hat{a}_{j},\hat{a}_{k}]=0$ and $[\hat{a}_{j},\hat{a}_{k}^{\dagger}]=\delta_{jk}$, the set of commutators (4) can be used and the operator set should be now expanded as $\{A\}^{\rm T}=\{\hat{c}_{j},\hat{c}_{j}^{\dagger},\hat{n}_{j},\hat{d}_{k},\hat{d}_{k}^{\dagger},\hat{m}_{k};j\in[1,N];k\in[1,M]\}$ and $\{A\}^{\rm T}=\{\hat{c}_{j},\hat{c}_{j}^{\dagger},\hat{n}_{j},\hat{d}_{k},\hat{d}_{k}^{\dagger},\hat{m}_{k},\hat{c}_{j}\hat{d}_{k},\hat{c}_{j}\hat{d}_{k}^{\dagger},\hat{c}_{j}\hat{m}_{k},\hat{c}_{j}^{\dagger}\hat{d}_{k},\hat{c}_{j}^{\dagger}\hat{d}_{k}^{\dagger},\hat{c}_{j}^{\dagger}\hat{m}_{k},\hat{n}_{j}\hat{d}_{k},\hat{n}_{j}\hat{d}_{k}^{\dagger},\hat{n}_{j}\hat{m}_{k};j\in[1,N];k\in[1,M]\}$ respectively for first and second order perturbations. Hence, the corresponding dimensions will be now respectively either $3(N+M)$ or $3(N+M+3NM)$. Higher-order commutators (19) and (20) can be still used again by only addition of appropriate photonic $j$ and phononic $k$ mode indices to the respective operators contained in the expanded operator basis set $\{A\}$.

IV.3 Noise Spectra

The required noise spectra 39b of cross operators is clearly a product of each of the individual terms, since the nature of particles are different. However, the noise spectra of quadratic operators themselves need to be appropriately expressed. For instance, $\hat{d}_{\textrm{in}}$ actually corresponds to the spectral input noise of the square operator $\hat{d}=\hat{b}\hat{b}/2\sqrt{\Gamma}$ from (3), which clearly satisfies $\hat{d}_{\textrm{in}}(t)=\frac{1}{2}\hat{b}_{\textrm{in}}(t)\hat{b}_{\textrm{in}}(t)/\sqrt{\Gamma}$, or $\hat{d}_{\textrm{in}}(w)=\frac{1}{2}\hat{b}_{\textrm{in}}(w)\ast\hat{b}_{\textrm{in}}(w)/\sqrt{\Gamma}$ in the frequency domain, where $\ast$ merely represents the convolution operation. Therefore, once $\hat{a}_{\textrm{in}}(w)$ and $\hat{b}_{\textrm{in}}(w)$ are known, all relevant remaining input noise spectra could be obtained accordingly using simple convolutions or products in frequency domain.

As a result, the corresponding spectral density of the noise input terms to the cross operators can be determined from the relevant vacuum noise fluctuations and performing a Fourier transform. For instance, we have $S_{CDCD}[w]=S_{CC}[w]S_{DD}[w]$ where $S_{CC}[w]=\frac{1}{4}S_{A^{2}A^{2}}[w]$ and $S_{DD}[w]=\frac{1}{4}S_{B^{2}B^{2}}[w]$. Then Isserlis-Wick theorem 32 ; 40 could be exploited to yield the desired expressions. If we assume

[TABLE]

where the dimensionless correlation integrator runs on phase, instead of time, as

[TABLE]

then the functions $\zeta(\cdot),$ $\psi(\cdot),$ and the operator $\hat{v}(\cdot)$ should be all having the dimension of $\hat{f}^{2}(\cdot)$ as well. That means if $\hat{f}$ is dimensionless, which is the case for the choice of ladder operators, then $\zeta(\cdot),$ $\psi(\cdot),$ and $\hat{v}(\cdot)$ become dimensionless, too. The functions $\zeta(\cdot)$ and $\psi(\cdot)$ together can cause squeezing or thermal states if appropriately defined 24 ; 26a . By Isserlis-Wick theorem applied to scalars we have $\left<x_{1}x_{2}x_{3}x_{4}\right>=\left<x_{1}x_{2}\right>\left<x_{3}x_{4}\right>+\left<x_{1}x_{3}\right>\left<x_{2}x_{4}\right>+\left<x_{1}x_{4}\right>\left<x_{3}x_{4}\right>$. This gives for the operators

[TABLE]

Hence, for a given stochastic process where $\left<\hat{f}(t)\hat{f}(\tau)\right>=0$, $\left<\hat{f}(t)\hat{f}^{\dagger}(\tau)\right>=\Psi(t-\tau)$, and having the scalar commutator $[\hat{f}(t),\hat{f}^{\dagger}(\tau)]=\Upsilon(t-\tau)$, we simply get

[TABLE]

Now, suppose that we have a coherent field of photons at the angular frequency $\omega$ with an initial Gaussian distribution, in which $\Psi(t)=\exp(-\chi^{2}\omega^{2}t^{2}/2)\exp(-i\omega t)$ and $\Upsilon(t)=\Psi(t)$, while having the linewidth $\Delta f=\frac{1}{2\pi}\chi\omega$. Clearly, $\chi$ is a dimensionless and positive real number. In the limit of $\chi\rightarrow 0^{+}$, the expected relationship $\Psi(t)=\sqrt{2\pi}\delta(\omega t)/\chi$ is easily recovered.

This particular definition of the correlating function $\Psi(t)$ ensures that the corresponding spectral density is appropriately normalized, that is

[TABLE]

Hence, one may obtain the following spectral density

[TABLE]

which is centered at the doubled frequency $2\omega$, has a linewidth of $\sqrt{2}\Delta f$, and satisfies the property

[TABLE]

Once the spectral densities of input noise terms are found, spectral densities of all output fields immediately follows (8,9) as $\{A[w]\}_{\rm out}=[\textbf{S}^{\dagger}(w)\textbf{S}(w)]\{A[w]\}_{\rm in}$, in which $[\textbf{S}^{\dagger}(w)\textbf{S}(w)]=[|S_{ij}(w)|^{2}]$, $\{A[w]\}_{\rm in}$ is an array containing the spectral densities of inputs, and similarly $\{A[w]\}_{\rm out}$ is the array of spectral densities at each of the output fields.

IV.4 Estimation of $g^{(2)}(0)$

Many of the important features of an interacting quantum system is given by its second-order correlation function $g^{(2)}(0)$ at zero time-delay 36a ; 36b ; 36c defined as

[TABLE]

It is fairly easy to estimate this function once the spectral densities of all higher order operators of the nonlinear system are calculated. For this purpose, we may first employ the definition (3) to rewrite

[TABLE]

Estimation of the average within brackets can be done by having $S_{CC}[w]=\frac{1}{4}S_{A^{2}A^{2}}[w]$ corresponding to the higher-order operator $\hat{c}$. This can be assumed to has been already found from knowledge of the scattering matrix $[\textbf{S}(w)]$, spectral densities of input fields $\{A[w]\}_{\rm in}$, and subsequent derivation of spectral density array of output fields $\{A[w]\}_{\rm out}$. Then, $S_{CC}[w]$ will be simply an element of the vector $\{A[w]\}_{\rm out}$. Using (24), this results in a fairly brief representation

[TABLE]

With the assumptions above for an ideal initial Gaussian distribution, we have $\Psi(0)=\Upsilon(0)=1$ and thus $g^{(2)}(0)=4(\frac{1}{2}-\bar{n})/\bar{n}^{2}$. One should have in mind that this relationship cannot be readily used for a coherent radiation, since for a practical laser the true statistics is Poissonian and not Gaussian. This analysis thus reveals that the cavity occupation number of such an ideal laser with the threshold defined as $g^{(2)}(0)=1$ is exactly $\bar{n}=\sqrt{6}-2\approx 0.450$. This is in contrast to the widely used assumption of quantum threshold condition $\bar{n}=1$ 40a ; 40b ; 40b1 ; 40b2 ; 40b3 ; 40b4 . Interestingly, a new study 40c of photon statistics in weakly nonlinear optical cavities based on extensive density matrix calculations 40d ; 40e yields the value $\bar{n}=0.4172$, which is in reasonable agreement to our estimate. An earlier investigation on quantum-dot photonic crystal cavity lasers 40f ; 40g also gives the value $\bar{n}=0.485$.

V Anharmonic Oscillator

The quantum anharmonic oscillator appears in many nonlinear systems including quadratic optomechanics 40h ; 40i , where our method here is applicable. The anharmonic Kerr Hamiltonian is 41 ; 42

[TABLE]

in which $\zeta$ is a constant. It is well known that in case of $\zeta>2\omega$ this system exhibits an effective bistable potential, and is otherwise monostable. However, we are here much interested in a slightly different but more complicated form given by 2DM

[TABLE]

which is monostable or bistable if both $\omega$ and $\zeta$ are respectively positive or negative. This type of nonlinearity is of particular importance in fourth-order analysis of qubits 43 ; 44 ; 45 ; 46 ; 47 ; 49 ; 50 . While the Hamiltonian (32) is for a single-mode field, the case of multi-mode electromagnetic field could be easily devised following the existing interaction Hamiltonians 2DM and the presented method in this article. Nevertheless, the above expression after some algebraic manipulations can be put into the form

[TABLE]

where a trivial constant term $\hbar\zeta$ is dropped. Here, we may proceed with the 8-dimensional basis operator set $\{A\}^{\rm T}=\{\hat{c},\hat{c}^{\dagger},\hat{n},\hat{n}^{2},\hat{c}^{2},\hat{c}^{\dagger 2},\hat{n}\hat{c},\hat{c}^{\dagger}\hat{n}\}$, resulting in a second order perturbation accuracy.

Treating this problem using the Langevin equation (10), regardless of the values of $\zeta$ and $\omega$, is possible, only if the following non-trivial exact commutators

[TABLE]

are known, which may be found after significant algebra. The rest of required commutators which are not conjugates of those in the above, can either directly or after factorization of a common term be easily found from (4). Again, the set of commutators (34) does not yet satisfy the closedness property within $\{S\}={\rm span}(\{A\}\cup\{\hat{1}\})$, unless the approximate linearization

[TABLE]

is employed. The rest of the process is identical to the one described under (20). Construction of the respective noise terms is also possible by iterated use of the results in §IV.3 and so on.

V.1 The Husimi-Kano Q-functions

It is mostly appropriate that moments of operators are known, which are scalar functions and much easier to work with. The particular choice of $Q-$functions 50a is preferred when dealing with ladder operators, and are obtained by taking the expectation value of density operator with respect to a complex coherent state $\ket{\alpha}$ and dividing by $\pi$. This definition leads to a non-negative real valued function $Q(\alpha)=Q(\Re[\alpha],\Im[\alpha])$ of $\ket{\alpha}$. Then, obtaining $Q-$function moments of any expression containing the ladder operators would be straightforward 50a . However, it must be antinormally ordered, with creators be moved to the right. In $\{A\}^{\rm T}$ above all operators are actually in the normal form, except $\hat{n}^{2}$. It is possible to put the nontrivial members of $\{A\}$ in the antinormal order

[TABLE]

While evaluating $Q-$function moments, $\hat{a}$ and $\hat{a}^{\dagger}$ are replaced with $\alpha$ and $\alpha^{\ast}$ respectively as

[TABLE]

All remains now is to redefine the array of $Q-$functions bases, using common terms as $\{\braket{A}\}^{\rm T}=\{\alpha^{2},\alpha^{\ast 2},|\alpha|^{2},|\alpha|^{4},\alpha^{4}$, $\alpha^{\ast 4},\alpha^{2}|\alpha|^{2},\alpha^{\ast 2}|\alpha|^{2}\}$ from which the original $Q-$functions could be readily restored. This translates into a set of scalar differential equations which conveniently could be solved. Fluctuations of noise terms also vanish while taking the expectation values, and only their average values survive. To illustrate this, suppose that the system is driven by a coherent field $\hat{a}_{\rm in}$ with the normalized electric field amplitude $\beta=\alpha/\sqrt{2}$ and at the frequency $\omega$. Then, the $Q-$function moments of the input fields after defining the loss rates $\Gamma_{3}=2\Gamma_{2}=4\Gamma_{1}$ become $\braket{\hat{a}_{\rm in}}=\sqrt{2\Gamma_{1}}\beta$, $\braket{\hat{c}_{\rm in}}=\sqrt{\Gamma_{2}}\beta$, $\braket{\hat{n}_{\rm in}}=\sqrt{\Gamma_{2}}(2|\beta|^{2}+1)$, $\braket{\hat{c}^{2}_{\rm in}}=\sqrt{\Gamma_{3}}\beta^{2}$, and $\braket{\hat{n}_{\rm in}\hat{c}_{\rm in}}=\sqrt{\Gamma_{3}}\beta^{2}(2|\beta|^{2}+3)$.

V.2 Quantum Limited Amplifiers

The same method can be extended to the quantum limited amplifiers, which in the general form coincides with the expression (31), but is usually solved using a zeroth-order perturbation 36e . For the single-mode degenerate quantum limited amplifier 5 ; 36e ; 36e1 , the corresponding Hamiltonian is slightly different given by $\mathbb{H}=\hbar\omega\hat{n}+\hbar(g\hat{c}+g^{\ast}\hat{c}^{\dagger})$, with the 3-dimensional basis $\{A\}^{\rm T}=\{\hat{n},\hat{c},\hat{c}^{\dagger}\}$ which satisfies closedness. Then, the second-order accurate Langevin equations with inclusion of the self-energy $\mathbb{H}_{\rm self}=\hbar\omega\hat{n}$ can be shown to be unconditionally stable with $\Re\{{\rm eig}[{\bf M}]\}<0$, given by

[TABLE]

In presence of Kerr nonlinearity 36f as $\mathbb{H}=\hbar\omega\hat{n}+\hbar(g\hat{c}+g^{\ast}\hat{c}^{\dagger})+\hbar\gamma\hat{c}^{\dagger}\hat{c}$, one may use $4\hat{c}^{\dagger}\hat{c}=\hat{n}^{2}-\hat{n}$, $[\hat{n}^{2},\hat{c}]\approx-\frac{1}{2}(6\bar{n}+7)\hat{c}$, and the basis $\{A\}^{\rm T}=\{\hat{n},\hat{n}^{2},\hat{c},\hat{c}^{\dagger}\}$ to construct a set of $4\times 4$ integrable Langevin equations. The rest of necessary commutators are already found in (4), (34), and (35).

V.3 Quantum Nondemolition Measurements

Quantum nondemolition measurements of states require a cross-Kerr nonlinear interaction of the type $\mathbb{H}=\hbar\omega\hat{a}^{\dagger}\hat{a}+\hbar\Omega\hat{b}^{\dagger}\hat{b}+\hbar\chi\hat{a}^{\dagger}\hat{a}\hat{b}^{\dagger}\hat{b}=\hbar\omega\hat{n}+\hbar\Omega\hat{m}+\hbar\chi\hat{n}\hat{m}$, in which $\hat{a}$ and $\hat{b}$ fields respectively correspond to the probe and signal 36g ; 36h . This system can be conveniently analyzed by the preferred choice 36g of the higher-order operators $\{A\}^{\rm T}=\{\hat{n},\hat{m},\hat{C},\hat{S}\}$, where

[TABLE]

are quadratures of the readout observable. It is straightforward to show by induction that $[f(\hat{a}^{\dagger}),\hat{a}]=-f^{\prime}(\hat{a}^{\dagger})$ and $[\hat{a}^{\dagger},f(\hat{a})]=-f^{\prime}(\hat{a})$ with $f(\cdot):\mathcal{R}\mapsto\mathcal{R}$ being a real function of its argument. Now, the non-zero commutators of the basis $\{A\}^{\rm T}$ can be found after some algebra as $[\hat{n},\hat{C}]=-i\hat{S}$, $[\hat{n},\hat{S}]=i\hat{C}$, and $[\hat{C},\hat{S}]=\frac{1}{2}i(\hat{n}+2)^{-1}$. All remains to construct the Langevin equations now, is to linearize the last commutators as $[\hat{C},\hat{S}]\approx\frac{1}{2}i(\bar{n}+2)^{-1}$, by which the basis $\{A\}^{\rm T}=\{\hat{n},\hat{m},\hat{C},\hat{S}\}$ would satisfy closedness. Input noise terms to the operators $\hat{C}$ and $\hat{S}$ should be constructed by linear combinations of $\hat{a}_{\rm in}$ and $\hat{a}_{\rm in}^{\dagger}$ while replacing the multiplier term $1/\sqrt{\hat{n}+1}$ with the linearized form $1/\sqrt{\bar{n}+1}$.

In the end, it has to be mentioned that under external drive, periodicity, or dynamical control $[\textbf{M}(t)]$ in (6) is time-dependent 36f ; 51 . For instance, the ultimate optomechanical cooling limit is a function of system dynamics 56 . Then, integration should be done numerically, since exact analytical solutions without infinite perturbations exist only for very restricted cases. This is, however, beyond the scope of the current study.

VI Conclusions

A new method was described to solve quadratic quantum interactions using perturbative truncation schemes, by including higher-order operators in the solution space. Spectral densities of higher-order operators, calculation of the second-order correlation function, as well as the quantum anharmonic oscillator and transformation to scalar forms using $Q-$functions were discussed. Finally, applications of the presented approach to quantum limited amplifiers, and nondemolition measurements were demonstrated.

Acknowledgements.

This work been supported by Laboratory of Photonics and Quantum Measurements at École Polytechnique Fédérale de Lausanne and Research Deputy of Sharif University of Technology. The author thanks Prof. Franco Nori, Prof. Vincenzo Savona, Dr. Alexey Feofanov, Dr. Christophe Galland, Dr. Sahar Sahebdivan, as well as Liu Qiu and Amir H. Ghadimi for comments and/or discussions. The author is highly indebted to Dr. Hiwa Mahmoudi at Institute of Electrodynamics, Microwave and Circuit Engineering in Technische Universität Wien, and in particular, the Laboratory for Quantum Foundations and Quantum Information on the Nano- and Microscale in Vienna Center for Quantum Science and Technology (VCQ) at Universität Wien for their warm and receptive hospitality during which the numerical computations and final revisions took place. The huge effort needed in improving the presentation of this article has not been possible without support and encouragement of the celebrated artist, Anastasia Huppmann.

Appendices

We show by numerical solution of a stochastic nonlinear operator differential equation, and also taking its expectation values within the Mean Field Approximation, that the proposed analytical scheme in the above referenced manuscript works very well, is convergent, and uniformly converges to the accurate solution.

We also demonstrate the second-order exact solution to the standard optomechanical Hamiltonian, and show that the proposed method can reproduce the side-band asymmetry in quantum optomechanics surprisingly well, which is a rigorous proof of the quantum mechanical capability of the proposed approach.

Appendix A Nonlinear First-order Circuit

First, consider the infinitely-ordered nonlinear operator equation

[TABLE]

which models the voltage operator of an RC circuit shunted by a nonlinear ideal diode, driven by a sinusoidal voltage source $v\left(t\right)=V_{0}e^{-\alpha t}\sin\left(\omega t\right)$, and stochastic noise $\hat{n}\left(t\right)$. We suppose that the noise $\hat{n}\left(t\right)$ is governed by a Weiner process. Here, and without loss of generality, both $\kappa$ and $\mu$ are taken to be positive real parameters. Hence, this model does not include an oscillating part due to an imaginary $\mu$, which could have been otherwise absorbed into $\hat{u}\left(t\right)$ by a rotating frame transformation. This particular choice also eliminates the imaginary part of $\hat{u}\left(t\right)$. We also assume here, for the illustrative purpose of this example, that $\mu=V_{0}=\tau=1$, $\omega=2\pi\alpha,$ and $\omega=2\pi\times\text{1kHz}$.

The reason for choosing this particular differential equation is that it is nonlinear to the infinite order, and also the extended basis of higher order operators all commute and therefore trivially form a closed basis.

Using the proposed method in the paper under consideration, this above operator equation can be first put into the infinitely-ordered linear system of ordinary differential equations as

[TABLE]

Subsequently, the input terms can be linearized using the proposed method in the paper. Doing this results in

[TABLE]

in which $\bar{u}=\frac{1}{T}\int^{T}_{0}{\left\langle\hat{u}\left(t\right)\right\rangle dt}$ is the time-average of the input. Now, the above system of equations can be exactly integrated, after truncation to a finite-order.

Using an extensive code written in Mathematica, the above system of linear stochastic equations can be treated and integrated as an Itô process, and the results for various orders of truncation between 2 and 6 versus the numerically exact solution are displayed in Fig. 1.

It is still not quite clear that the method is convergent to the exact solution, since the Itô integration of a Weiner process every time is carried over a different sequence of random numbers. This difficulty cannot be avoided in principle, since there is no way to reset the numerically random sequence.

Therefore, as a double check, we take the expectation values, which discards the noise term, and transform a mean-field approximation to reach a similar system of differential equations, however, expressed in terms of the expectation value function $\left\langle\hat{u}\left(t\right)\right\rangle$ and its higher orders. This is equivalent to solving the nonlinear differential equation

[TABLE]

given the fact that $\left\langle\hat{n}\left(t\right)\right\rangle=0$.

Doing this immediately reveals the convergence property of our proposed method, illustrated in Fig. 2. As it can be clearly verified, the numerical solutions are so rapidly and accurately converging to the exact solution, that they are practically indistinguishable beyond the two lowest truncation orders.

Appendix B Optomechanical Hamiltonian

The standard optomechanical Hamiltonian is not quadratic, but can be solved within the second-order accuracy using the method described in the article. Here, we are able to show that the proposed approach actually can reproduce fundamental quantum mechanical properties of a typical nonlinear system, such as side-band asymmetry.

Its standard optomechanical Hamiltonian reads

[TABLE]

In order to form a closed basis of operators, we may choose

[TABLE]

which forms a $6\times 6$ system of Langevin equations. It is easy to verify that this system is exactly closed, by calculation of all possible commutation pairs between the elements. Out of the $6!$ commutators, the non-zero ones are

[TABLE]

This is obviously a closed basis. Now, one may proceed with composition of the Langevin equations. They are given by

[TABLE]

where $\gamma=\kappa+\Gamma$, we have set $\hat{x}=\hat{a}$ in (10) in all equations, ${\hat{n}}_{\text{in}}={\hat{a}}^{\dagger}{\hat{a}}_{\text{in}}+\hat{a}{{\hat{a}}^{\dagger}}_{\text{in}}$ and ${\hat{c}}_{\text{in}}=2\hat{a}{\hat{a}}_{\text{in}}$. This subsequently can be linearized to get the second-order accurate optomechanical system of equations as

[TABLE]

in which $F=g_{0}\sqrt{\bar{n}}$, $K=4\bar{n}\kappa$, $G=g_{0}\bar{n}$, $L=G+g_{0}\left(\bar{m}+1\right)$, and can be further approximated by $L=G$ under normal experimental conditions of an ultracold cavity. The average population value $\bar{m}={1}/{\left[\exp\left({\hbar\Omega}/{k_{\text{B}}T}\right)\right]}$, while $\bar{n}$ can be obtained from the steady state solution of the first row by replacements of input noise term $\sqrt{\kappa}{\hat{a}}_{\text{in}}\to\alpha$, where $\alpha$ is the input photon flux.

Doing this gives the third order equation in terms of $\sqrt{\bar{n}}$ as

[TABLE]

Solution of this third-order algebraic equation gives the three distinct solutions

[TABLE]

in which

[TABLE]

is a complex number. These three roots clearly become two distinct solutions in the limit of lossless cavity with $\gamma=0$ as long as $27g_{0}C^{2}>4B_{r}^{2}$, and very well reproduce the expected bistable behavior of cavity photon number $\bar{n}$ with respect to various parameters such as detuning $\Delta$ and input photon flux $\alpha$. This has been shown for typical normalized variables in Figs. 3 and 4.

For the noise terms, we may insert

[TABLE]

Here, the spectral density of ${{\hat{a}}_{\text{in}}}^{2}$ has already been calculated in Section IV.3. As opposed to the second-order accurate optomechanical Langevin equations, the first-order accurate equations are simply obtained by truncating (150) which recovers the well-known $3\times 3$ system with the obviously closed basis $A^{\text{T}}=\left\{\hat{a},\hat{b},{\hat{b}}^{\dagger}\right\}$ after simple algebraic manipulation

[TABLE]

Here, we have used the further replacements $g_{0}\hat{n}\to G$ and $g_{0}\hat{a}\hat{b}\to F\hat{b}$. However, if we had made the replacement $g_{0}\hat{n}\to F\hat{a}$ we would get

[TABLE]

Finally, we are able to show that this proposed approach can reproduce the side-band asymmetry of an optomechanical cavity, which is a fundamental quantum mechanical property.

In the context of quantum optomechanics, it is general practice to solve either (198) or (199) with asymmetric noise spectral densities such as $S_{BB}(+\omega)\neq S_{BB}(-\omega)$. While being a purely nonlinear quantum mechanical effect, this is the only known way at this moment to reproduce the side-band asymmetry. Here, we can show indeed that asymmetric spectral noise terms are not needed to obtain the much expected side-band asymmetry.

This result becomes possible by considering that in the system (183) both of the terms $\hat{a}\hat{b}$ and $\hat{a}\hat{b}^{\dagger}$ are present in the second-order accurate basis. The first term and its conjugate $\{\hat{a}\hat{b},\hat{a}^{\dagger}\hat{b}^{\dagger}\}$ represent the parametric interactions on the blue side, while the second term and its conjugate $\{\hat{a}\hat{b}^{\dagger},\hat{a}\hat{b}^{\dagger}\}$ represent the hopping interactions on the red side. By assumption of some test values to the input parameters as well as symmetric input noise terms such as $S_{BB}(+\omega)=S_{BB}(-\omega)$, the spectral density of the collected output light from the cavity can be calculated.

This has been done for the two cases of pumping on the red-side and blue-side shown in Fig. 5. More interestingly and as the ultimate verification of the approach, the pumping can be tried exactly on the cavity resonance. Carrying out this procedure allows calculation of cavity sidebands at blue and red, the ratio of which is here plotted and shown in Fig. 6. This figure reveals that there actually exists an easily observable asymmetry, which is hallmark of nonlinear quantum mechanical effects in cavity optomechanics. This should have rigorously proven, while leaving no further doubt, that the second-order accurate system of equations (183) can reproduce the side-band asymmetry in a satisfactory manner.

Detailed numerical simulation of standard and non-standard optomechanical problems using the proposed approach in paper, and verification against experiments remains as a subject of a future study.