1. Introduction
Optomechanical systems (OMS) have raised a rapidly growing interest, since it not only produces intriguing phenomena such as mode conversion [1–3], optomechanically induced transparency [4, 5] and optomechanical entanglement [6–9], but have extensive applications ranging from optical feedback cooling [10–12], noise squeezing of the light beam [13–15] to quantum information processing [16–18]. Recently, the unprecedented experimental progress like mechanical squeezing control [19], single-phonon level quantum control [20–23] and fabrication of entanglement between light and mechanics [7] has greatly promoted the development in OMS, such as the increasing coupling strengths between light and mechanical motion in optomechanical setups, which have improved by orders of magnitude [24–26]. There are also several theoretical studies on OMS in strong-coupling regimes, for example, the optical nonlinearities like photon blockade [27–30] and nonlinear dynamics [31] in OMS have been studied, and another study found that the steady state mechanical Wigner density contains strong negative parts at strong optomechanical coupling, signaling stable nonclassical states [32].
One of the most fundamental quantum properties of emitting light is the statistics, which can be characterized by correlation functions [33, 34] measured via a Hanbury–Brown–Twiss experiment [35]. Similarly, the emitting phonons can also exhibit some quantum statistical properties, e.g. phonon anti-bunching behaviors and phonon blockade effects [36, 37], which can also be described by the phonon correlation functions. The statistics of photons in OMS have been extensively studied [27, 38, 39], however, to our best knowledge, there have been few studies on that of phonons. For example, the phonon sub-Poissonian statistics and phonon blockade effects in quadratically coupled OMS have been studied by solving the master equation numerically [40–43] or analytically via the phase space method [44]. Although the phonon statistical properties are explored, a comprehensive understanding of different scattering processes is required to reveal the mechanism of phonon (anti-)bunching behaviors.
In this paper, we study the quantum statistics of phonon radiations in linearly coupled OMS at low temperatures, with the aim to establish the relation of the phonon (anti-)bunching behavior on relevant system parameters. In addition to providing the numerical results based on the standard master equation approach, we also calculate the second-order correlation functions analytically using the scattering formalism in the weak driving and coupling limit. With the analytical result, we can identify three dominant scattering processes, from which an unambiguous mechanism of the phonon pair generation naturally arises. Both the numerical and analytical results show that the anti-bunching behavior of two phonons requires a large detuning between the frequencies of the incident laser and the cavity.
The paper is organized as follows. In section 2 , we illustrate OMS with the Hamiltonian, and derive the formal expression of correlation functions. In section 3 , we evaluate correlation functions numerically using the master equation, which agree with that from perturbative expansion with high accuracy. Inspired by the correspondence between the scattering theory and the master equation, in section 4 , we also use the field theoretical approach to obtain the correlation functions analytically, which allows us to find out the condition of phonon (anti-)bunching behaviors and provide an intuitive picture of phonon pair generations. Finally, the results are summarized with an outlook in section 5 .
2. Model
We consider the OMS consist of a photon mode with frequency ωa and a phonon mode with frequency ωb as shown in figure 1. The photon and phonon modes decay into two independent reservoirs respectively. The Hamiltonian is H = Hc + Hr + Hcr + Hd, where the cavity is described by
$\begin{eqnarray}{H}_{c}={\omega }_{a}{a}^{\dagger }a+{\omega }_{b}{b}^{\dagger }b+{{ga}}^{\dagger }a(b+{b}^{\dagger }),\end{eqnarray}$
where a†(a) is the creation (annihilation) operator of the photon mode, b†(b) is that of the phonon mode and g is the single-photon coupling strength. The propagation of photons and phonons in reservoirs is described by the Hamiltonian $\begin{eqnarray}{H}_{r}=\int {\rm{d}}{kk}{\alpha }_{k}^{\dagger }{\alpha }_{k}+\int {\rm{d}}{qq}{\beta }_{q}^{\dagger }{\beta }_{q},\end{eqnarray}$
where ${\alpha }_{k}^{\dagger }({\alpha }_{k})$ is the creation (annihilation) operator of photons in the reservoir and ${\beta }_{q}^{\dagger }({\beta }_{q})$ is that of phonons in the reservoir. The decay of photons and phonons into reservoirs is depicted by $\begin{eqnarray}\begin{array}{rcl}{H}_{\mathrm{cr}} & = & {\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{a}}{2\pi }}\displaystyle \int {\rm{d}}k({\alpha }_{k}^{\dagger }a-{a}^{\dagger }{\alpha }_{k})\\ & & +{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{b}}{2\pi }}\displaystyle \int {\rm{d}}q({\beta }_{q}^{\dagger }b-{b}^{\dagger }{\beta }_{q})\end{array}\end{eqnarray}$
with the decay rates γa(γb) of photons (phonons). The cavity is driven by a laser field, which is characterized by $\begin{eqnarray}{H}_{d}={\rm{\Omega }}({a}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\omega }_{d}t}+{\rm{h}}.{\rm{c}}.).\end{eqnarray}$
Figure 1. Schematic of OMS in which the optical cavity contains a photon mode, whose creation (annihilation) operator is a†(a), coupling to a phonon mode, whose creation (annihilation) operator is b†(b). The two cavity modes decay into the reservoirs respectively. |
In the rotation frame with respect to ωd, where $\rho \to {{\rm{e}}}^{{\rm{i}}{\omega }_{d}{a}^{\dagger }{at}}\rho {{\rm{e}}}^{-{\rm{i}}{\omega }_{d}{a}^{\dagger }{at}}$, ωa → Δ = ωa − ωd, Hd → ω(a† + a), by eliminating the reservoir degrees of freedom, we obtain the master equation [45]
$\begin{eqnarray}{\dot{\rho }}_{c}=-{\rm{i}}\left[{H}_{c}+{H}_{d},{\rho }_{c}\right]+{D}_{a}({\rho }_{c})+{D}_{b}({\rho }_{c})\equiv { \mathcal L }[{\rho }_{c}].\end{eqnarray}$
Here, ${\rho }_{c}={\mathrm{Tr}}_{\alpha ,\beta }\left\{\rho \right\}$ is the reduced density matrix and Da(ρ) =γa(2aρa† − a†aρ − ρa†a)/2, Db(ρ) = γb(2bρb† −b†bρ −ρb†b)/2.To characterize the quantum statistics in reservoirs, we calculate the second-order correlation function of phonons in reservoir ${G}_{\beta }^{(2)}(x)={\mathrm{lim}}_{T\to +\infty }\langle {\tilde{\beta }}_{T}^{\dagger }{\tilde{\beta }}_{T\,+\,x}^{\dagger }{\tilde{\beta }}_{T\,+\,x}{\tilde{\beta }}_{T}\rangle $ where ${\tilde{\beta }}_{x}=\int {\rm{d}}k{\beta }_{k}{{\rm{e}}}^{{\rm{i}}{kx}}/\sqrt{2\pi }$ is the Fourier transform of βk. The correlation function follows from the input–output theory and the quantum regression theorem as [45]
$\begin{eqnarray}\begin{array}{rcl}{G}_{\beta }^{(2)}(x) & \approx & {\mathrm{lim}}_{T\to +\infty }{\gamma }_{b}^{2}\mathrm{Tr}\left\{b{{\rm{e}}}^{{ \mathcal L }x}\left[b{{\rm{e}}}^{{ \mathcal L }T}\left[{\rho }_{0}\right]{b}^{\dagger }\right]{b}^{\dagger }\right\}\\ & \equiv & {\mathrm{lim}}_{T\to +\infty }{\gamma }_{b}^{2}\mathrm{Tr}\left\{b{{\rm{e}}}^{{ \mathcal L }x}[b{\rho }_{{ss}}{b}^{\dagger }]{b}^{\dagger }\right\},\end{array}\end{eqnarray}$
where ρ0 = ∣0⟩a⟨0∣ ⨂ ∣0⟩b⟨0∣ is the initial state, and ρss is the steady state satisfying ${ \mathcal L }[{\rho }_{{ss}}]=0$.3. Quantum statistical character
In the numerical calculation, we choose the cut-offs na and nb for the excitation numbers of cavity photon and phonon modes, respectively. The master equation, equation (5 ), can be numerically evaluated in the superspace [46], as shown in appendix A . The steady state ρss determines the second-order correlation function ${G}_{\beta }^{(2)}(x)$ through equation (6 ). The numerical results of the normalized correlation function ${g}_{\beta }^{(2)}(x)={G}_{\beta }^{(2)}(x)/\max \left\{{G}_{\beta }^{(2)}(x^{\prime} )\right\}$ under experimental parameters are shown in figure 2, and we also exhibit some results under non-experimental parameters in appendix B . It can be seen that ${G}_{\beta }^{(2)}(x)$ is always bunching for small ∣Δ∣. As ∣Δ∣ increases, ${G}_{\beta }^{(2)}(x)$ has the tendency to be anti-bunching. The correlation function is indeed anti-bunching for $| {\rm{\Delta }}| \gg \max \left\{{\omega }_{b},{\gamma }_{a},{\gamma }_{b},g\right\}$, however, the larger ∣Δ∣ is, the smaller the correlation function is, making the anti-bunching phenomenon difficult to observe, as we will discuss later.
Figure 2. The numerical results of the normalized correlation function ${g}_{\beta }^{(2)}(x)={G}_{\beta }^{(2)}(x)/\max \{{G}_{\beta }^{(2)}(x^{\prime} )\}$. The parameters are ωb = 1, ω = 10−7 and (a) g = 10−5, γa = 0.13, γb = 10−5; (b) g = 10−5, γa = 0.13, γb = 0.1; (c) g = 10−5, γa = 1, γb = 10−5; (d) g = 0.1, γa = 0.13, γb = 10−5. |
The approximate result in the weak driving intensity limit can be obtained analytically via the Dyson expansion9 )–(11 ) are discussed in the next section via the scattering theory.
$\begin{eqnarray}\begin{array}{rcl}{{\rm{e}}}^{{ \mathcal L }T} & \approx & {{\rm{e}}}^{{{ \mathcal L }}_{0}T}+{\displaystyle \int }_{0}^{T}{\rm{d}}{t}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{0}(T-{t}_{1})}{{ \mathcal L }}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{0}{t}_{1}}\\ & & +{\displaystyle \int }_{0}^{T}{\rm{d}}{t}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{0}(T-{t}_{1})}{\displaystyle \int }_{0}^{{t}_{1}}{\rm{d}}{t}_{2}{{ \mathcal L }}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{0}({t}_{1}-{t}_{2})}{{ \mathcal L }}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{0}{t}_{2}}\end{array}\end{eqnarray}$
to second order of ${{ \mathcal L }}_{1}$, where ${ \mathcal L }={{ \mathcal L }}_{0}+{{ \mathcal L }}_{1}$ is determined by ${{ \mathcal L }}_{0}={{ \mathcal L }}_{\mathrm{eff}}+{{ \mathcal L }}_{a}+{{ \mathcal L }}_{b}$ and the driving term ${{ \mathcal L }}_{1}\rho =-{\rm{i}}[{H}_{d},\rho ]$. ${{ \mathcal L }}_{\mathrm{eff}}\rho =-{\rm{i}}({H}_{\mathrm{eff}}\rho -\rho {H}_{\mathrm{eff}}^{\dagger })$ describes the free evolution by Heff = Hc − iγaa†a/2 − iγbb†b/2. ${{ \mathcal L }}_{a}\rho ={\gamma }_{a}a\rho {a}^{\dagger }$ and ${{ \mathcal L }}_{b}\rho \,={\gamma }_{b}b\rho {b}^{\dagger }$ characterize the emission of a photon and a phonon, respectively. The Dyson-expansion of ${{\rm{e}}}^{{{ \mathcal L }}_{0}x}$ to the first order of ${{ \mathcal L }}_{a}$ $\begin{eqnarray}{{\rm{e}}}^{{{ \mathcal L }}_{0}x}\approx {{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}x}+{\int }_{0}^{x}{\rm{d}}{t}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}(x-{t}_{1})}{{ \mathcal L }}_{a}{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{1}}\end{eqnarray}$
enable us to obtain ${G}_{\beta }^{(2)}(x)\approx {G}_{1}(x)+{G}_{2}(x)+{G}_{3}(x)$ approximately in the weak coupling limit $\left|g/({\omega }_{b}-{\rm{i}}{\gamma }_{b}/2)\right|\ll 1$. Here, $\begin{eqnarray}\begin{array}{rcl}{G}_{1}(x) & = & {\gamma }_{b}^{2}{{\rm{\Omega }}}^{2}{}_{b}\langle 0| {}_{a}\langle 0| {\displaystyle \int }_{0}^{+\infty }{\rm{d}}{t}_{1}{\rm{d}}{t}_{2}{\rm{d}}{t}_{3}b{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}x}\left[b{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{1}}\right.\\ & \times & \left.{{ \mathcal L }}_{a}{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{2}}\left[{{ \mathcal L }}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{3}}{{ \mathcal L }}_{1}\left[{\rho }_{0}\right]\right]{b}^{\dagger }\right]{b}^{\dagger }| 0{\rangle }_{a}| 0{\rangle }_{b},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{G}_{2}(x) & = & {\gamma }_{b}^{2}{{\rm{\Omega }}}^{2}{}_{b}\langle 0| {}_{a}\langle 0| {\displaystyle \int }_{0}^{x}{\rm{d}}\tau {\displaystyle \int }_{0}^{+\infty }{\rm{d}}{t}_{1}{\rm{d}}{t}_{2}b{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}(x-\tau )}{{ \mathcal L }}_{a}\\ & & \quad \times {{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}\tau }\left[b{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{1}}\left[{{ \mathcal L }}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{2}}{{ \mathcal L }}_{1}\left[{\rho }_{0}\right]\right]{b}^{\dagger }\right]{b}^{\dagger }| 0{\rangle }_{a}| 0{\rangle }_{b},\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{rcl}{G}_{3}(x) & = & {\gamma }_{b}^{2}{{\rm{\Omega }}}^{2}{}_{b}\langle 0| {}_{a}\langle 1| {\displaystyle \int }_{0}^{+\infty }{\rm{d}}{t}_{1}{\rm{d}}{t}_{2}b{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}x}\\ & & \quad \times \left[b{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{1}}\left[{{ \mathcal L }}_{1}{{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{eff}}{t}_{2}}{{ \mathcal L }}_{1}\left[{\rho }_{0}\right]\right]{b}^{\dagger }\right]{b}^{\dagger }| 1{\rangle }_{a}| 0{\rangle }_{b}\end{array}\end{eqnarray}$
describe three different processes respectively: (a) the photon is emitted before the emissions of two phonons, (b) the photon is emitted between the emissions of two phonons and (c) the photon is emitted after the emissions of two phonons. The analytical results of equations (4. Analytical results by scattering theory
Under the weak driving and coupling limit, the close connection between the master equation and the scattering theory results in the correlation function proportional to the probability distribution of two phonons in the coordinate space [47, 48]
$\begin{eqnarray}{G}_{\beta }^{(2)}(x)\propto {{\rm{\Omega }}}^{2}\int {\rm{d}}y{\left|{\rm{\Psi }}(y;{x}_{0},{x}_{0}+x)\right|}^{2},\end{eqnarray}$
where the scattering amplitude $\begin{eqnarray}{\rm{\Psi }}(y;{x}_{1},{x}_{2})=\langle 0| {\tilde{\alpha }}_{y}{\tilde{\beta }}_{{x}_{1}}{\tilde{\beta }}_{{x}_{2}}S{\alpha }_{{k}_{{\rm{i}}}}| 0\rangle \end{eqnarray}$
is determined by the Fourier transform ${\tilde{\alpha }}_{y}=\int {\rm{d}}k{{\rm{e}}}^{{\rm{i}}{ky}}{\alpha }_{k}/\sqrt{2\pi }$ of αk. The initial momentum ki is given by the frequency ωd = cki and S is the S-matrix defined as $\begin{eqnarray}\begin{array}{rcl}S & = & \mathop{\mathrm{lim}}\limits_{T\to +\infty }{{\rm{e}}}^{{\rm{i}}{H}_{r}T}{{\rm{e}}}^{-2{\rm{i}}({H}_{c}+{H}_{r}+{H}_{\mathrm{cr}})T}{{\rm{e}}}^{{\rm{i}}{H}_{r}T}\\ & = & \mathop{\mathrm{lim}}\limits_{T\to +\infty }{ \mathcal T }\left[{{\rm{e}}}^{-{\rm{i}}{\displaystyle \int }_{-T}^{+T}\tilde{H}(t){\rm{d}}t}\right],\end{array}\end{eqnarray}$
where $\tilde{H}$ is the Hamiltonian in the interaction picture $\begin{eqnarray}\begin{array}{rcl}\tilde{H}(t) & = & {H}_{{\rm{c}}}+{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{a}}{2\pi }}\displaystyle \int {\rm{d}}k({\alpha }_{k}^{\dagger }{{\rm{e}}}^{{\rm{i}}{kt}}a-{\rm{h}}.{\rm{c}}.)\\ & & +{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{b}}{2\pi }}\displaystyle \int {\rm{d}}q({\beta }_{q}^{\dagger }{{\rm{e}}}^{{\rm{i}}{qt}}b-{\rm{h}}.{\rm{c}}.).\end{array}\end{eqnarray}$
4.1. Field theoretical approach
We calculate the scattering amplitude equation (13 ) using the path integral approach, which can also be achieved by an alternative method, i.e. quantum regression theorem, as shown in appendix C [49]. The wavefunction (13 ) is determined by its Fourier transform, i.e.17 ) gives rise to the scattering amplitude ${ \mathcal A }={{ \mathcal F }}_{\mathrm{out}}^{* }{{ \mathcal F }}_{\mathrm{in}}{{ \mathcal A }}_{J}$, where ${{ \mathcal A }}_{J}$ is defined as18 ), we have written ${{ \mathcal A }}_{J}$ using the path integral formalism, where the actions of OMS and reservoirs are
$\begin{eqnarray}{ \mathcal A }=\langle 0| {\alpha }_{{k}_{f}}{\beta }_{{q}_{1}}{\beta }_{{q}_{2}}S{\alpha }_{{k}_{i}}^{\dagger }| 0\rangle .\end{eqnarray}$
The initial and final states $\begin{eqnarray}{\alpha }_{{k}_{i}}^{\dagger }| 0\rangle ={{ \mathcal F }}_{\mathrm{in}}{{\rm{e}}}^{\int {\rm{d}}{{kJ}}_{k}{\alpha }_{k}^{\dagger }+\int {\rm{d}}{{qL}}_{q}{\beta }_{q}^{\dagger }}| 0\rangle ,\end{eqnarray}$
$\begin{eqnarray}\langle 0| {\alpha }_{{k}_{f}}{\beta }_{{q}_{1}}{\beta }_{{q}_{2}}={{ \mathcal F }}_{\mathrm{out}}^{* }\langle 0| {{\rm{e}}}^{\int {\rm{d}}{{kJ}}_{k}^{* }{\alpha }_{k}+\int {\rm{d}}{{qL}}_{q}^{* }{\beta }_{q}}\end{eqnarray}$
can be written in terms of functional derivatives ${{ \mathcal F }}_{\mathrm{in}}={\mathrm{lim}}_{\left\{J,L\right\}\to 0}\delta /\delta {J}_{{k}_{i}}$ and ${{ \mathcal F }}_{\mathrm{out}}={\mathrm{lim}}_{\left\{J,L\right\}\to 0}{\delta }^{3}/$ $\delta {J}_{{k}_{f}}\delta {L}_{{q}_{1}}\delta {L}_{{q}_{2}}$. The equation ( $\begin{eqnarray}\begin{array}{l}{{ \mathcal A }}_{J}=\langle 0| {{\rm{e}}}^{\displaystyle \int {\rm{d}}{{kJ}}_{k}^{* }{\alpha }_{k}+\displaystyle \int {\rm{d}}{{qL}}_{q}^{* }{\beta }_{q}}S{{\rm{e}}}^{\displaystyle \int {\rm{d}}{{kJ}}_{k}{\alpha }_{k}^{\dagger }+\displaystyle \int {\rm{d}}{{qL}}_{q}{\beta }_{q}^{\dagger }}| 0\rangle \\ =\mathop{\mathrm{lim}}\limits_{T\to +\infty }\displaystyle \int { \mathcal D }[a,{a}^{* };b,{b}^{* }]{{\rm{e}}}^{{\rm{i}}{S}_{c}}\\ \times \displaystyle \int { \mathcal D }[\alpha ,{\alpha }^{* };\beta ,{\beta }^{* }]{{\rm{e}}}^{\sum _{k}{J}_{k}^{* }{\alpha }_{k}(T)+\sum _{q}{L}_{q}^{* }{\beta }_{q}(T)+{\rm{i}}{S}_{r}}.\end{array}\end{eqnarray}$
In the second equality in equation ( $\begin{eqnarray}\begin{array}{rcl}{S}_{c} & = & \mathop{\mathrm{lim}}\limits_{T\to +\infty }{\displaystyle \int }_{-T}^{T}{\rm{d}}t\left[{\rm{i}}{a}^{* }\dot{a}+{\rm{i}}{b}^{* }\dot{b}\right.\\ & & \left.-{\omega }_{a}{a}^{* }a-{\omega }_{b}{b}^{* }b-{{ga}}^{* }a(b+{b}^{* })\right]\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{rcl}{S}_{r} & = & \mathop{\mathrm{lim}}\limits_{T\to +\infty }{\displaystyle \int }_{-T}^{T}{\rm{d}}t\left\{\displaystyle \int {\rm{d}}k\left[\Space{0ex}{3.0ex}{0ex}{\rm{i}}{\alpha }_{k}^{* }{\dot{\alpha }}_{k}\right.\right.\\ & & \left.-{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{a}}{2\pi }}({\alpha }_{k}^{* }{{\rm{e}}}^{{\rm{i}}{kt}}a-{\rm{h}}.{\rm{c}}.)\right]\\ & & \left.+\displaystyle \int {\rm{d}}q\left[{\rm{i}}{\beta }_{q}^{* }{\dot{\beta }}_{q}-{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{b}}{2\pi }}({\beta }_{q}^{* }{{\rm{e}}}^{{\rm{i}}{qt}}b-{\rm{h}}.{\rm{c}}.)\right]\right\},\end{array}\end{eqnarray}$
respectively.The saddle point solutions of bath degrees of freedom22 ), eventually leads to13 ) can be obtained straightforwardly from the Fourier transform of ${ \mathcal A }$ as [50]
$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{k,\mathrm{cl}}(t) & = & {J}_{k}+\sqrt{\displaystyle \frac{{\gamma }_{a}}{2\pi }}{\displaystyle \int }_{-T}^{T}{\rm{d}}{ua}(u){{\rm{e}}}^{{\rm{i}}{ku}},\\ {\beta }_{q,\mathrm{cl}}(t) & = & {L}_{q}+\sqrt{\displaystyle \frac{{\gamma }_{b}}{2\pi }}{\displaystyle \int }_{-T}^{T}{\rm{d}}{ub}(u){{\rm{e}}}^{{\rm{i}}{qu}}\end{array}\end{eqnarray}$
are determined by $\delta {S}_{r}/\delta {\alpha }_{k}^{* }=\delta {S}_{r}/\delta {\beta }_{q}^{* }=0$. Following the saddle point method, we expand the fields around the classical trajectories as αk = αk,cl + δαk and βq = βq,cl + δβq, where δαk and δβq are quantum fluctuations. By integrating out fluctuation fields, we obtain $\begin{eqnarray}{{ \mathcal A }}_{J}={{\rm{e}}}^{\int {\rm{d}}k| {J}_{k}{| }^{2}+\int {\rm{d}}q| {L}_{q}{| }^{2}}\int { \mathcal D }[a,{a}^{* };b,{b}^{* }]{{\rm{e}}}^{{\rm{i}}{S}_{\mathrm{eff}}}{{\rm{e}}}^{{\rm{i}}{S}_{J}},\end{eqnarray}$
where the effective action $\begin{eqnarray}\begin{array}{rcl}{S}_{\mathrm{eff}} & = & \mathop{\mathrm{lim}}\limits_{T\to +\infty }{\displaystyle \int }_{-T}^{T}{\rm{d}}t\left[{\rm{i}}{a}^{* }\dot{a}+{\rm{i}}{b}^{* }\dot{b}-\left({\omega }_{a}-{\rm{i}}\displaystyle \frac{{\gamma }_{a}}{2}\right){a}^{* }a\right.\\ & & \left.-\left({\omega }_{b}-{\rm{i}}\displaystyle \frac{{\gamma }_{b}}{2}\right){b}^{* }b-{{ga}}^{* }a(b+{b}^{* })\right],\end{array}\end{eqnarray}$
and the source terms $\begin{eqnarray}\begin{array}{rcl}{S}_{J} & = & -\mathop{\mathrm{lim}}\limits_{T\to +\infty }{\displaystyle \int }_{-T}^{T}{\rm{d}}t\left[{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{a}}{2\pi }}\displaystyle \int {\rm{d}}k({J}_{k}^{* }a{{\rm{e}}}^{{\rm{i}}{kt}}-{\rm{h}}.{\rm{c}}.)\right.\\ & & \left.+{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{b}}{2\pi }}\displaystyle \int {\rm{d}}q({L}_{q}^{* }b{{\rm{e}}}^{{\rm{i}}{qt}}-{\rm{h}}.{\rm{c}}.)\right].\end{array}\end{eqnarray}$
The functional derivatives ${{ \mathcal F }}_{\mathrm{out}}$ and ${{ \mathcal F }}_{\mathrm{in}}$ to ${{ \mathcal A }}_{J}$, i.e. equation ( $\begin{eqnarray}\begin{array}{rcl}{ \mathcal A } & = & -\displaystyle \frac{{\gamma }_{a}{\gamma }_{b}}{4{\pi }^{2}}\displaystyle \int {\rm{d}}{t}_{{k}_{i}}{\rm{d}}{t}_{{k}_{f}}{\rm{d}}{t}_{{q}_{1}}{\rm{d}}{t}_{{q}_{2}}{{\rm{e}}}^{{\rm{i}}{k}_{f}{t}_{{k}_{f}}+{\rm{i}}{q}_{1}{t}_{{q}_{1}}+{\rm{i}}{q}_{2}{t}_{{q}_{2}}-{\rm{i}}{k}_{i}{t}_{{k}_{i}}}\\ & & \times \langle 0| { \mathcal T }{\left[a({t}_{{k}_{f}})b({t}_{{q}_{1}})b({t}_{{q}_{2}}){a}^{\dagger }({t}_{{k}_{i}})\right]}_{\mathrm{eff}}| 0\rangle ,\end{array}\end{eqnarray}$
where the footnote ‘eff’ denotes the evolution of operators ${\left[A(t)\right]}_{\mathrm{eff}}={{\rm{e}}}^{{\rm{i}}{H}_{\mathrm{eff}}t}A{{\rm{e}}}^{-{\rm{i}}{H}_{\mathrm{eff}}t}$ governed by the effective Hamiltonian Heff = Hc − iγaa†a/2 − iγbb†b/2. The scattering amplitude in equation ( $\begin{eqnarray}{\rm{\Psi }}(y;{x}_{1},{x}_{2})=\displaystyle \frac{1}{{\left(2\pi \right)}^{\tfrac{3}{2}}}\int {\rm{d}}{k}_{f}{\rm{d}}{q}_{1}{\rm{d}}{q}_{2}{{\rm{e}}}^{{\rm{i}}{{yk}}_{f}+{\rm{i}}{x}_{1}{q}_{1}+{\rm{i}}{x}_{2}{q}_{2}}{ \mathcal A }.\end{eqnarray}$
Due to the identicality of two emitting phonons, there are only 3 time orders in equation (25 ), which correspond to the scattering processes (a), (b) and (c) in section 3 and give rise to the same results in equations (9 )–(11 ) up to a global factor γa/2πω2. The detailed calculations are left in appendix C . In figure 3, we present the relative differences $({\left[{g}_{\beta }^{(2)}(x)\right]}_{\mathrm{numerical}}-{\left[{g}_{\beta }^{(2)}(x)\right]}_{\mathrm{analytical}})/{\left[{g}_{\beta }^{(2)}(x)\right]}_{\mathrm{numerical}}$ between the results in figure 2 and those based on the scattering theory, which show a good agreement between the numerical and analytical results. To show the relationship between the photon emission and the phonon emission clearly, the representative transitions among the energy levels of the system corresponding to the three scattering processes are depicted in figure 4.
Figure 4. Representative transitions among the energy levels of the system corresponding to the three scattering processes, where the subfigures (a)–(c) correspond to processes (a)–(c) respectively. We have omitted the higher-order transitions of g due to the weak coupling condition. |
4.2. Results
The intuitive picture of the quantum statistical characters is shown in the following. G1(x) displays the bunching behavior, since it represents the two-phonon cascade decays after the emission of one photon; G2(x) displays the anti-bunching behavior since it represents the emissions of the two phonons are separated by the emission of the single photon. G3(x) displays the bunching behavior for small ∣Δ∣ and the anti-bunching behavior for large ∣Δ∣. The intriguing behavior of G3(x) can be interpreted as follows. In the process represented by G3(x), two phonons are emitted before the emission of the photon. As result, the presence of the photon in the cavity induces the destructive interference of two phonons at x = 0 for large ∣Δ∣.
We note that the larger ∣Δ∣ is, the larger proportion of processes (b) and (c) occupy among (a)–(c). To illustrate it, we define ${g}_{i}={\int }_{0}^{+\infty }{\rm{d}}{{xG}}_{i}(x)$, and show the ratios ri = gi/∑jgj in figure 5, where r1 reaches its peak while r2 and r3 reach their minima at Δ ≈ ωb.
Figure 5. The relative proportion of ${\int }_{0}^{+\infty }{\rm{d}}{{xG}}_{1\sim 3}(x)$ to their sum. The proportion of process (a)–(c) are given by ${\int }_{0}^{+\infty }{\rm{d}}{{xG}}_{1\sim 3}(x)/{\sum }_{n=1}^{3}{\int }_{0}^{+\infty }{\rm{d}}{{xG}}_{n}(x)$, respectively. The parameters are ωb = 1, g = 10−5, γa = 0.13, γb = 10−5. |
Because of the large proportion of the process (b) in the large ∣Δ∣ limit, the two emitting phonons show anti-bunching behavior, and the large ∣Δ∣ also indicates an inresonant scattering of the incident photon with a small emission rate (lower than 10−6) of two phonons.
For small ∣Δ∣, the large proportion of the process (a) indicates the bunching behavior of two emitting phonons. Besides, due to the small detuning, the driving field resonates with the system, making the emission rate of two phonons relatively larger but still not high (lower than 10−3).
If the weak coupling condition is released, e.g. ∣g/(ω − iγb/2)∣ ∼ 1, it turns out that the emissivity of two bunched phonons can be enhanced by one to two orders of magnitudes, i.e. 10−2–10−1, while the emission rate of two anti-bunched phonons can reach 10−4.
5. Conclusion and outlook
In summary, we have identified the quantum statistical behavior of two phonons in the weak coupled OMS under a weak driving field. We calculate the second-order correlation functions by using the standard master equation approach and the scattering formalism. We obtain the analytic results for the correlation functions by identifying three dominant physical processes. Our results show that the small detuning induces the bunching behavior of phonons with a relatively large emission rate. The anti-bunching behavior of phonons requires large detuning, which leads to a small two-phonon emissivity. To make the manifest anti-bunching phenomenon two feasible schemes, increasing the coupling strength and the driving intensity are proposed, which require more studies on multiphonon and multiphoton scatterings in future work.
Acknowledgments
We are grateful to Junqiao Pan, Yuqi Wang and Yaoqi Tian for useful discussions. T.Shi acknowledges the support from NSFC Grant No. 11 974 363.
Appendix A. Mapping of operators to superspace
Consider a formula in the form of AρB where A, B are operators and ρ is a density matrix. Make the following mappings5 ) can be mapped toA3 ) enables us to evaluate the master equation numerically.
$\begin{eqnarray}\rho =\sum _{m,n}{\rho }_{{mn}}| m\rangle \langle n| \to \sum _{m,n}{\rho }_{{mn}}| m\rangle | n\rangle \equiv | \rho \rangle ,\end{eqnarray}$
and $\begin{eqnarray}A(\cdot )B\to {B}^{{\rm{T}}}\otimes A(\cdot ).\end{eqnarray}$
So the master equation ( $\begin{eqnarray}\begin{array}{l}{\partial }_{t}| \rho \rangle \\ =[-{\rm{i}}(I\otimes {H}_{e}-{H}_{e}^{\dagger }\otimes I)+{\gamma }_{a}{a}^{* }\otimes a+{\gamma }_{b}{b}^{* }\otimes b]| \rho \rangle ,\end{array}\end{eqnarray}$
where He = Hc + Hd − iγaa†a/2 − iγbb†b/2. Equation (Appendix B. Numerical results with non-experimental parameters
We have shown the numerical results of the normalized correlation functions under experimental parameters in figure 2, and here we give some other results under non-experimental parameters in figure B1 in order to show the validity of our conclusion in a wider range of parameters. In figure B1 we can see that the phonons also show an anti-bunching behavior as ∣Δ∣ increases to greater than other parameters.
Figure B1. Some other numerical results of the normalized correlation function ${g}_{\beta }^{(2)}(x)={G}_{\beta }^{(2)}(x)/\max \{{G}_{\beta }^{(2)}(x^{\prime} )\}$ with non-experimental parameters. The parameters are ωb = 1, ω = 10−7 and (a) g = 0.01, γa = 0.02, γb = 0.002; (b) g = 0.3, γa = 20, γb = 5; (c) g = 0.05, γa = 1, γb = 0.05; (d) g = 0.3, γa = 2, γb = 5. |
Appendix C. Scattering amplitude by quantum regression theorem
We introduce another method by using the quantum regression theorem to calculate the scattering amplitude equation (13 ). Let’s start with equation (17 ), which can be written in another form14 ) and (C1 ) into (16 ) and coming back from the rotating frame with respect to Hr, we can write the scattering amplitude asC4 ) give a Green function, and we get
$\begin{eqnarray}{\alpha }_{{k}_{i}}^{\dagger }| 0\rangle ={{ \mathcal F }}_{\mathrm{in}}{{\rm{e}}}^{\tfrac{1}{2}(\int {\rm{d}}k| {J}_{k}{| }^{2}+{\rm{d}}q| {L}_{q}{| }^{2})}D(\{J,L\})| 0\rangle ,\end{eqnarray}$
$\begin{eqnarray}\langle 0| {\alpha }_{{k}_{f}}{\beta }_{{q}_{1}}{\beta }_{{q}_{2}}={{ \mathcal F }}_{\mathrm{out}}^{* }\langle 0| {D}^{\dagger }(\{J,L\}){{\rm{e}}}^{\tfrac{1}{2}(\int {\rm{d}}k| {J}_{k}{| }^{2}+{\rm{d}}q| {L}_{q}{| }^{2})},\end{eqnarray}$
where the translation operator is $\begin{eqnarray}D(\{J,L\})={{\rm{e}}}^{\int {\rm{d}}k({J}_{k}{\alpha }_{k}^{\dagger }-{J}_{k}^{* }{\alpha }_{k})}{{\rm{e}}}^{\int {\rm{d}}q({L}_{q}{\beta }_{q}^{\dagger }-{L}_{q}^{* }{\beta }_{q})},\end{eqnarray}$
which can transform the Hamiltonian into ${D}^{\dagger }(\{J,L\})\tilde{H}D(\{J,L\})=\tilde{H}+{H}_{D}(t)$ and $\begin{eqnarray}\begin{array}{rcl}{H}_{D}(t) & = & {\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{a}}{2\pi }}\displaystyle \int {\rm{d}}k\left({J}_{k}^{* }{{\rm{e}}}^{{\rm{i}}{kt}}a-{J}_{k}{{\rm{e}}}^{-{\rm{i}}{kt}}{a}^{\dagger }\right)\\ & & +\,{\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{b}}{2\pi }}\displaystyle \int {\rm{d}}q\left({L}_{q}^{* }{{\rm{e}}}^{{\rm{i}}{qt}}b-{L}_{q}{{\rm{e}}}^{-{\rm{i}}{qt}}{b}^{\dagger }\right)\end{array}\end{eqnarray}$
is the driving term. Substituting equations ( $\begin{eqnarray}\begin{array}{l}{ \mathcal A }=\mathop{\mathrm{lim}}\limits_{\displaystyle \genfrac{}{}{0em}{}{T\to +\infty }{\left\{J,{J}^{* },L,{L}^{* }\right\}\to 0}}\displaystyle \frac{{\delta }^{3}}{\delta {J}_{{k}_{f}}^{* }\delta {L}_{{q}_{1}}^{* }\delta {L}_{{q}_{2}}^{* }}\displaystyle \frac{\delta }{\delta {J}_{{k}_{i}}}\\ \times \langle 0| { \mathcal T }\left[{{\rm{e}}}^{-{\rm{i}}{\displaystyle \int }_{-T}^{T}[H+{H}_{D}(t)]{\rm{d}}t}\right]| 0\rangle {{\rm{e}}}^{\displaystyle \int {\rm{d}}k| {J}_{k}{| }^{2}+\displaystyle \int {\rm{d}}q| {L}_{q}{| }^{2}}.\end{array}\end{eqnarray}$
The functional derivatives in equation ( $\begin{eqnarray}\begin{array}{rcl}{ \mathcal A } & = & -\displaystyle \frac{{\gamma }_{a}{\gamma }_{b}}{4{\pi }^{2}}\displaystyle \int {\rm{d}}{t}_{{k}_{i}}{\rm{d}}{t}_{{k}_{f}}{\rm{d}}{t}_{{q}_{1}}{\rm{d}}{t}_{{q}_{2}}{{\rm{e}}}^{{\rm{i}}{k}_{f}{t}_{{k}_{f}}+{\rm{i}}{q}_{1}{t}_{{q}_{1}}+{\rm{i}}{q}_{2}{t}_{{q}_{2}}-{\rm{i}}{k}_{i}{t}_{{k}_{i}}}\\ & & \times \langle 0| { \mathcal T }[a({t}_{{k}_{f}})b({t}_{{q}_{1}})b({t}_{{q}_{2}}){a}^{\dagger }({t}_{{k}_{i}})]| 0\rangle ,\end{array}\end{eqnarray}$
where the operators in Heisenberg picture A(t) = eiHtAe−iHt are governed by H. Then let’s study the Green function $\langle 0| { \mathcal T }[a({t}_{{k}_{f}})b({t}_{{q}_{1}})b({t}_{{q}_{2}}){a}^{\dagger }({t}_{{k}_{i}})]| 0\rangle $. Now that $b({t}_{{q}_{1}})$ and $b({t}_{{q}_{2}})$ are commutative symmetric and a, b vanish the vaccum state, there are 3 nonvanishing orders for the 4 operators in the time-ordering bracket. We take ${t}_{{k}_{f}}\gt {t}_{{q}_{1}}\gt {t}_{{q}_{2}}\gt {t}_{{k}_{i}}$ for example. the Green function turns to be $\begin{eqnarray}\begin{array}{l}\langle 0| a({t}_{{k}_{f}})b({t}_{{q}_{1}})b({t}_{{q}_{2}}){a}^{\dagger }({t}_{{k}_{i}})| 0\rangle \\ =\,\langle 0| a{{\rm{e}}}^{-{\rm{i}}H({t}_{{k}_{f}}-{t}_{{q}_{1}})}b{{\rm{e}}}^{-{\rm{i}}H({t}_{{q}_{1}}-{t}_{{q}_{2}})}b{{\rm{e}}}^{-{\rm{i}}H({t}_{{q}_{2}}-{t}_{{k}_{i}})}{a}^{\dagger }| 0\rangle \\ =\,\mathrm{Tr}\left\{a{{\rm{e}}}^{-{\rm{i}}H({t}_{{k}_{f}}-{t}_{{q}_{1}})}b{{\rm{e}}}^{-{\rm{i}}H({t}_{{q}_{1}}-{t}_{{q}_{2}})}b{{\rm{e}}}^{-{\rm{i}}H({t}_{{q}_{2}}-{t}_{{k}_{i}})}{a}^{\dagger }\right.\\ \left.\times \,| 0\rangle \langle 0| {{\rm{e}}}^{{\rm{i}}H({t}_{{q}_{2}}-{t}_{{k}_{i}})}{{\rm{e}}}^{{\rm{i}}H({t}_{{q}_{1}}-{t}_{{q}_{2}})}{{\rm{e}}}^{{\rm{i}}H({t}_{{k}_{{\rm{f}}}}-{t}_{{q}_{1}})}\right\}\\ =\,{\mathrm{Tr}}_{c}\left\{a{{\rm{e}}}^{{ \mathcal L }^{\prime} ({t}_{{k}_{f}}-{t}_{{q}_{1}})}b{{\rm{e}}}^{{ \mathcal L }^{\prime} ({t}_{{q}_{1}}-{t}_{{q}_{2}})}b{{\rm{e}}}^{{ \mathcal L }^{\prime} ({t}_{{q}_{2}}-{t}_{{k}_{i}})}{a}^{\dagger }| 0\rangle \langle 0| \right\},\end{array}\end{eqnarray}$
where in the last equality we trace out the reservoir modes and apply the master equation, which gives $\begin{eqnarray}{ \mathcal L }^{\prime} (\cdot )=-{\rm{i}}[{H}_{\mathrm{eff}}(\cdot )-(\cdot ){H}_{\mathrm{eff}}^{\dagger }]+{\gamma }_{a}a(\cdot ){a}^{\dagger }+{\gamma }_{b}(\cdot ){b}^{\dagger },\end{eqnarray}$
and we restate that the non-Hermitian effective Hamiltonian is defined as ${H}_{\mathrm{eff}}={H}_{c}-{\rm{i}}\displaystyle \frac{{\gamma }_{a}}{2}{a}^{\dagger }a-{\rm{i}}\tfrac{{\gamma }_{b}}{2}{b}^{\dagger }b$.The vacuum state ⟨0∣ is annihilated by the operators a† and b†, which prevents any contribution from the jumps induced by the last two terms in equation (C7 ). So we can rewrite equation (C6 ) asC5 ) can be rewritten as25 ). The scattering amplitude in equation (13 ) is directly obtained from the Fourier transformC10 ) divides the possibility distribution of phonons into 3 parts as $\int {\rm{d}}y| {\rm{\Psi }}(y;{x}_{1},{x}_{2}){| }^{2}={\sum }_{n=1}^{3}{\tilde{G}}_{n}({x}_{1},{x}_{2})$, where
$\begin{eqnarray}\begin{array}{l}\langle 0| a({t}_{{k}_{f}})b({t}_{{q}_{1}})b({t}_{{q}_{2}}){a}^{\dagger }({t}_{{k}_{i}})| 0\rangle \\ =\langle 0| a{{\rm{e}}}^{-{\rm{i}}{H}_{\mathrm{eff}}({t}_{{k}_{f}}-{t}_{{q}_{1}})}b{{\rm{e}}}^{-{\rm{i}}{H}_{\mathrm{eff}}({t}_{{q}_{1}}-{t}_{{q}_{2}})}b{{\rm{e}}}^{-{\rm{i}}{H}_{\mathrm{eff}}({t}_{{q}_{2}}-{t}_{{k}_{i}})}{a}^{\dagger }| 0\rangle \\ =\langle 0| {\left[a({t}_{{k}_{f}})b({t}_{{q}_{1}})b({t}_{{q}_{2}}){a}^{\dagger }({t}_{{k}_{i}})\right]}_{\mathrm{eff}}| 0\rangle ,\end{array}\end{eqnarray}$
where the footnote ‘eff’ denotes the operators are governed by the effective Hamiltonian. Following the same steps, equation ( $\begin{eqnarray}\begin{array}{rcl}{ \mathcal A } & = & -\displaystyle \frac{{\gamma }_{a}{\gamma }_{b}}{4{\pi }^{2}}\displaystyle \int {\rm{d}}{t}_{{k}_{i}}{\rm{d}}{t}_{{k}_{f}}{\rm{d}}{t}_{{q}_{1}}{\rm{d}}{t}_{{q}_{2}}{{\rm{e}}}^{{\rm{i}}{k}_{f}{t}_{{k}_{f}}+{\rm{i}}{q}_{1}{t}_{{q}_{1}}+{\rm{i}}{q}_{2}{t}_{{q}_{2}}-{\rm{i}}{k}_{i}{t}_{{k}_{i}}}\\ & & \times \langle 0| { \mathcal T }\,{\left[a({t}_{{k}_{f}})b({t}_{{q}_{1}})b({t}_{{q}_{2}}){a}^{\dagger }({t}_{{k}_{i}})\right]}_{\mathrm{eff}}| 0\rangle ,\end{array}\end{eqnarray}$
which is exactly equation ( $\begin{eqnarray}\begin{array}{rcl}{\rm{\Psi }}(y;{x}_{1},{x}_{2}) & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{\tfrac{3}{2}}}\displaystyle \int {\rm{d}}{k}_{f}{\rm{d}}{q}_{1}{\rm{d}}{q}_{2}{{\rm{e}}}^{{\rm{i}}{{yk}}_{f}+{\rm{i}}{x}_{1}{q}_{1}+{\rm{i}}{x}_{2}{q}_{2}}{ \mathcal A }\\ & = & -\displaystyle \frac{{\gamma }_{a}{\gamma }_{b}}{\sqrt{2\pi }}\displaystyle \int {\rm{d}}{t}_{{k}_{i}}{{\rm{e}}}^{-{\rm{i}}{k}_{i}{t}_{{k}_{i}}}\\ & & \times \langle 0| { \mathcal T }{\left[a(-y)b(-{x}_{1})b(-{x}_{2}){a}^{\dagger }({t}_{{k}_{i}})\right]}_{\mathrm{eff}}| 0\rangle .\end{array}\end{eqnarray}$
The possibility distribution of the two phonons is given by ∫dy∣ψ(y; x1, x2)∣2. Omitting the exchange of b(−x1) and b(−x2), the 3 nonvanishing time order in equation ( $\begin{eqnarray}{\tilde{G}}_{1}(x)={\int }_{{x}_{0}+| x| }^{+\infty }{\rm{d}}y| {\rm{\Psi }}(y;{x}_{0},{x}_{0}+| x| ){| }^{2},\end{eqnarray}$
$\begin{eqnarray}{\tilde{G}}_{2}(x)={\int }_{{x}_{0}}^{{x}_{0}+| x| }{\rm{d}}y| {\rm{\Psi }}(y;{x}_{0},{x}_{0}+| x| ){| }^{2},\end{eqnarray}$
$\begin{eqnarray}{\tilde{G}}_{3}(x)={\int }_{-\infty }^{{x}_{0}}{\rm{d}}y| {\rm{\Psi }}(y;{x}_{0},{x}_{0}+| x| ){| }^{2}.\end{eqnarray}$
In the above equations, we denote x1 and x2 as x0 and x0 + ∣x∣ for convenience, and omit the parameter x0 in ${\tilde{G}}_{n}(x)$ due to the translation invariance of scattering amplitude. The three terms respectively represent 3 scattering processes: ${\tilde{G}}_{1}(x)$ implies the photon is emitted before the two phonons, ${\tilde{G}}_{2}(x)$ implies the photon is emitted between the two phonons and ${\tilde{G}}_{3}(x)$ implies the photon is emitted after the two phonons. The three processes are corresponding to that in equations (9 )–(11 ). In fact, if we substitute equation (C10 ) into equation (C11 ) and change the integral variables and integral sequences appropriately, we can easily check that ${\tilde{G}}_{n}(x)={\gamma }_{a}{G}_{n}(x)/2\pi {{\rm{\Omega }}}^{2}$.
Finally, we directly give the expression of equation (C11 )
$\begin{eqnarray}{\tilde{G}}_{1}(x)=\displaystyle \frac{{\gamma }_{a}^{2}{\gamma }_{b}}{4\pi }{{\rm{e}}}^{-{\gamma }_{b}| x| }{\left|{}_{b}\langle 0| {b}^{2}\displaystyle \frac{-{\rm{i}}}{{H}_{1a}-{k}_{i}}| 0{\rangle }_{b}\right|}^{2},\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{\tilde{G}}_{2}(x)=\displaystyle \frac{{\gamma }_{a}^{2}{\gamma }_{b}^{2}}{2\pi }\displaystyle \sum _{{j}_{1},{j}_{2}}\displaystyle \frac{-{\rm{i}}\left[{{\rm{e}}}^{{\rm{i}}({E}_{{j}_{1}}^{* }-{E}_{{j}_{2}})| x| }-{{\rm{e}}}^{-{\gamma }_{b}| x| }\right]}{{E}_{{j}_{1}}^{* }-{E}_{{j}_{2}}-{\rm{i}}{\gamma }_{b}}\\ \times {\left({}_{b}\langle 0| {bV}| {j}_{1}\rangle \langle {j}_{1}| {V}^{-1}b\displaystyle \frac{-{\rm{i}}}{{H}_{1a}-{k}_{i}}| 0{\rangle }_{b}\right)}^{* }\\ \times \left({}_{b}\langle 0| {bV}| {j}_{2}\rangle \langle {j}_{2}| {V}^{-1}b\displaystyle \frac{-{\rm{i}}}{{H}_{1a}-{k}_{i}}| 0{\rangle }_{b}\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{\tilde{G}}_{3}(x)=\displaystyle \frac{{\gamma }_{a}^{2}{\gamma }_{b}^{2}}{2\pi }\displaystyle \sum _{{j}_{1},{j}_{2}}\displaystyle \sum _{{d}_{1},{d}_{2}=0}^{1}\displaystyle \frac{{\mathrm{ie}}^{{\rm{i}}\left({E}_{{j}_{1}+{d}_{1}}^{* }-{E}_{{j}_{2}+{d}_{2}}\right)| x| }}{{E}_{{j}_{1}}^{* }-{E}_{{j}_{2}}}\\ \times \left({}_{b}\langle 0| V| {j}_{1}{\rangle }_{b}\langle {j}_{1}| \left(b-\displaystyle \frac{2g}{2{\omega }_{b}-{\rm{i}}{\gamma }_{b}}\right)| {j}_{1}+{d}_{1}{\rangle }_{b}\langle {j}_{1}+{d}_{1}| \right.\\ {\left.\times {V}^{-1}b\displaystyle \frac{-{\rm{i}}}{{H}_{1a}-{k}_{i}}| 0{\rangle }_{b}\right)}^{* }\left({}_{b}\langle 0| V| {j}_{2}{\rangle }_{b}\langle {j}_{2}| \left(b-\displaystyle \frac{2g}{2{\omega }_{b}-{\rm{i}}{\gamma }_{b}}\right)\right.\\ \left.\times | {j}_{2}+{d}_{2}{\rangle }_{b}{\left\langle {j}_{2}+{d}_{2}| {V}^{-1}b\displaystyle \frac{-{\rm{i}}}{{H}_{1a}-{k}_{i}}| 0\right\rangle }_{b}\right),\end{array}\end{eqnarray}$
where $\begin{eqnarray}\begin{array}{rcl}{H}_{1a} & = & {}_{a}\langle 1| {H}_{\mathrm{eff}}| 1{\rangle }_{a}\\ & = & \left({\omega }_{b}-{\rm{i}}\displaystyle \frac{{\gamma }_{b}}{2}\right){b}^{\dagger }b+g(b+{b}^{\dagger })+{\omega }_{a}-{\rm{i}}\displaystyle \frac{{\gamma }_{a}}{2},\end{array}\end{eqnarray}$
and $\begin{eqnarray}V={{\rm{e}}}^{-\tfrac{2g}{2{\omega }_{b}-{\rm{i}}{\gamma }_{b}}{b}^{\dagger }+\displaystyle \frac{2g}{2{\omega }_{b}-{\rm{i}}{\gamma }_{b}}b}\end{eqnarray}$
can diagonalize H1a as $\begin{eqnarray}\begin{array}{rcl}{V}^{-1}{H}_{1a}V & = & \left({\omega }_{b}-{\rm{i}}\displaystyle \frac{{\gamma }_{b}}{2}\right){b}^{\dagger }b+{\omega }_{a}-\displaystyle \frac{4{g}^{2}{\omega }_{b}}{4{\omega }_{b}^{2}+{\gamma }_{b}^{2}}\\ & & -{\rm{i}}\left(\displaystyle \frac{{\gamma }_{a}}{2}+\displaystyle \frac{2{g}^{2}{\gamma }_{b}}{4{\omega }_{b}^{2}+{\gamma }_{b}^{2}}\right)\end{array}\end{eqnarray}$
whose eigenvalue is $\begin{eqnarray}\begin{array}{l}{E}_{j}=j\left({\omega }_{b}-{\rm{i}}\displaystyle \frac{{\gamma }_{b}}{2}\right)+{\omega }_{a}\\ -\displaystyle \frac{4{g}^{2}{\omega }_{b}}{4{\omega }_{b}^{2}+{\gamma }_{b}^{2}}-{\rm{i}}\left(\displaystyle \frac{{\gamma }_{a}}{2}+\displaystyle \frac{2{g}^{2}{\gamma }_{b}}{4{\omega }_{b}^{2}+{\gamma }_{b}^{2}}\right).\end{array}\end{eqnarray}$