Skip to content
BY 4.0 license Open Access Published by De Gruyter September 15, 2020

Generation and dynamics of entangled fermion–photon–phonon states in nanocavities

  • Mikhail Tokman , Maria Erukhimova , Yongrui Wang , Qianfan Chen and Alexey Belyanin ORCID logo EMAIL logo
From the journal Nanophotonics

Abstract

We develop the analytic theory describing the formation and evolution of entangled quantum states for a fermionic quantum emitter coupled simultaneously to a quantized electromagnetic field in a nanocavity and quantized phonon or mechanical vibrational modes. The theory is applicable to a broad range of cavity quantum optomechanics problems and emerging research on plasmonic nanocavities coupled to single molecules and other quantum emitters. The optimal conditions for a tripartite entanglement are realized near the parametric resonances in a coupled system. The model includes dissipation and decoherence effects due to coupling of the fermion, photon, and phonon subsystems to their dissipative reservoirs within the stochastic evolution approach, which is derived from the Heisenberg–Langevin formalism. Our theory provides analytic expressions for the time evolution of the quantum state and observables and the emission spectra. The limit of a classical acoustic pumping and the interplay between parametric and standard one-photon resonances are analyzed.

1 Introduction

There is a lot of recent interest in the quantum dynamics of fermion systems coupled to both an electromagnetic (EM) mode in a cavity and quantum or classical mechanical/acoustic oscillations or phonon vibrations. This problem is related to the burgeoning fields of cavity optomechanics [1], [2], [3] and quantum acoustics [4], [5], [6]. Another example where this situation can be realized is a molecule placed in a plasmonic nanocavity [7], [8]. In this case, the fermion system may comprise two or more electron states forming an optical transition, whereas the phonon field is simply a vibrational mode of a molecule. One can also imagine a situation where a quantum emitter such as a quantum dot (QD) or an optically active defect in a solid matrix is coupled to the quantized phonon modes of a crystal lattice, which would be an extension of an extremely active field of research on phonon–polaritons or plasmon–phonon–polaritons [9], [10] into a fully quantum regime.

Apart from the fundamental interest, the studies of such systems are motivated by quantum information applications. Indeed, the presence of a classical or quantized acoustic mode provides an extra handle to control the quantum state of a coupled fermion–boson quantum system. In the extreme quantum limit in which the fermionic degree of freedom and all bosonic degrees of freedom (both photons and phonons) are quantized, a strong enough coupling between them leads to an entangled fermion–photon–phonon state, which is a complex enough system to implement basic gates for quantum computation or other applications. Such a system has not been realized experimentally. However, many ingredients have been already demonstrated, such as strong coupling between a nanocavity mode and a single molecule [11], numerous examples of strong coupling between nanocavity modes and a single fermionic quantum emitter such as a color center [12] or a QD (see e.g., the studies by Yoshie et al. [13] and Reithmaier et al. [14] for semiconductor cavity–QD systems and the studies by Leng et al. [15], Bitton et al. [16], and Park et al. [17] for plasmonic cavities), strong coupling and entanglement of acoustic phonons [18], [19], resolving the energy levels of a nanomechanical oscillator [6], or cooling a macroscopic system into its motional ground state [20].

Interaction of three or more modes of oscillations, whether they are classical or quantized, is strongly enhanced close to the parametric resonance, which is therefore the most interesting region to study. Fortunately for theorists, the analysis near the parametric resonance is greatly simplified because some form of a slowly varying amplitude method for classical systems [21], [22] or the rotating wave approximation (RWA) for quantum systems [23] can be applied. The use of RWA restricts the coupling strength to the values much lower than the characteristic energies in the system, such as the optical transition or vibrational energy. The emerging studies of the so-called ultra-strong coupling regime [24] have to go beyond the RWA. Nevertheless, for the vast majority of experiments, including nonperturbative strong coupling dynamics and entanglement, the RWA is adequate and provides some crucial simplifications that allow one to obtain analytic solutions.

In particular, within Schrödinger’s description, the equations of motion for the components of an infinitely dimensional state vector |Ψ that describes a coupled fermion–boson system can be split into the blocks of low dimensions if the RWA is applied. This is true even if the dynamics of the fermion subsystem is nonperturbative, e.g., the effects of saturation are important. Note that there is no such simplification in the Heisenberg representation, i.e., when solving the equations of motion

(1)ddtgˆ=i[Hˆ,gˆ],

where gˆ is the Heisenberg operator of a certain physical observable g and Hˆ is the Hamiltonian of the system. Operator-valued Eq. (1) is generally impossible to split into smaller blocks, even within the RWA. This happens because some matrix elements gAB(t) of the Heisenberg operator are determined by states |A, |B which belong to different blocks that evolve independently in the Schrödinger picture. The simplification could only be possible for specially selected initial conditions in which the Heisenberg operator is determined on a “truncated” basis belonging to only one of the independent blocks. The Schrödinger’s approach also leads to fewer equations for the state vector components than the approach based on the von Neumann master equation for the elements of the density matrix.

Obviously, the Schrödinger equation in its standard form cannot be applied to describe open systems coupled to a dissipative reservoir. In this case, the stochastic versions of the equation of evolution for the state vector have been developed, e.g., the method of quantum jumps [23], [25]. This method is optimal for numerical analysis in the Monte-Carlo type schemes. Here, we formulate the stochastic equation for the state vector derived from the Heisenberg–Langevin approach which is more conducive to the analytic treatment. Its key element is an assumption that there exists the operator of evolution Uˆ, which is determined unambiguously not only by the parameters of the dynamical system but also by the statistical properties of a dissipative reservoir.

There were a number of theoretical studies of tripartite entanglement in open optomechanical systems, either for purely bosonic field modes or involving the atomic degree of freedom (see, e.g., [26], [27], [28], [29], [30], [31], [32]). The existing studies are based on either the Heisenberg–Langevin approach or the master equation. In these cases, the analytic solution is possible only after some drastic approximations such as adiabatic elimination of some degrees of freedom or within the linear perturbation theory, when the atomic populations are unperturbed. The work closest to our model is the study by Liao et al. [31], which considers tripartite entanglement in the vicinity of a parametric resonance by numerically solving the Lindblad master equation. The present paper is different in several important aspects. First, our work develops a new formalism based on the stochastic evolution of the state vector, which allowed us to obtain explicit analytic expressions for the evolution of the quantum states and all relevant observables: field and atom energies, emission spectra, etc. We were also able to derive analytic criteria for the separation of resonances, which is important for real experimental situations in which the frequency of the phonons or mechanical oscillations is much lower than the photon or electron transition frequencies, so the parametric and one-photon resonances can easily overlap.

Second, the study by Liao et al. [31] takes into account only the dissipation in the atomic subsystem, whereas our work includes relaxation and fluctuations in all subsystems: electronic, EM, and phonon/mechanical. We were able to write explicit expressions for the relaxation and noise terms which include all relevant relaxation channels and ensure that the system goes into a physically meaningful equilibrium in the absence of external excitation. Moreover, our analysis takes into account both inelastic processes (dissipation) and purely elastic decoherence processes.

The paper is structured as follows. Section II formulates the model and the Hamiltonian for coupled quantized fermion, photon, and phonon fields in a nanocavity. Section III derives the solution for the quantum states of a closed system in the vicinity of a parametric resonance and analyzes its properties. In Section IV, we provide the stochastic equation describing the evolution of quantum states of an open system in contact with a dissipative reservoir and describe the observables. In Section V, we consider the case of a classical acoustic pumping. Section VI describes the interplay of parametric and standard one-photon resonances and provides the conditions under which these resonances can be separated. Section VII gives an example of manipulating entangled electron–photon states by an acoustic pumping. Appendix contains the derivation of the stochastic equation of evolution from the Heisenberg–Langevin approach and compares with Lindblad density matrix formalism.

The focus of the paper is to provide analytic solutions for the quantum dynamics in systems of coupled electron, photon, and phonon excitations including dissipation and decoherence effects. Here, we emphasize “analytic” which means that we provide the expressions for the time evolution of the state vector and observable quantities in the form which shows explicitly the dependence on all experimental parameters: transition energies and frequencies, matrix elements of the optical transitions, the spatial structure of the field modes, relaxation rates for all constituent subsystems, ambient temperatures, etc. That is why we believe that the results obtained in this paper will be useful for the experimentalists working on a broad range of nanophotonic systems.

2 A coupled quantized electron–photon–phonon system: the model

Consider a quantized electron system coupled to the quantum EM field of a nanocavity and classical or quantized vibrational (phonon) modes, see Figure 1 which sketches two out of many possible scenarios.

Figure 1: (a) A sketch of a molecule in a nanocavity created by a metallic nanotip and a substrate; (b) A sketch of a quantum dot coupled to optical and mechanical vibrational modes in a nanocavity.
Figure 1:

(a) A sketch of a molecule in a nanocavity created by a metallic nanotip and a substrate; (b) A sketch of a quantum dot coupled to optical and mechanical vibrational modes in a nanocavity.

Here, the electron transition energy is W, the photon and phonon mode frequencies are ω and Ω, respectively. The decay constants γ, μω, and μΩ of the electron, photon, and phonon subsystems due to couplings to their respective dissipative reservoirs are also indicated. Figure 1a sketches a single molecule in a nanocavity formed by a nanotip and a metal substrate. Here, the EM field is coupled to a transition between electron states, and this coupling is modulated by molecular vibrations. Figure 1b shows an electron transition in a QD coupled to the EM field. The figure implies that it is a QD which experiences mechanical or acoustic vibrations, but our treatment below works for any mechanism of relative displacement between the electron system and the field of an EM cavity mode, including the situations where it is the wall of a nanocavity which experiences oscillations.

We start from writing down a general Hamiltonian for a coupled quantized electron–photon–phonon system and derive its various approximate forms: the RWA, small-amplitude acoustic oscillations, classical versus quantum phonon mode, etc.

2.1 The fermion subsystem

Consider the simplest version of the fermion subsystem: two electron states |0 and |1 with energies 0 and W, respectively. We will call it an “atom” for brevity, although it can be electron states of a molecule, a QD, or any other electron system. Introduce creation and annihilation operators of the excited state |1, σˆ=|01|, σˆ=|10|, which satisfy standard commutation relations for fermions:

σˆ|0=|1,σˆ|1=|0,σˆσˆ=σˆσˆ=0;[σˆ,σˆ]+=σˆσˆ+σˆσˆ=1.

The Hamiltonian of an atom is given as follows:

(2)Hˆa=Wσˆσˆ.

We will also need the dipole moment operator,

(3)dˆ=d(σˆ+σˆ),

where d=1|dˆ|0 is a real vector. For a finite motion, we can always choose the coordinate representation of stationary states in terms of real functions.

2.2 Quantized EM modes of a cavity

We use a standard representation for the electric field operator in a cavity:

(4)Eˆ=i[Ei(r)cˆi+Ei*(r)cˆi],

where cˆi,cˆi are creation and annihilation operators for photons at frequency ωi ; the functions Ei(r) describe the spatial structure of the EM modes in a cavity. The functions Ei(r) and the relation between the modal frequency ωi  and Ei(r) can be found by solving the boundary value problem of the classical electrodynamics [23]. The normalization conditions [33]

(5)V[ωi2ε(ωi,r)]ωiωiEi*(r)Ei(r)d3r=4πωi

ensure correct bosonic commutators [cˆi,cˆi]=δij and the field Hamiltonian in the form

(6)Hˆem=iωi(cˆicˆi+12).

here, V is a quantization volume and ε(ω,r) is the dielectric function of a dispersive medium that fills the cavity. Eq. (5) is true for any fields satisfying Maxwell’s equations as long as intracavity losses can be neglected and the flux of the Poynting vector through the total cavity surface is zero (see e.g., Refs. [33], [34], [35], [36]). Of course the photon losses are always important when calculating the decoherence rates and fluctuations. What matters for Eq. (5) is that the effect of losses on the spatial structure of the cavity modes is insignificant. The latter is true as long as it makes sense to talk about cavity modes at all, which means in practice that the cavity Q-factor is at least around 10 or greater.

2.3 The quantized phonon field

We assume that our two-level atom is dressed by a phonon field which can be described by the displacement operator:

(7)qˆ=iqˆi;qˆi=Qi(r)bˆi+Qi*(r)bˆi

here, bˆi and bˆi are annihilation and creation operators of phonons, and the functions Qi(r) determine the spatial structure of oscillations at frequencies Ωi. Expression (7) can be used when the amplitude of oscillations is small enough. One can always choose the normalization of functions Qi(r) corresponding to standard commutation relations for bosons, [bˆi,bˆj]=δij and a standard form for the Hamiltonian of mechanical oscillations:

(8)Hˆp=iΩi(bˆibˆi+12).

2.4 An atom coupled to quantized EM and phonon fields

Now, we can combine all ingredients into a coupled quantized system. Adding the interaction Hamiltonian with a EM cavity mode in the electric dipole approximation, dˆEˆ, the Hamiltonian of an atom coupled to a single mode EM field can be written as follows:

(9)Hˆ=Hˆem+Hˆad(σˆ+σˆ)[E(r)cˆ+E*(r)cˆ]r=ra,

where r=ra denotes the position of an atom inside the cavity. The effect of “dressing” of the coupled atom–EM field system by mechanical oscillations in its most general form can be included by adding the Hamiltonian of phonon modes Hˆp and substituting rara+qˆ in Eq. (9). This will work for an arbitrary relative displacement of an atom with respect to the EM cavity mode. Keeping only one phonon mode for simplicity, in which

(10)qˆ=Q(r)bˆ+Q*(r)bˆ,

and expanding in Taylor series in the vicinity of r=ra, we obtain the total Hamiltonian,

(11)Hˆ=Hˆem+Hˆa+Hˆp(χσˆcˆ+χ*σˆcˆ+χσˆcˆ+χ*σˆcˆ)(η1σˆcˆbˆ+η1*σˆcˆbˆ+η2σˆcˆbˆ+η2*σˆcˆbˆ+η1σˆcˆbˆ+η1*σˆcˆbˆ+η2σˆcˆbˆ+η2*σˆcˆbˆ)

where

χ=dEr=ra,η1=dQEr=ra,η2=dQ*Er=ra.

Our model corresponds to the situation when the amplitude of phonon (acoustic) or mechanical oscillations is much larger than the size of an atom. In this case, an acoustic field shifts the potential well for electrons as a whole, rather than deforming it. It is possible to have an opposite situation when the acoustic field deforms the potential well for electrons, thus modulating the dipole moment of the optical transition. This will change the expression for the effective constant of the parametric coupling but will not change the resulting Hamiltonian.

We can always take the functions E(r) and Q(r) to be real at the position of an atom, but we cannot keep the derivatives real at the same time if the modal structure eik⋅r. However, for ideal cavity modes, the latter is possible. As we will see below, the best conditions for electron–photon–phonon entanglement are reached in the vicinity of the parametric resonance:

(12)Wω±Ω.

When the upper sign is chosen in Eq. (12), the RWA applied to the Hamiltonian (11) yields the following equation:

(13)Hˆ=Hˆem+Hˆa+Hˆp(ησˆcˆbˆ+η*σˆcˆbˆ)

where ηη1. For the lower sign in Eq. (12), the RWA Hamiltonian is as follows:

(14)Hˆ=Hˆem+Hˆa+Hˆp(ησˆcˆbˆ+η*σˆcˆbˆ)

where ηη2.

2.5 An atom coupled to the quantized EM field and dressed by a classical acoustic field

For classical acoustic oscillations, the operator qˆ=Q(r)bˆ+Q*(r)bˆ in Eq. (10) becomes a classical function

(15)q=Q(r)eiΩt+Q*(r)eiΩt

where Q is a coordinate-dependent complex amplitude of classical oscillations. Near the parametric resonance (ω+ΩW), the RWA Hamiltonian takes the following form:

(16)Hˆ=Hˆem+Hˆa(σˆcˆeiΩt+*σˆcˆeiΩt).

where =[d(Q)E]r=ra. The value of the acoustic frequency Ω in Eq. (16) can be of either sign, corresponding to the choice “±” in the parametric resonance condition Eq. (12); when the sign of Ω changes from positive to negative, one should replace Q with Q* in the above expression for .

Qualitatively, Hamiltonian (13) corresponds to the decay of the fermionic excitation into a photon and phonon; Hamiltonian (14) corresponds to the decay of a photon into a phonon and fermionic excitation, whereas Hamiltonian (16) describes parametric decay of a photon into an atomic excitation and back, mediated by classical acoustic oscillations.

3 Parametric resonance in a closed system

When the system is closed and there is no dissipation, the general analytic solution to the dynamics of coupled fermions, photons, and phonons can be obtained in the RWA. We write the state vector as follows:

(17)Ψ=α,n=0(Cαn0|α|n|0+Cαn1|α|n|1).

here, Greek letters denote phonon states, Latin letters denote photon states, and numbers 0, 1 describe fermion states. We will keep the same sequence of symbols throughout the paper:

Cphononphotonfermion|phonon|photon|fermion.

Next, we substitute Eq. (17) into the Schrödinger equation,

(18)it|Ψ=Hˆ|Ψ

where Hˆ is the RWA Hamiltonian. For definiteness, we consider the vicinity of the parametric resonance with a plus sign, ω+ΩW, which corresponds to the Hamiltonian (13). In this case, the equations for the coefficients in Eq. (17) can be separated into the pairs of coupled equations

(19)ddt(Cαn0C(α1)(n1)1)+(iωα,niΩR(α,n)*iΩR(α,n)iωα,niΔ)(Cαn0C(α1)(n1)1)=0,

and a separate equation for the lowest energy state:

(20)C000+iω0,0C000=0,

where

ΩR(α,n)=ηαn,ωα,n=Ω(α+12)+ω(n+12),Δ=Ω+ωW.

Note that approximate Eqs. (19) and (20) preserve the norm exactly:

|C000|2+α=1,n=1,(|Cαn0|2+|C(α1)(n1)1|2)=α=0,n=0,(|Cαn0|2+|Cαn1|2)=const.

The solution to Eq. (20) is trivial: C000(t)=C000(0)exp(iω0,0t). The solution to Eq. (19) takes the following form:

(21)(Cαn0C(α1)(n1)1)=AeΛ1(α,n)t(1a1(α,n))+BeΛ2(α,n)t(1a2(α,n)),

where the constants A and B are determined from initial conditions. Here, the eigenvalues Λ1,2(α,n) and eigenvectors (1a1,2(α,n)) of the matrix of coefficients in Eq. (19) are given as follows:

(22)Λ1,2(α,n)=iωα,niδ1,2(α,n),a1,2(α,n)=δ1,2(α,n)ΩR(α,n)*,

where

(23)δ1,2(α,n)=Δ2±Δ24+|ΩR(α,n)|2.

Figure 2 shows the eigenfrequencies of the system given by Eq. (22) with α=n=1, shifted by ω1,1|Δ=0. One can see the anticrossing with splitting by 2ΩR(1,1) at the parametric resonance.

Figure 2: Frequency eigenvalues of the coupled electron–photon–phonon quantum system as a function of detuning from the parametric resonance Wℏ=Ω+ω$\frac{W}{\hslash }={\Omega}+\omega $. All frequencies are in units of the generalized Rabi frequency ΩR(1,1)${{\Omega}}_{R}^{\left(1,1\right)}$. The values of eigenfrequencies are shifted vertically by ω1,1|Δ=0${\omega }_{1,1}{\vert }_{{\Delta}=0}$.
Figure 2:

Frequency eigenvalues of the coupled electron–photon–phonon quantum system as a function of detuning from the parametric resonance W=Ω+ω. All frequencies are in units of the generalized Rabi frequency ΩR(1,1). The values of eigenfrequencies are shifted vertically by ω1,1|Δ=0.

As an example, consider an exact parametric resonance W=Ω+ω and the simplest initial state Ψ0=|0|0|1 corresponding to the initially excited atom in a cavity. In this case, the only nonzero amplitudes are C001 and C110:

(24)(C110C001)=12ei(ω1,1|ΩR(1,1)|)t(eiθ1)+12ei(ω1,1+|ΩR(1,1)|)t(eiθ1),

where

ω1,1=Ω(1+12)+ω(1+12),ΩR(1,1)=η=|ΩR(1,1)|eiθ.

The resulting state vector is given as follows:

(25)Ψ=eiω1,1t[ieiθsin(|ΩR(1,1)|t)|1|1|0+cos(|ΩR(1,1)|t)|0|0|1].

This is clearly an entangled electron–photon–phonon state, which is not surprising. In the absence of dissipation, any coupling between these subsystems leads to entanglement.

State (25) is a tripartite entangled state which belongs to the family of Greenberger–Horne–Zeilinger (GHZ) states. It can be reduced to a standard GHZ state by local operations [37], [38], e.g., by rotations on the Bloch sphere of each qubit. In most cases discussed in the literature, the GHZ states are made of identical subsystems, e.g., photons [39], [40], which determine the way how they can be manipulated and used. In our case, each subsystem is of different nature: a fermionic electron system, a bosonic EM field, and a bosonic phonon field, so we envision at least two interesting applications. One is to determine the statistics of phonons or atomic excitations by measuring the statistics of photons. The latter is relatively easy to do, whereas phonon counting or direct measurement of quantized mechanical oscillations is extremely difficult due to their low energy and the lack of suitable detectors. The second application is control of a bipartite entangled state of two subsystems, say an atom and a photon mode, by using the state of the third subsystem, say phonons or mechanical oscillations, as a control handle. One example is given in Section 7 of the paper. One can come up with other combinations involving more complex qubits consisting of coupled fermion–boson subsystems, for example, an entanglement of the atom–phonon system by a classical EM field, etc.

The dynamics of the corresponding physical observables, such as the energy of the field and the atom, is Rabi oscillations at the frequency which generalizes a standard Rabi frequency to the case of a parametric photon–phonon–atom resonance and which depends on both the spatial structure of the photon and phonon fields and their occupation numbers:

(26)Ψ|Eˆ2|Ψ=|E(r)|2[2cos(2|ΩR(1,1)|t)]
(27)Ψ|Hˆa|Ψ=W1+cos(2|ΩR(1,1)|t)2

It is illustrated in Figure 3 which shows the normalized EM field energy density and energy of an atom as a function of time. Note that the EM field energy never reaches zero because of the presence of zero-point vacuum energy. With detuning from the parametric resonance, the amplitude of the oscillations will decrease.

Figure 3: (a) Normalized field intensity, 〈Ψ|Eˆ2|Ψ〉/|E(r)|2$\langle {\Psi}\vert {\hat{\mathbf{E}}}^{2}\vert {\Psi}\rangle /{\vert \mathbf{E}\left(\mathbf{r}\right)\vert }^{2}$, and (b) normalized atom energy 〈Ψ|Hˆa|Ψ〉/W$\langle {\Psi}\vert {\hat{H}}_{a}\vert {\Psi}\rangle /W$ as a function of time in units of the generalized Rabi frequency ΩR(1,1)${{\Omega}}_{R}^{\left(1,1\right)}$.
Figure 3:

(a) Normalized field intensity, Ψ|Eˆ2|Ψ/|E(r)|2, and (b) normalized atom energy Ψ|Hˆa|Ψ/W as a function of time in units of the generalized Rabi frequency ΩR(1,1).

4 Dynamics of an open electron–photon–phonon system

4.1 Stochastic evolution equation

Now we include the processes of relaxation and decoherence in an open system, which is (weakly) coupled to a dissipative reservoir. We will use the stochastic equation of evolution for the state vector, which is derived in Appendix. This is basically the Schrödinger equation modified by adding a linear relaxation operator and the noise source term with appropriate correlation properties. The latter are related to the parameters of the relaxation operator, which is a manifestation of the fluctuation–dissipation theorem [41].

Within our approach, the system is described by a state vector which has a fluctuating component: |Ψ=|Ψ¯+|Ψ˜, where the straight bar means averaging over the statistics of noise, and the wavy bar denotes the fluctuating component. This state vector is of course very different from the state vector obtained by solving a standard Schrödinger equation for a closed system. In fact, coupling to a dissipative reservoir leads to the formation of a mixed state, which can be described by a density matrix ρˆ=|Ψ¯Ψ|¯+|Ψ˜Ψ|˜¯. In Appendix, we derived the general form of the stochastic equation of evolution from the Heisenberg–Langevin equations [23], [34], [42] and showed how physically reasonable constraints on the observables determine the properties of the noise sources. We also demonstrated the relationship between our approach and the Lindblad method of solving the master equation.

One can view the stochastic equation approach as a convenient formalism for calculating physical observables which allows one to obtain analytic solutions for the evolution of a coupled system even in the presence of dissipation and decoherence. In this section, we use the stochastic equation for the state vector given by Eqs. (A9) and (A10). The effective Hamiltonian in Eqs. (A9) and (A10) is determined by an approximation the user wants. If one wants the Markovian approximation, the Hamiltonian is obtained simply by summing up partial Lindbladians for all subsystems, whatever they are (in our case, these are a fermion emitter, an EM cavity mode, and a phonon mode). Then, the noise source term is determined unambiguously by conservation of the norm of the state vector and the requirement that the system should approach thermal equilibrium when the external perturbation is turned off. This immediately gives Eqs. (28) and (29) below.

Following the derivation in Appendix, Eqs. (19) and (20) are modified as follows:

(28)ddt(Cαn0C(α1)(n1)1)+(iωα,n+γαn0iΩR(α,n)*iΩR(α,n)iωα,niΔ+γ(α1)(n1)1)(Cαn0C(α1)(n1)1)=i(Rαn0R(α1)(n1)1),
(29)C˙000+(iω0,0+γ000)C000=iR000.

Coupling to a reservoir introduces two main differences to Eqs. (28) and (29) as compared to Eqs. (19) and (20) for a closed system. First, eigenfrequencies acquire imaginary parts which describe relaxation:

ωα,nωα,niγαn0,ωα,nΔωα,nΔiγ(α1)(n1)1,ω0,0ω0,0iγ000.

The relaxation constants are determined by the properties of all subsystems. They are derived in Appendix, and their explicit form is given in the end of this section.

Second, the right-hand side of Eqs. (28) and (29) contain noise sources iRαn0, iR(α1)(n1)1, and iR000. They are equal to 0 after averaging over the noise statistics: Rαn0¯=R(α1)(n1)1¯=R000¯. The averages of the quadratic combinations of noise source terms are nonzero. Including the noise sources is crucial for consistency of the formalism: it ensures the conservation of the norm of the state vector and leads to a physically meaningful equilibrium state. Note that the Weisskopf–Wigner theory does not enforce the conservation of the norm.

4.2 Evolution of the state amplitudes and observables

The solution to Eq. (29) is given as follows:

(30)C000=e(iω0,0+γ000)t[C000(0)i0te(iω0,0+γ000)τR000(τ)dτ].

The solution to Eq. (28) is determined again by the eigenvalues and eigenvectors of the matrix of coefficients, which are now modified by relaxation rates:

(31)Λ1,2(α,n)=iωα,niδ1,2(α,n),a1,2(α,n)=δ1,2(α,n)iγαn0ΩR(α,n)*,

where

(32)δ1,2(α,n)=Δ2+iγαn0+γ(α1)(n1)12±[Δ+i(γ(α1)(n1)1γαn0)]24+|ΩR(α,n)|2.

The solution to Eq. (28) takes the following form:

(33)(Cαn0C(α1)(n1)1)=eΛ1(α,n)t(1a1(α,n))(Ai0teΛ1(α,n)τRαn0(τ)a2(α,n)R(α1)(n1)1(τ)a2(α,n)a1(α,n)dτ)+eΛ2(α,n)t(1a2(α,n))(Bi0teΛ2(α,n)τR(α1)(n1)1(τ)Rαn0(τ)a1(α,n)a2(α,n)a1(α,n)dτ)

where the constants A and B are determined by initial conditions.

As an example, we consider the reservoir at low temperatures, when the steady-state population should go to the ground state |0|0|0. Also, we will neglect purely elastic dephasing processes which lead to atomic decoherence without changing the populations. The elastic processes will be added later. In this case, we can take γ000=0, as shown below. We will also assume that the only nonzero correlator of noise is delta-correlated in time:

(34)R000(t+ξ)R000*(t)¯=2δ(ξ)D000.

Then, Eqs. (29) and (30) yield the following:

(35)ddt|C000|2¯=D000,

whereas Eq. (28) gives the following:

(36)ddt(|Cαn0|2¯+|C(α1)(n1)1|2¯)=2(γαn0|Cαn0|2¯+γ(α1)(n1)1|C(α1)(n1)1|2¯).

This equation guarantees that the system occupies the ground state at t.

The noise intensity is determined by the condition that the norm of the state vector be conserved. This gives the following equation:

(37)D000=2α=1,n=1,(γαn0|Cαn0|2¯+γ(α1)(n1)1|C(α1)(n1)1|2¯).

In Appendix, we discuss in detail the dependence of the noise correlator on the averaged dyadic components of the state vector. We also show how to find the correlators which ensure that the system approaches thermal distribution at a finite temperature.

The above formalism allows us to obtain analytic solutions to the state vector and observables at any temperatures and detunings from the parametric resonance, while still within the RWA limits. However, the resulting expressions are very cumbersome, and they are better to visualize in the plots. Let us give an example of the solution at zero reservoir temperature and exactly at the parametric resonance W=Ω+ω, when the expressions are more manageable. Consider the initial state Ψ0=|0|0|1 when an atom is excited and boson modes are in the ground state. In this case, the only nonzero amplitudes are C000, C001, and C110. To make the algebra a bit simpler, we assume that the dissipation is weak enough and its effect on the eigenvectors (1a1,2(α,n)) can be neglected. As a result, we obtain the following equation:

(38)Ψ=e(iω1,1+γ001+γ1102)t[ieiθsin(|Ω˜R(1,1)|t)|1|1|0+cos(|Ω˜R(1,1)|t)|0|0|1]+C000|0|0|0,

where

|C000|2¯=1e(γ110+γ001)t,Ω˜R(1,1)=|ΩR(1,1)|2(γ001γ110)24,θ=Arg[ΩR(1,1)].

As we see, dissipation leads not only to the relaxation of the entangled part of the state vector but also to the frequency shift of the Rabi oscillations. This shift is absent if γ001=γ110.

The resulting expressions for the observables such as the EM field intensity and the energy of the atomic excitation are given as follows:

(39)Ψ|Eˆ2|Ψ=|E(r)|2[1+e(γ110+γ001)tcos(2|Ω˜R(1,1)|t)e(γ110+γ001)t],
(40)Ψ|Hˆa|Ψ=W1+cos(2|Ω˜R(1,1)|t)2e(γ110+γ001)t

Figure 4 illustrates the dynamics of observables in Eqs. (39) and (40).

Figure 4: (a) Normalized field intensity, 〈Ψ|Eˆ2|Ψ〉/|E(r)|2$\langle {\Psi}\vert {\hat{\mathbf{E}}}^{2}\vert {\Psi}\rangle /{\vert \mathbf{E}\left(\mathbf{r}\right)\vert }^{2}$, and (b) normalized atom energy 〈Ψ|Hˆa|Ψ〉/W$\langle {\Psi}\vert {\hat{H}}_{a}\vert {\Psi}\rangle /W$ a function of time in units of the generalized Rabi frequency ΩR(1,1)${{\Omega}}_{R}^{\left(1,1\right)}$. Here, γ110+γ001=0.3ΩR(1,1)${\gamma }_{110}+{\gamma }_{001}=0.3{{\Omega}}_{R}^{\left(1,1\right)}$.
Figure 4:

(a) Normalized field intensity, Ψ|Eˆ2|Ψ/|E(r)|2, and (b) normalized atom energy Ψ|Hˆa|Ψ/W a function of time in units of the generalized Rabi frequency ΩR(1,1). Here, γ110+γ001=0.3ΩR(1,1).

Note that the Weisskopf–Wigner theory would give the same expression (40) for the atomic energy but a wrong expression for the EM field intensity:

Ψ|Eˆ2|Ψ=|E(r)|2[2cos(2|Ω˜R(1,1)|t)]e(γ110+γ001)t,

which does not approach the correct vacuum state.

4.3 Emission spectra

According to the study by Scully and Zubairy [23], the power spectrum of the emission is given as follows:

(41)S(r,ν)=1πRe0dτG(1)(r,r;τ)eiντ,

where G(1)(r,r;τ) is the field autocorrelation function at the position r of the detector:

(42)G(1)(r,r;τ)=|E(r)|20dtcˆd(t)cˆd(t+τ)¯.

where cˆd(t),cˆd(t) are annihilation and creation operators for the photons which interact with the detector, and the Heisenberg picture is used. We will assume that the coupling between the photons and the detector is weak, so the photon detection does not affect the dynamics of the intracavity photons. According to the study by Madsen and Lodahl [45], cˆd(t)cˆ(t) for a nanocavity, so we can calculate the G(1)(r,r;τ) using operators for the cavity field cˆ(t),cˆ(t), up to a constant factor in the result. Note that the lower limit of the integral over t is set to be t=0, which requires that no photons exist before t=0.

In the Heisenberg–Langevin approach, an operator in the Heisenberg picture can be expressed through Schrödinger’s operators using the effective Hamiltonian Hˆeff, which contains the anti-Hermitian part, see the Appendix. At the same time, the inhomogeneous term proportional to the noise sources should be added. Including these noise terms in the solution for the field operators when calculating the emission spectra is equivalent to taking into account the detection of thermal radiation which seeps into the cavity from outside and spontaneous emission resulting from thermal excitation of an atom. We assume that the reservoir temperature in energy units is much lower than W and ω so that the contribution of these noise terms to the emission spectra can be neglected (although noise is still needed to preserve the norm).

Then, the average correlator cˆ(t)cˆ(t+τ)¯ is expressed as follows:

(43)cˆ(t)cˆ(t+τ)¯=Ψ(t=0)|eiHˆefft/cˆeiHˆefft/eiHˆeff(t+τ)/cˆeiHˆeff(t+τ)/|Ψ(t=0)¯=Ψ(t)|cˆeiHˆefft/eiHˆeff(t+τ)/cˆ|Ψ(t+τ)¯,

where |Ψ(t) is the state vector of the system which we found in the previous subsection. It can be written as |Ψ(t)=n=0Cn(t)|n|Ψnα,e(t), where |Ψnα,e(t) is the part describing phonons and electrons. Therefore,

(44)cˆ(t)cˆ(t+τ)=(n=0Cn*(t)n|Ψnα,e(t)|)cˆeiHˆefft/eiHˆeff(t+τ)/cˆ(n=0Cn(t+τ)|n|Ψnα,e(t+τ))=(n=0nCn*(t)n1|Ψnα,e(t)|)eiHˆefft/eiHˆeff(t+τ)/(n=0nCn(t+τ)|n1|Ψnα,e(t+τ)).

Consider a simple example when the initial state is |0|0|1. Within the RWA, the system can only reach states |0|0|1, |1|1|0, and |0|0|0. After acting with cˆ on a state of the system, a new state |1|0|0 can also appear, but it cannot evolve into other states. So, in this case, we have the following:

(45)cˆ(t)cˆ(t+τ)¯=(C1*(t)0|Ψ1α,e(t)|)eiHˆefft/eiHˆeff(t+τ)/(C1(t+τ)|0|Ψ1α,e(t+τ))¯=(C110*(t)1|0|0|)eiHˆefft/eiHˆeff(t+τ)/(C110(t+τ)|1|0|0)¯=C110*(t)C110(t+τ)exp[iω1,0τγ100(2t+τ)],

where we used Eqs. (A32) and (A33) and assumed that the noise for state |1|0|0 has zero correlator. Since

(46)C110(t)=isin(|Ω˜R(1,1)|t)exp[iω1,1tγ110+γ0012t],

we obtain

(47)cˆ(t)cˆ(t+τ)¯=sin(|Ω˜R(1,1)|t)sin(|Ω˜R(1,1)|(t+τ))exp[iωτ]exp[γac(2t+τ)],

where we introduced the notation γacγ100+γ110+γ0012. Then, the power spectrum is found to be

(48)S(r,ν)1π|E(r)|2|Ω˜R(1,1)|24γac(|Ω˜R(1,1)|2+γac2)Re[2γaci(νω)[γaci(νω)]2+|Ω˜R(1,1)|2].

The normalized power spectra are shown in Figure 5 for various values of |Ω˜R(1,1)|/γac. For |Ω˜R(1,1)|<γac, the spectrum has a single maximum at zero detuning ν=ω. For |Ω˜R(1,1)|>γac, the spectra are split and their maxima (same value for all spectra) are reached at detunings given by (νω)2=|Ω˜R(1,1)|2γac2. Therefore, to reach the strong coupling regime, the Rabi frequency |Ω˜R(1,1)| has to exceed the combination of the decoherence rates denoted by γac.

Figure 5: The emission spectra for |Ω˜R(1,1)|/γac$\vert {\tilde {{\Omega}}}_{R}^{\left(1,1\right)}\vert /{\gamma }_{\text{ac}}$ equal to 0.5, 1, 2, and 5. All spectra are normalized by the same constant.
Figure 5:

The emission spectra for |Ω˜R(1,1)|/γac equal to 0.5, 1, 2, and 5. All spectra are normalized by the same constant.

Note that in a standard Lindblad formalism, γ parameters are inverse decay times for the populations or field energies, i.e., they are related to the quantities which are quadratic with respect to the state vector within the Schrödinger’s approach. That is why our γ’s in the stochastic equation of evolution, which is linear with respect to the state vector, correspond to the Lindblad γ’s divided by 2.

4.4 Relaxation rates

Finally, we give explicit expressions for the relaxation constants γαn0 and γαn1. They were derived in Appendix using the Lindblad master equation approach and assuming statistical independence of “partial” dissipative reservoirs for the atomic, EM, and phonon subsystems. Within the model which neglects purely elastic dephasing processes, the result is as follows:

(49)γαn0=γ2N1Ta+μω2[n¯ωTem(n+1)+(n¯ωTem+1)n]+μΩ2[n¯ΩTp(α+1)+(n¯ΩTp+1)α],
(50)γαn1=γ2N0Ta+μω2[n¯ωTem(n+1)+(n¯ωTem+1)n]+μΩ2[n¯ΩTp(α+1)+(n¯ΩTp+1)α],

where γ, μω, and μΩ are partial relaxation rates of the atomic, photon, and phonon subsystems, respectively; N0Ta=11+eWTa, N1Ta=eWTa1+eWTa, n¯ωTem=1eωTem1, n¯ΩTp=1eΩTp1 are their occupation numbers at thermal equilibrium; Ta,em,p are temperatures of partial atom, photon, and phonon reservoirs in energy units, respectively. As a reminder, the atom energy is equal to 0 in state |0 and W in state |1.

If all reservoirs are at zero temperature, we obtain the following:

(51)γαn0=μω2n+μΩ2α,γαn1=γ2+μω2n+μΩ2α.
Eq. (51) shows that γ000=0, validating our choice earlier in this section. We also obtain physically intuitive expressions for γ110 and γ001: γ110=μω2+μΩ2, γ001=γ2.

Purely elastic processes which lead to the atomic decoherence with characteristic time 1/γel can be taken into account by adding γel to the relaxation constant γ001. At the same time, one has to modify the noise correlator R001(t+ξ)R001*(t)¯ by adding to it the quantity 2γel2|C001|2¯. This prescription will lead to correct dynamics of the observables; see the last section in the Appendix. Note that the population relaxation times will not depend on γel as it should be.

The above expressions allow one to quickly estimate the feasibility of reaching strong coupling regime and quantum entanglement when all fields are quantized. In a semiconductor dielectric cavity at near-infrared wavelengths ∼1 μm for the refractive index ∼3.5 and a typical dipole matrix element of the interband optical transition d/e∼0.5 nm, the maximum vacuum Rabi frequency is of the order of 100–200 μeV [13], [14], diffraction limited by the cavity size. This sets the upper limit for the sum of relaxation rates γacγ100+γ110+γ0012 introduced above. In the experiments with single QDs [13], [14], the phonons were not involved and the sum of relaxation rates in the electron and cavity subsystems was kept below 100 μeV at low temperatures, allowing them to reach the strong coupling regime.

In plasmonic cavities, a sub-nm field localization can be achieved, leading to Rabi frequencies of the order of 100–200 meV for the same order of the transition dipole moments. However, the combined relaxation rate is much higher, up to ∼100 meV, typically dominated by cavity losses. Here, the strong coupling of plasmons to a single quantum emitter has been achieved in multiple experiments, as discussed in the Section 1. However, reaching the quantum regime for the EM and especially the phonon fields remains a challenge. It could be beneficial to consider longer wavelength emitters with the optical transition at the midinfrared and even terahertz wavelengths. Indeed, with increasing wavelength, the plasmon losses go down and the matrix element of a dipole-allowed transition increases, whereas the plasmon localization stays largely the same.

5 Classical acoustic or mechanical oscillations

Quantization of acoustic or mechanical oscillations is very difficult to achieve because of their low energy of the quantum. In most experiments, they remain classical and therefore cease to be an independent degree of freedom. Instead, their amplitude becomes an external time-dependent parameter, like an external pumping. In this case, the RWA Hamiltonian is given by Eq. (16). It depends only on quantum operators σˆ,σˆ and cˆ,cˆ; therefore, the state vector has to be expanded over the basis states |n|0 and |n|1:

(52)Ψ=n=0(Cn0|n|0+Cn1|n|1).

Substituting Eq. (52) in the Schrödinger equation with the Hamiltonian (16), we again get separation into a block of two equations,

(53)C˙n0=iωnCn0+i*eiΩtC(n1)1n,
(54)C˙(n1)1=i(ωn1+W)C(n1)1+ieiΩtCn0n,

and a separate equation for the amplitude of the ground state |0|0 of the system:

(55)C˙00=iω0C00,

where ωn=ω(n+12) and =[d(Q)E]r=ra (see Eq. (16)). After making the substitution C(n1)1=G(n1)1eiΩt, Eqs. (53) and (54) give the equations similar in form to Eq. (19):

(56)ddt(Cn0G(n1)1)+(iωniΩR(n)*iΩR(n)iωniΔ)(Cn0G(n1)1)=0,

where

ΩR(n)=n,Δ=Ω+ωW,ωnΔ=ωn1+W.
Eqs. (55) and (56) are different from Eqs. (19) and (20) only in one aspect: they do not contain the index of the quantum state of the phonon field, whereas the Rabi frequency depends on the amplitude of classical acoustic oscillations Q(ra), see Section 2.5. Obviously, the solution to Eqs. (55) and (56) will have the same form and the expressions (26), (27) for the observables will remain the same, after dropping the index of the quantum phonon state and redefining the Rabi frequency.

Dissipation due to coupling to a reservoir can be included using the stochastic equation of evolution of the state vector, see the Appendix. The corresponding equations are again similar to those for a fully quantum problem given by Eqs. (28) and (29):

(57)C˙00+i(ω0+γ00)C00=iR00,
(58)ddt(Cn0C(n1)1)+(iωn+γn0iΩR(n)*iΩR(n)iωniΔ+γ(n1)1)(Cn0C(n1)1)=i(Rn0R(n1)1).

Since the acoustic field is now a given external pumping, the relaxation constants should not depend on the parameters of a phonon reservoir. They can be obtained after obvious simplification of Eqs. (49) and (50):

(59)γn0=γ2N1Ta+μω2[n¯ωTem(n+1)+(n¯ωTem+1)n],
(60)γn1=γ2N0Ta+μω2[n¯ωTem(n+1)+(n¯ωTem+1)n],

All expressions for the state vector and observables can be obtained from the corresponding expressions in Section 4 after dropping the index α of the quantum state of the phonon field and redefining the frequency of Rabi oscillations.

6 Separation and interplay of the parametric and one-photon resonance

For an electron system coupled to a EM cavity mode and dressed by a phonon field, the phonon frequency Ω can be much lower than the optical frequency. In this case, the overlap of the parametric (three-wave) resonance ω±ΩW and the one-photon (two-wave) resonance ωW can be an issue.

Here, we derive the analytic criteria for the separation of these resonances and show numerically what happens when they overlap. Throughout this section, we neglect losses. Since the parametric and one-photon resonances are separated by the phonon frequency Ω, the losses can be neglected when the value of Ω exceeds the sum of the spectral widths of the EM cavity mode and the electron transition which originate from dissipation and decoherence. If this condition is violated, resonances overlap strongly and their separation is impossible anyway.

The separation criterion imposes certain restrictions on the Rabi frequencies of the two resonances. To derive these restrictions, we neglect dissipation and retain in the Hamiltonian (11) both the RWA terms near the parametric resonance ω±ΩW and the terms near a one-photon resonance ωW. Since the result will be almost the same whether the phonon field is quantized or classical, we will consider the classical phonon field to keep the expressions a bit shorter. The resulting Hamiltonian is given as follows:

(61)Hˆ=ω(cˆcˆ+12)+Wσˆσˆ(χ+eiΩt)σˆcˆ(χ*+*eiΩt)σˆcˆ,

where χ=(dE)r=ra and was defined in Eq. (16); Q is now a complex-valued amplitude of classical phonon oscillations. The value of Ω in Eq. (61) can be both positive and negative, corresponding to the choice of an upper or lower sign in the parametric resonance condition ω±ΩW. The change of sign in Ω corresponds to replacing Q with Q* in the expression for .

The state vector should be sought in the form of Eq. (52). After substituting it into the Schrödinger equation, we obtain coupled equations for the amplitudes of basis states |n|0, |n1|1:

(62)C˙n0+iωnCn0i(χ*+*eiΩt)nC(n1)1=0,
(63)C˙(n1)1+i(ωn1+W)C(n1)1i(χ+eiΩt)nCn0=0,

and

(64)C˙00+iω0C00=0,

where ωn=ω(n+12). To compare these equations with Eqs. (53) and (54), it is convenient to assume that the system is exactly at one of the resonances and study the behavior of the solution with increasing the detuning from another resonance. For example, we assume an exact parametric resonance ω+Ω=W. In this case, the detuning from the two-wave resonance is Wω=Ω. After the substitution Cn0=Gn0eiωnt and C(n1)1=G(n1)1ei(ωn1+W)t, we obtain from Eqs. (62) and (63) that

(65)G˙n0i(χ*+*eiΩt)nG(n1)1ei(Wω)t=0,
(66)G˙(n1)1i(χ+eiΩt)nGn0ei(Wω)t=0.

If we neglect at first the perturbation of the system in the vicinity of the two-wave resonance, the solution to Eqs. (65) and (66) at χ=0 is

(67)(Gn0G(n1)1)=AeiΩR(3)t(11)+BeiΩR(3)t(11),

where ΩR(3)=1n is the Rabi frequency of the parametric resonance and A and B are arbitrary constants. The state described by Eq. (67) is obviously entangled.

To write the formal solution to Eqs. (65) and (66), we make another substitution of variables: Gn0±G(n1)1=G±. The result is

G˙±iΩR(3)G±=iΩR(2)*eiΩtGn0±iΩR(2)eiΩtG(n1)1,

where ΩR(2)=1χn is the Rabi frequency corresponding to the one-photon (two-wave) resonance. The solution to the last equation is given as follows:

(68)G±=2(A,B)e±iΩR(3)t+ie±iΩR(3)t0teiΩR(3)τ[ΩR(2)*eiΩτGn0(τ)±ΩR(2)eiΩτG(n1)1(τ)]dτ.

Considering the terms proportional to ΩR(2) as perturbation, we seek the solution as

(Gn0G(n1)1)=AeiΩR(3)t(11)+BeiΩR(3)t(11)+(δGn0δG(n1)1).

To estimate the magnitude of the perturbation, we substitute Eq. (67) into Eq. (68). After some algebra, we obtain that under the condition ΩR(3)Ω, the magnitude of the perturbation is given as follows:

δGn0,(n1)1|ΩR(2)Ω|Gn0,(n1)1,

whereas if ΩR(3)Ω, the magnitude of the perturbation is given as follows:

δGn0,(n1)1|ΩR(2)ΩR(3)|Gn0,(n1)1.

To summarize this part, if both Rabi frequencies ΩR(3), ΩR(2)Ω, the two resonances can be treated independently for any relationship between the magnitudes of ΩR(3) and ΩR(2). If the above inequality is violated, one can neglect one of the resonances only if its associated Rabi frequency is much lower than the Rabi frequency of another resonance. These restrictions are obvious from qualitative physical reasoning: either the magnitudes of the Rabi splittings are much smaller than the distance between resonances or one of the splittings is much weaker than another one.

When the effect of the neighboring resonance is nonnegligible, it can still be taken into account in the solution. Indeed, consider the solution to Eqs. (62) and (63), taking into account only the two-wave resonance, i.e., taking =0. After obvious substitutions, we arrive at

(69)(Cn0C(n1)1)=Aei(ωnΩ2+Ω24+|ΩR(2)|2)t×(1Ω2+Ω24+|ΩR(2)|2ΩR(2))+Bei(ωn1+W+Ω2Ω24+|ΩR(2)|2)t×(ΩR(2)Ω2Ω24+|ΩR(2)|21);

In the limit ΩΩR(2), we obtain the following:

(70)(Cn0C(n1)1)Aei(ωn+|ΩR(2)|2Ω)t(1|ΩR(2)|Ω)+Bei(ωn1+W|ΩR(2)|2Ω)t(|ΩR(2)|Ω1)

It is clear from Eq. (70) that the entanglement of states described by Cn0 and C(n1)1 is determined by a small parameter |ΩR(2)|Ω, whereas at exact resonance, the entanglement is always stronger, see Eq. (67). Therefore, when ΩR(2)Ω, we can neglect the contribution of the two-photon resonance to the entanglement of states |n|0 and |n1|1. However, it follows from Eq. (69) that the two-wave resonance shifts the eigenfrequencies of the system. Qualitatively, these shifts can be included by putting χ=0 in Eqs. (62) and (63) but replacing the eigenfrequencies ωn and ωn1 according to Eq. (70):

(71)ωnωn+|ΩR(2)|2Ω,ωn1+Wωn1+W|ΩR(2)|2Ω.

If ΩR(3)|ΩR(2)|2Ω, these shifts can be significant in order to interpret the spectra near the three-wave parametric resonance.

The same reasoning can be carried out to analyze the effect of a detuned three-wave resonance on the solution near the two-wave resonance.

These results can be verified by an exact numerical solution of Eqs. (62) and (63) for given initial conditions. After that, we can obtain the spectra of Cn0 and C(n1)1. Since they are oscillating functions, their spectra form discrete lines at frequencies which we denote as ωosc.

As an example, we select the case of n=1, set |ΩR(2)|=|ΩR(3)|=0.1Ω, and choose the initial condition as Cn0(0)=0 and C(n1)1(0)=1. The frequencies ωosc of the spectral lines for Cn0 and C(n1)1 are shown in Figure 6. Their values are shifted by ωosc,0=ω1|ω=W/. The area of the dot for each spectral line is proportional to the square of its amplitude. If a marker is not visible, it means the corresponding line is very weak and can be neglected. The anticrossing can be seen at both the one-photon resonance and parametric resonance.

Figure 6: The frequencies ωosc${\omega }_{\text{osc}}$ of the spectral lines for Cn0${C}_{n0}$ (left panel) and C(n−1)1${C}_{\left(n-1\right)1}$ (right panel), with n=1$n=1$, as functions of the photon frequency ω. The photon frequencies are shifted by W/ℏ$W/\hslash $, and the positions of spectral lines ωosc${\omega }_{\text{osc}}$ are shifted by ωosc,0=ω1|ω=W/ℏ${\omega }_{osc,0}={\omega }_{1}{\vert }_{\omega =W/\hslash }$. The area of a marker is proportional to the amplitude squared of the spectral line. Both axes are in units of Ω${\Omega}$. The parameters are |ΩR(2)|=|ΩR(3)|=0.1 Ω$\vert {{\Omega}}_{R}^{\left(2\right)}\vert =\vert {{\Omega}}_{R}^{\left(3\right)}\vert =0.1\,{\Omega}$, and the initial condition is Cn0(0)=0${C}_{n0}\left(0\right)=0$ and C(n−1)1(0)=1${C}_{\left(n-1\right)1}\left(0\right)=1$.
Figure 6:

The frequencies ωosc of the spectral lines for Cn0 (left panel) and C(n1)1 (right panel), with n=1, as functions of the photon frequency ω. The photon frequencies are shifted by W/, and the positions of spectral lines ωosc are shifted by ωosc,0=ω1|ω=W/. The area of a marker is proportional to the amplitude squared of the spectral line. Both axes are in units of Ω. The parameters are |ΩR(2)|=|ΩR(3)|=0.1Ω, and the initial condition is Cn0(0)=0 and C(n1)1(0)=1.

As an illustration of the violation of the condition for resonance separation, we show the oscillation frequencies for |ΩR(2)|=|ΩR(3)|=0.5Ω in Figure 7. Here, the anticrossing picture of isolated resonances is smeared and cannot be observed.

Figure 7: The frequencies ωosc${\omega }_{\text{osc}}$ of the spectral lines for Cn0${C}_{n0}$ (left panel) and C(n−1)1${C}_{\left(n-1\right)1}$ (right panel), with n=1$n=1$. The notations are the same as in Figure 6. The parameters are |ΩR(2)|=|ΩR(3)|=0.5 Ω$\vert {{\Omega}}_{R}^{\left(2\right)}\vert =\vert {{\Omega}}_{R}^{\left(3\right)}\vert =0.5\,{\Omega}$, and the initial condition is Cn0(0)=0${C}_{n0}\left(0\right)=0$ and C(n−1)1(0)=1${C}_{\left(n-1\right)1}\left(0\right)=1$.
Figure 7:

The frequencies ωosc of the spectral lines for Cn0 (left panel) and C(n1)1 (right panel), with n=1. The notations are the same as in Figure 6. The parameters are |ΩR(2)|=|ΩR(3)|=0.5Ω, and the initial condition is Cn0(0)=0 and C(n1)1(0)=1.

7 Control of entangled states

In order to control the quantum state of the system, turn the entanglement on/off, read or write information into a qubit, or implement a logic gate, one has to vary the parameters of a system, for example, the detuning from resonance, the field amplitude of the EM mode at the atom position, or the intensity of a classical acoustic pumping. The analytic results obtained in previous sections can be readily generalized when the variation of a parameter is adiabatic, i.e., slower than the optical frequencies ω or W. Since the space is limited, the time-dependent problem will be considered elsewhere. Here, we consider just one example, namely turning on/off of a classical acoustic pumping q=Q(r)eiΩt+Q*(r)eiΩt.

For maximal control, it is beneficial to place an atom at the point where E(r=ra)0, whereas (Q)Er=ra is maximized. The equations of motion for quantum state amplitudes were derived in Section 5, see Eqs. (53)–(55).

Consider an exact parametric resonance ω+Ω=W for simplicity, when

ωn=ωn1+W.

The solution to Eqs. (53)–(55) when the acoustic pumping is turned off is given as follows:

Ψ=C00(0)eiω0t|0|0+n=1(Cn0(0)eiωnt|n|0+C(n1)1(0)ei(ωn1+W)t|n1|1)

The solution when the acoustic pumping is turned on is given as follows:

(72)Ψ=C00(0)eiω0t|0|0n=1[(Anei|ΩR(n)|t+Bnei|ΩR(n)|t)eiωnt|n|0+(Anei|ΩR(n)|t+Bnei|ΩR(n)|t)eiθi(ωn1+W)t|n1|1]

Assume that the initial quantum state before the pumping was turned on was not entangled, for example, an atom was in an excited state and there were no photons:

Ψ=ei(ω0+W)t|0|1.

If the acoustic pumping is turned on at t=0, the quantum state becomes entangled:

(73)Ψ=ieiω1tiθsin(|ΩR(1)|t)|1|0+ei(ω0+W)tcos(|ΩR(1)|t)|0|1,

Then the acoustic pumping can be turned off. Depending on the turnoff moment of time, one can obtain various entangled photon–atom states, e.g., Bell states, etc. The above reasoning is valid when the turn-on/off rate is slower than the optical frequencies and the detuning from the two-wave resonance ω=W.

8 Conclusions

In conclusion, we showed how the entanglement in a system of a fermionic quantum emitter coupled to a quantized EM field in a nanocavity and quantized phonon or mechanical vibrational modes emerges in the vicinity of a parametric resonance in the system. We developed analytic theory describing the formation and evolution of entangled quantum states, which can be applied to a broad range of cavity quantum optomechanics problems and emerging nanocavity strong coupling experiments. The model includes decoherence effects due to coupling of the fermion, photon, and phonon subsystems to their dissipative reservoirs within the stochastic evolution approach, which is derived from the Heisenberg–Langevin formalism. We showed that our approach provided the results for physical observables equivalent to those obtained from the density matrix equations with the relaxation operator in Lindblad form. We derived analytic expressions for the time evolution of the quantum state and observables and the emission spectra. The limit of a classical acoustic pumping, the control of entangled states, and the interplay between parametric and standard two-wave resonances were discussed.


Corresponding author: Alexey Belyanin, Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843, USA, E-mail:

Award Identifier / Grant number: 20-02-00100

Funding source: Texas A and M University

Funding source: Air Force Office for Scientific Research

Award Identifier / Grant number: FA9550-17-1-0341

Funding source: Federal Research Center Institute of Applied Physics of the Russian Academy of Sciences

Award Identifier / Grant number: 0035-2019-004

Award Identifier / Grant number: 1936276

  1. Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.

  2. Research funding: This work has been supported in part by the Air Force Office for Scientific Research through Grant No. FA9550-17-1-0341, National Science Foundation Award No. 1936276, and Texas A&M University through STRP, X-grant and T3-grant programs. M.T. acknowledges the support from the Russian Foundation for Basic Research Grant No. 20-02-00100, and M.E. acknowledges the support from the Federal Research Center Institute of Applied Physics of the Russian Academy of Sciences (Project No. 0035-2019-004).

  3. Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

Appendix A The stochastic equation of evolution for the state vector

The description of open quantum systems within the stochastic equation of evolution for the state vector is usually formulated for a Monte Carlo–type numerical scheme, e.g., the method of quantum jumps [23], [25]. We developed an approach suitable for analytic derivations. Our stochastic equation of evolution is basically the Schrödinger equation modified by adding a linear relaxation operator and the noise source term with appropriate correlation properties. The latter are related to the parameters of the relaxation operator in such a way that the expressions for the statistically averaged quantities satisfy certain physically meaningful conditions.

The protocol of introducing the relaxation operator with a corresponding noise source term to the quantum dynamics is well known in the Heisenberg picture, where it is called the Heisenberg–Langevin method [23], [34], [42]. We develop a conceptually similar approach for the Schrödinger equation. Here, we derive the general form of the stochastic equation of evolution from the Heisenberg–Langevin equations and track how certain physically reasonable constraints on the observables determine the correlation properties of the noise sources.

1. From Heisenberg–Langevin equations to the stochastic equation for the state vector

The Heisenberg–Langevin equation for the operator gˆ of a certain observable quantity takes the following form [23], [34], [42]:

(A1)ddtgˆ=i[Hˆ,gˆ]+Rˆ(gˆ)+Lˆg(t),

where Rˆ(gˆ) is the relaxation operator, Lˆg(t) is the Langevin noise source satisfying Lˆg(t)¯=0, where the bar means statistical averaging. For given commutation relations of the two operators, [gˆ1,gˆ2]=C, where C is a constant, correct Langevin sources should ensure the conservation of commutation relations at any moment of time, despite the presence of the relaxation operator in Eq. (A1), see [34], [43], [44].

The group of terms i[Hˆ,gˆ]+Rˆ(gˆ) can often be written as follows:

(A2)i[Hˆ,gˆ]+Rˆ(gˆ)=i(HˆeffgˆgˆHˆeff),

where Hˆeff is a non-Hermitian operator. For example, if the relaxation operator describes dissipation with relaxation constant γ so that gˆ¯eγt, then Hˆeff=Hˆiγ21ˆ, where 1ˆ is a unit operator. Note that in the master equation for the density matrix, the relaxation is often introduced in a conceptually similar way [25], [Hˆ,ρˆ]HˆeffρˆρˆHˆeff, which is however slightly different from the form used in Eq. (A2): [Hˆ,gˆ]HˆeffgˆgˆHˆeff. The difference is because the commutator of an unknown operator with Hamiltonian enters with opposite sign in the master equation as compared to the Heisenberg equation.

Now consider the transition from the Heisenberg–Langevin equation to the stochastic equation for the state vector. The key point is to assume that there exists the operator of evolution Uˆ(t), which is determined not only by the system parameters but also by the properties of the reservoir. This operator determines the evolution of the state vector:

(A3)|Ψ(t)=Uˆ(t)|Ψ0,Ψ(t)|=Ψ0|Uˆ(t),

where Ψ0=Ψ(0). Hereafter, we will denote the operators in the Schrödinger picture with index “s” to distinguish them from the Heisenberg operators. An observable can be calculated as follows:

g(t)=Ψ(t)|gˆS|Ψ(t)=Ψ0|gˆ(t)|Ψ0

which leads to

(A4)gˆ(t)=Uˆ(t)gˆSUˆ(t),

Since the substitution of Eqs. (A3) and (A4) into the standard Heisenberg equation leads to the standard Schrödinger equation, it makes sense to apply the same procedure to the Heisenberg–Langevin equation in order to obtain the “stochastic variant” of the Schrödinger equation. The solution of the latter should yield the expression for an observable,

g(t)=Ψ(t)|gˆS|Ψ(t)¯,

which is different from the standard expression by additional averaging over the noise statistics.

Note that an open system interacting with a reservoir is generally in a mixed state and should be described by the density matrix. We are describing the state of the system with a state vector which has a fluctuating component. For example, in a certain basis |α, the state vector will be Cα(t)=Cα¯+Cα˜, where the fluctuating component is denoted with a wavy bar. The elements of the density matrix of the corresponding mixed state are ραβ=CαCβ*¯=Cα¯Cβ*¯+Cα˜Cβ˜*¯.

The solution to the Heisenberg–Langevin equation can be expressed through the evolution operator Uˆ(t) using Eq. (A4). The noise source terms should be chosen to ensure the conservation of commutation relations at any moment of time, despite the presence of the relaxation operator. Since commutation relations between any two operators are conserved if and only if the evolution operator Uˆ(t) is unitary, a correct noise source in the Heisenberg–Langevin equation will automatically ensure the condition UˆUˆ=1ˆ.

We implement the above protocol. Substituting Eq. (A4) together with Hˆeff=UˆHˆeff,SUˆ and Hˆeff=UˆHˆeff,SUˆ into Eqs. (A1) and (A2) and using UˆUˆ=1ˆ, we arrive at the following:

(A5)(ddtUˆiUˆHˆeff,S)gˆSUˆ+UˆgˆS(ddtUˆ+iHˆeff,SUˆ)=Lˆg

Next, we introduce the operator Fˆ, defined as follows:

(A6)Lˆg=2UˆgˆSFˆ

For the operator Lˆg, Eq. (A6) gives Lˆg=2Fˆ(UˆgˆS)=2FˆgˆSUˆ. Since gˆ and gˆS are Hermitian operators, Lˆg has to be Hermitian too. (One can develop the Heisenberg–Langevin formalism for non-Hermitian operators too, for example, creation or annihilation operators, but the derivation becomes longer.) Then the operator Lˆg can be “split” between the two terms on the left-hand side of Eq. (A5) using the following relationship:

(A7)Lˆg=UˆgˆSFˆ+FˆgˆSUˆ

Substituting the latter into Eq. (A5), we obtain the following:

(ddtUˆiUˆHˆeff,SFˆ)gˆSUˆ+UˆgˆS(ddtUˆ+iHˆeff,SUˆFˆ)=0.

For simplicity, we will assume operator Hˆeff to be constant with time, i.e., we will not differentiate between Hˆeff and Hˆeff,S.

The last equation is satisfied for sure if

(A8)ddtUˆ=iHˆeffUˆ+Fˆ,ddtUˆ=iUˆHˆeff+Fˆ.

Multiplying Eq. (A8) by the initial state vector |Ψ0 from the right and from the left, we obtain the stochastic equation for the state vector and its Hermitian conjugate:

(A9)ddt|Ψ=iHˆeff|Ψi|R(t)
(A10)ddtΨ|=iΨ|Hˆeff+iR(t)|

where we introduced the notations iFˆ|Ψ0|R(t), iΨ0|FˆR(t)|. We will also need Eqs. (A9) and (A10) in a particular basis |α:

(A11)ddtCα=iν(Hˆeff)ανCνiRα,
(A12)ddtCα*=ivCν*(Hˆeff)να+iRα*,

where Rα=α|R, (Hˆeff)αβ=α|Hˆeff|β.

Applying the same procedure to the standard Heisenberg Eq. (1), we obtain that in Eqs. (A9) and (A10): HˆeffHˆeff=Hˆ and |Rt0, which corresponds to the standard Schrödinger equation and its Hermitian conjugate.

Note that intermediate relations (A8) for the evolution operator and in particular operator Fˆ should not depend on the choice of a particular physical observable g in the original Heisenberg–Langevin Eq. (A1). We assume that the Langevin operators in the original equation do not contradict this physically reasonable requirement.

In general, statistical properties of noise that ensure certain physically meaningful requirements impose certain constraints on the noise source |R which enters the right-hand side of the stochastic equation for the state vector. In particular, it is natural to require that the statistically averaged quantity |R¯=0. We will also require that the noise source |R has the correlation properties that preserve the norm of the state vector averaged over the reservoir statistics:

(A13)Ψ(t)|Ψ(t)¯=1.

2. Noise correlator

The solution to Eqs. (A9) and (A10) can be formally written as follows:

(A14)|Ψ=eiHˆefft|Ψ0i0teiHˆeff(τt)|R(τ)dτ,
(A15)Ψ|=Ψ0|eiHˆefft+i0tR(τ)|eiHˆeff(τt)dτ,

In the basis |α, Eqs. (A14) and (A15) can be transformed into the following equations:

(A16)Cα=α|eiHˆefft|Ψ0i0tα|eiHˆeff(τt)|R(τ)dτ,
(A17)Cα*=Ψ0|eiHˆefft|α+i0tR(τ)|eiHˆeff(τt)|αdτ.

In order to calculate the observables, we need to know the expressions for the averaged dyadic combinations of the amplitudes. We can find them using Eqs. (A11) and (A12):

(A18)ddtCαCβ*¯=iv(Hαν(h)CνCβ*¯CαCν*¯Hνβ(h))iv(Hαν(ah)CνCβ*¯+CαCν*¯Hνβ(ah))+(iCβ*Rα¯+iRβ*Cα¯),

where we separated the Hermitian and anti-Hermitian components of the effective Hamiltonian: α|Hˆeff|β=Hαβ(h)+Hαβ(ah). Substituting Eqs. (A16) and (A17) into the last term in Eq. (A18), we obtain the following:

iCβ*Rα¯+iCαRβ*¯=12t0R(t+ξ)|eiHˆeffξ|βα|R(t)¯dξ
+12t0R(t)|βα|eiHˆeffξ|R(t+ξ)¯dξ.

To proceed further with analytical results, we need to evaluate these integrals. The simplest situation is when the noise source terms are delta-correlated in time (Markovian). In this case, only the point ξ=0 contributes to the integrals. As a result, Eq. (A18) is transformed to the following equation:

(A19)ddtCαCβ*¯=iv(Hαν(h)CνCβ*¯CαCν*¯Hνβ(h))iv(Hαν(ah)CνCβ*¯+CαCν*¯Hνβ(ah))+Dαβ,

where the correlator Dαβ is defined as follows:

(A20)Rβ*(t+ξ)Rα(t)¯=Rβ*(t)Rα(t+ξ)¯=2δ(ξ)Dαβ

The time derivative of the norm of the state vector is given as follows:

(A21)ddtα|Cα|2¯=α[iν(Hαν(ah)CνCα*¯+CαCν*¯Hνα(ah))Dαα]

Clearly, the components Dαα of the noise correlator need to compensate the decrease in the norm due to the anti-Hermitian component of the effective Hamiltonian. Therefore, the expressions for Hαβ(ah) and Dαα have to be mutually consistent. This is the manifestation of the fluctuation–dissipation theorem [41].

Note that the noise correlator could depend on the averaged combinations (e.g., dyadics) of the components of the state vector. This is because the noise source term |R introduced above depends on the initial state |Ψ0 and the evolution operator Uˆ, and these quantities form the state vector components at any given time. Of course, what we call a “state vector” is the solution of the stochastic equation of motion, which is very different from the solution of the conventional Schrödinger equation for a closed system. In particular, we postulated the existence of the evolution operator Uˆ determined not only by the parameters of the dynamical system but also by the properties of a dissipative reservoir, although we did not specify any particular expression for Uˆ.

As an example, consider a simple diagonal anti-Hermitian operator Hαν(ah):

(A22)Hαν(ah)=iγαδαν

and introduce the following models:

  1. Populations relax much slower than coherences (expected for condensed matter systems). In this case, we can choose Dαβ=0, Dαα=2γα|Cα|2¯; within this model, the population at each state will be preserved.

  2. The state α=αdown has a minimal energy, while the reservoir temperature T=0. In this case, it is expected that all populations approach zero in equilibrium, whereas the occupation number of the ground state approaches 1, similar to the Weisskopf–Wigner model. The adequate choice of correlators is Dαβ=0, Dααδααdown, γαdown=0. The expression for the remaining nonzero correlator,

(A23)Dαdownαdown=ααdown2γα|Cα|2¯,

ensures the conservation of the norm:

ddtααdown|Cα|2¯=ααdown2γα|Cα|2¯=ddt|Cαdown|2¯.

This is an example of the correlator’s dependence on the state vector that we discussed before.

3. Comparison with the Lindblad method

One can choose the anti-Hermitian Hamiltonian Hαβ(ah) and correlators Dαβ in the stochastic equation of motion in such a way that Eq. (A19) for the dyadics CnCm*¯ corresponds exactly to the equations for the density matrix elements in the Lindblad approach. Indeed, the Lindblad form of the master equation has the following form [23], [25]:

(A24)ddtρˆ=i[Hˆ,ρˆ]+Lˆ(ρˆ)

where Lˆ(ρˆ) is the Lindbladian:

(A25)Lˆ(ρˆ)=12kγk(lˆklˆkρˆ+ρˆlˆklˆk2lˆkρˆlˆk),

Operators lˆk in Eq. (A25) and their number are determined by the model which describes the coupling of the dynamical system to the reservoir. The form of the relaxation operator given by Eq. (A25) preserves automatically the conservation of the trace of the density matrix, whereas the specific choice of relaxation constants ensures that the system approaches a proper steady state given by thermal equilibrium or supported by an incoherent pumping.

Eq. (A24) is convenient to represent in a slightly different form:

(A26)ddtρˆ=i(HˆeffρˆρˆHˆeff)+δLˆ(ρˆ)

where

(A27)Hˆeff=Hˆikγklˆklˆk,δLˆ(ρˆ)=kγklˆkρˆlˆk.

Writing the anti-Hermitian component of the Hamiltonian in Eqs. (A11) and (A12) as

(A28)Hαβ(ah)=iα|kγklˆklˆk|β,

and defining the corresponding correlator of the noise source as

(A29)Rβ*(t+ξ)Rα(t)¯=2δ(ξ)Dαβ,Dαβ=α|δLˆ(ρˆ)|βρmn=CnCm*¯,

We obtain the solution in which averaged over noise statistics dyadics CnCm*¯ correspond exactly to the elements of the density matrix within the Lindblad method.

Instead of deriving the stochastic equation of evolution of the state vector from the Heisenberg–Langevin equations, we could postulate it from the very beginning. After that, we could justify the choice of the effective Hamiltonian and noise correlators by ensuring that they lead to the same observables as the solution of the density matrix equations with the relaxation operator in the Lindblad form [25], [46], [47]. However, the demonstration of direct connection between the stochastic equation of evolution of the state vector and the Heisenberg–Langevin equation provides an important physical insight.

4. Relaxation rates for coupled subsystems interacting with a reservoir

Whenever we have several coupled subsystems (such as electrons, photon modes, and phonons in this paper), each coupled to its reservoir, the determination of relaxation rates of the whole system becomes nontrivial. The problem can be solved if we assume that these “partial” reservoirs are statistically independent. In this case, it is possible to add up partial Lindbladians and obtain the total effective Hamiltonian.

Consider the Hamiltonian (11) of the system formed by a two-level electron system coupled to an EM mode field and dressed by a phonon field:

(A30)Hˆ=Hˆem+Hˆa+Hˆp+Vˆ.

where, Hˆem=ω2(cˆcˆ+cˆcˆ) is the Hamiltonian for a single EM mode field, Hˆa=W1σˆσˆ+W0σˆσˆ is the Hamiltonian for a two-level “atom” with energy levels W0,1, Hˆp=Ω2(bˆbˆ+bˆbˆ) is the Hamiltonian for a phonon mode, Vˆ=Vˆ1+Vˆ2 is the interaction Hamiltonian, where Vˆ1,2 describe the atom–photon and atom–photon–phonon coupling, respectively:

Vˆ1=(χσˆcˆ+χ*σˆcˆ+χσˆcˆ+χ*σˆcˆ),Vˆ2=(η1σˆcˆbˆ+η1*σˆcˆbˆ+η2σˆcˆbˆ+η2*σˆcˆbˆ+η1σˆcˆbˆ+η1*σˆcˆbˆ+η2σˆcˆbˆ+η2*σˆcˆbˆ),

where χ, η1, η2 are coupling constants defined before.

Summing up the known (see e.g., [23], [25]) partial Lindbladians of two bosonic (infinite amount of energy levels) and one fermionic (two-level) subsystems, we obtain the following:

(A31)L(ρˆ)=γ2N1Ta(σˆσˆρˆ+ρˆσˆσˆ2σˆρˆσˆ)γ2N0Ta(σˆσˆρˆ+ρˆσˆσˆ2σˆρˆσˆ)μω2n¯ωTem(cˆcˆρˆ+ρˆcˆcˆ2cˆρˆcˆ)μω2(n¯ωTem+1)(cˆcˆρˆ+ρˆcˆcˆ2cˆρˆcˆ)μΩ2n¯ΩTp(bˆbˆρˆ+ρˆbˆbˆ2bˆρˆbˆ)μΩ2(n¯ΩTp+1)(bˆbˆρˆ+ρˆbˆbˆ2bˆρˆbˆ)

where γ, μω, and μΩ are partial relaxation rates of the systems,

N0,1Ta=(1+eW1W0Ta)1eW0,1W0Ta,n¯ωTem=(eωTem1)1,n¯ΩTp=(eΩTp1)1,
Ta,em,p are the temperatures of partial reservoirs. For the Lindblad master equation in the form Eq. (A26), we get the following:
(A32)Hˆeff=HˆiΓˆ,

where

(A33)Γˆ=2{γ(N1Taσˆσˆ+N0Taσˆσˆ)+μω[n¯ωTemcˆcˆ+(n¯ωTem+1)cˆcˆ]+μΩ[n¯ΩTpbˆbˆ+(n¯ΩTp+1)bˆbˆ]}.

Using the effective Hamiltonian given by Eqs. (A32) and (A33), we arrive at the stochastic equation for the state vector in the following form:

(A34)ddtCαn0=iW0+ω(n+12)+Ω(α+12)Cαn0iα|n|0|Vˆ|Ψγαn0Cαn0iRαn0,
(A35)ddtCαn1=iW1+ω(n+12)+Ω(α+12)Cαn1iα|n|1|Vˆ|Ψγαn1Cαn1iRαn1,

where

(A36)γαn0=γ2N1Ta+μω2[n¯ωTem(n+1)+(n¯ωTem+1)n]+μΩ2[n¯ΩTp(α+1)+(n¯ΩTp+1)α],
(A37)γαn1=γ2N0Ta+μω2[n¯ωTem(n+1)+(n¯ωTem+1)n]+μΩ2[n¯ΩTp(α+1)+(n¯ΩTp+1)α],
Eqs. (A36) and (A37) determine the rules of combining the “partial” relaxation rates for several coupled subsystems.

5. Including purely elastic dephasing processes

So far, we used the Lindbladian which includes only the dissipation and does not include purely elastic dephasing processes. In order to take them into account, we need to modify the Lindbladian to include the explicit dependence on the operators of populations. We will follow the prescription which can be found in the study by Fain and Khanin [48]. Using Eq. (A25) for the “partial” Lindbladian Lˆ2l of a two-level atom, where k=1,2,3, lˆ1=σˆ, lˆ2=σˆ, and lˆ3=σˆσˆσˆσˆ, we obtain the following:

(A38)L2l(ρˆ)=γ2N1T(σˆσˆρˆ+ρˆσˆσˆ2σˆρˆσˆ)γ2N0T(σˆσˆρˆ+ρˆσˆσˆ2σˆρˆσˆ)γel2ρˆ+δLˆel(ρˆ),

where

(A39)δLˆel(ρˆ)=γel2(σˆσˆσˆσˆ)ρˆ(σˆσˆσˆσˆ).

For the evolution of a two-level system, using the Lindbladian (A39) leads to standard density matrix equations with inverse relaxation times for the coherence, 1T2=γ2+γel, and populations, 1T1=γ.

Furthermore, using the scheme developed in Section 3 of the Appendix, we obtain that adding elastic processes to the stochastic equation of evolution for the state vector leads to the following modifications for the anti-Hermitian part of the effective Hamiltonian,

Hαβ(ah)Hαβ(ah)iδαβ4γel,

and the correlators of noise sources,

DαβDαβ+γel(δαβ|Cα|2¯12CαCβ*).

This is a general prescription. Since, in this work, we are only interested in the RWA dynamics of states |α|n|0 and |α1|n1|1, the same expressions for the observables can be obtained with a much simpler modification of the formalism. One can show that it is sufficient to modify the relaxation constants γ(α1)(n1)1 according to

γ(α1)(n1)1γ(α1)(n1)1+γel

and correlators D(α1)(n1)1;(α1)(n1)1 as

D(α1)(n1)1;(α1)(n1)1D(α1)(n1)1;(α1)(n1)1+2γel|C(α1)(n1)1|2¯.

References

[1] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, “Cavity optomechanics,” Rev. Mod. Phys., vol. 86, p. 1391, 2014, https://doi.org/10.1103/revmodphys.86.1391.Search in Google Scholar

[2] P. Meystre, “A short walk through quantum optomechanics,” Ann. Phys., vol. 525, p. 215, 2013, https://doi.org/10.1002/andp.201200226.Search in Google Scholar

[3] J.-M. Pirkkalainen, S.U. Cho, F. Massel, et al., “Cavity optomechanics mediated by a quantum two-level system,” Nat. Commun., vol. 6, p. 6981, 2015, https://doi.org/10.1038/ncomms7981.Search in Google Scholar PubMed PubMed Central

[4] Y. Chu, P. Kharel, W. H. Renninger, et al., “Quantum acoustics with superconducting qubits,” Science, vol. 358, p. 199, 2017, https://doi.org/10.1126/science.aao1511.Search in Google Scholar PubMed

[5] S. Hong, R. Riedinger, I. Marinkovic, et al., “Hanbury Brown and Twiss interferometry of single phonons from an optomechanical resonator,” Science, vol. 358, p. 203, 2017, https://doi.org/10.1126/science.aan7939.Search in Google Scholar PubMed

[6] P. Arrangoiz-Arriola, E. A. Wollack, Z. Wang, et al., “Resolving the energy levels of a nanomechanical oscillator,” Nature, vol. 571, p. 537, 2019, https://doi.org/10.1038/s41586-019-1386-x.Search in Google Scholar PubMed

[7] F. Benz, M. K. Schmidt, A. Dreismann, et al., “Single-molecule optomechanics in “picocavities”,” Science, vol. 354, p. 726, 2016, https://doi.org/10.1126/science.aah5243.Search in Google Scholar PubMed

[8] K.-D. Park, E. A. Muller, V. Kravtsov, et al., “Variable-temperature tip-enhanced raman spectroscopy of single-molecule fluctuations and dynamics,” Nano Lett., vol. 16, p. 479, 2016, https://doi.org/10.1021/acs.nanolett.5b04135.Search in Google Scholar PubMed

[9] M. S. Tame, K. R. McEnery, S. K. Ozdemir, J. Lee, S. A. Maier, and M. S. Kim, “Quantum plasmonics,” Nat. Phys., vol. 9, p. 329, 2013, https://doi.org/10.1038/nphys2615.Search in Google Scholar

[10] F. C. B. Maia, B. T. OCallahan, A. R. Cadore, et al., “Anisotropic flow control and gate modulation of hybrid phonon-polaritons,” Nano Lett., vol. 19, p. 708, 2019, https://doi.org/10.1021/acs.nanolett.8b03732.Search in Google Scholar PubMed

[11] R. Chikkaraddy, B. de Nijs, F. Benz, et al., “Single-molecule strong coupling at room temperature in plasmonic nanocavities,” Nature, vol. 535, p. 127, 2016, https://doi.org/10.1038/nature17974.Search in Google Scholar PubMed PubMed Central

[12] A. Sipahigil, R. E. Evans, D. D. Sukachev, et al., “An integrated diamond nanophotonics platform for quantum-optical networks,” Science, vol. 354, p. 847, 2016, https://doi.org/10.1126/science.aah6875.Search in Google Scholar PubMed

[13] T. Yoshie, A. Scherer, J. Hendrickson, et al., “Vacuum Rabi splitting with a single quantum dot in a photonic crystal nanocavity,” Nature, vol. 432, p. 200, 2004, https://doi.org/10.1038/nature03119.Search in Google Scholar PubMed

[14] J. P. Reithmaier, G. Sek, A. Loffler, et al., “Strong coupling in a single quantum dot-semiconductor microcavity system,” Nature, vol. 432, p. 197, 2004, https://doi.org/10.1038/nature02969.Search in Google Scholar PubMed

[15] H. Leng, B. Szychowski, M.-C. Daniel, and M. Pelton, “Strong coupling and induced transparency at room temperature with single quantum dots and gap plasmons,” Nat. Commun., vol. 9, p. 4012, 2018, https://doi.org/10.1038/s41467-018-06450-4.Search in Google Scholar PubMed PubMed Central

[16] O. Bitton, S. N. Gupta, and G. Haran, “Quantum dot plasmonics: from weak to strong coupling,” Nanophotonics, vol. 8, p. 559, 2019, https://doi.org/10.1515/nanoph-2018-0218.Search in Google Scholar

[17] K.-D. Park, M. A. May, H. Leng, et al., “Tip-enhanced strong coupling spectroscopy, imaging, and control of a single quantum emitter,” Sci. Adv., vol. 5, p. eaav5931, 2019, https://doi.org/10.1126/sciadv.aav5931.Search in Google Scholar PubMed PubMed Central

[18] K. J. Satzinger, Y. P. Zhong, H. Chang, et al., “Quantum control of surface acoustic-wave phonons,” Nature, vol. 563, p. 661665, 2018, https://doi.org/10.1038/s41586-018-0719-5.Search in Google Scholar PubMed

[19] A. Bienfait, Y. P. Zhong, H.-S. Chang, et al., “Quantum erasure using entangled surface acoustic phonons,” Phys. Rev. X, vol. 10, p. 021055, 2020, https://doi.org/10.1103/physrevx.10.021055.Search in Google Scholar

[20] U. Delic, M. Reisenbauer, K. Dare, et al., “Cooling of a levitated nanoparticle to the motional quantum ground state,” Science, vol. 367, p. 892, 2020, https://doi.org/10.1126/science.aba3993.Search in Google Scholar PubMed

[21] N. Bloembergen, Nonlinear Optics, Singapore, World Scientific, 1996.10.1142/3046Search in Google Scholar

[22] Y. A. Bogoliubov and N. N. Mitropolsky, Asymptotic Methods in the Theory of Nonlinear Oscillations, London, Gordon and Breach, 1961.Search in Google Scholar

[23] M. O. Scully and M. S. Zubairy, Quantum Optics, Cambridge, Cambridge University Press, 1997.10.1017/CBO9780511813993Search in Google Scholar

[24] P. Forn-Diaz, L. Lamata, E. Rico, J. Kono, and E. Solano, “Ultrastrong coupling regimes of light-matter interaction,” Rev. Mod. Phys., vol. 91, p. 025005, 2019, https://doi.org/10.1103/revmodphys.91.025005.Search in Google Scholar

[25] M. B. Plenio and P. L. Knight, “The quantum-jump approach to dissipative dynamics in quantum optics,” Rev. Mod. Phys., vol. 70, p. 101, 1998, https://doi.org/10.1103/revmodphys.70.101.Search in Google Scholar

[26] C. Genes, D. Vitali, and P. Tombesi, “Emergence of atom-light-mirror entanglement inside an optical cavity,” Phys. Rev. A, vol. 77, no. R, p. 050307, 2008, https://doi.org/10.1103/physreva.77.050307.Search in Google Scholar

[27] G. Chiara, M. Paternostro, and G. M. Palma, “Entanglement detection in hybrid optomechanical systems,” Phys. Rev. A, vol. 83, p. 052324, 2011, https://doi.org/10.1103/physreva.83.052324.Search in Google Scholar

[28] Y. Xiang, F. X. Sun, M. Wang, Q. H. Gong, and Q. Y. He, “Detection of genuine tripartite entanglement and steering in hybrid optomechanics,” Opt. Expr., vol. 23, p. 30104, 2015, https://doi.org/10.1364/oe.23.030104.Search in Google Scholar

[29] Y. Wang, S. Chesi, and A. A. Clerk, “Bipartite and tripartite output entanglement in three-mode optomechanical systems,” Phys. Rev. A, vol. 91, p. 013807, 2015, https://doi.org/10.1103/physreva.91.013807.Search in Google Scholar

[30] X. Yang, Y. Ling, X. Shao, and M. Xiao, “Generation of robust tripartite entanglement with a single-cavity optomechanical system,” Phys. Rev. A, vol. 85, p. 052303, 2017, https://doi.org/10.1103/physreva.95.052303.Search in Google Scholar

[31] Q. Liao, Y. Ye, P. Lin, N. Zhou, and W. Nie, “Tripartite entanglement in an atom-cavity-optomechanical system,” Int. J. Theor. Phys., vol. 57, p. 1319, 2018, https://doi.org/10.1007/s10773-017-3661-7.Search in Google Scholar

[32] C. Jiang, S. Tserkis, K. Collins, S. Onoe, Y. Li, and L. Tian, “Switchable bipartite and genuine tripartite entanglement via an optoelectromechanical interface,” Phys. Rev. A, vol. 101, p. 042320, 2020, https://doi.org/10.1103/physreva.101.042320.Search in Google Scholar

[33] M. Tokman, Y. Wang, I. Oladyshkin, A. R. Kutayiah, and A. Belyanin, “Laser-driven parametric instability and generation of entangled photon-plasmon states in graphene,” Phys. Rev. B, vol. 93, p. 235422, 2016, https://doi.org/10.1103/physrevb.93.235422.Search in Google Scholar

[34] M. Tokman, X. Yao, and A. Belyanin, “Generation of entangled photons in graphene in a strong magnetic field,” Phys. Rev. Lett., vol. 110, p. 077404, 2013, https://doi.org/10.1103/physrevlett.110.077404.Search in Google Scholar PubMed

[35] M. D. Tokman, M. A. Erukhimova, and V. V. Vdovin, “The features of a quantum description of radiation in an optically dense medium,” Ann. Phys., vol. 360, p. 571, 2015, https://doi.org/10.1016/j.aop.2015.05.030.Search in Google Scholar

[36] M. Tokman, Z. Long, S. AlMutairi, Y. Wang, M. Belkin, and A. Belyanin, “Enhancement of the spontaneous emission in subwavelength quasi-two-dimensional waveguides and resonators,” Phys. Rev. A, vol. 97, p. 043801, 2018, https://doi.org/10.1103/physreva.97.043801.Search in Google Scholar

[37] W. Dur, G. Vidal, and J. I. Cirac, “Three qubits can be entangled in two inequivalent ways,” Phys. Rev. A, vol. 62, p. 062314, 2000, https://doi.org/10.1103/physreva.62.062314.Search in Google Scholar

[38] M. M. Cunha, A. Fonseca, and E. O. Silva, arXiv:1909.00862v2.Search in Google Scholar

[39] L. K. Shalm, D. R. Hamel, Z. Yan, C. Simon, K. J. Resch, and T. Jennewein, “Three-photon energy-time entanglement,” Nat. Phys., vol. 9, p. 19, 2012, https://doi.org/10.1038/nphys2492.Search in Google Scholar

[40] A. Agusti, C. W. Sandbo Chang, F. Quijandria, G. Johansson, C. M. Wilson, and C. Sabin, “Tripartite genuine non-Gaussian entanglement in three-mode spontaneous parametric down-conversion,” Phys. Rev. Lett., vol. 125, p. 020502, 2020, https://doi.org/10.1103/physrevlett.125.020502.Search in Google Scholar

[41] L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 1, Oxford, Pergamon, 1965.Search in Google Scholar

[42] C. Gardiner and P. Zoller, Quantum Noise, Berlin, Heidelberg, Springer-Verlag, 2004.Search in Google Scholar

[43] M. Erukhimova and M. Tokman, “Squeezing of thermal fluctuations in four-wave mixing in a scheme,” Phys. Rev A, vol. 95, p. 013807, 2017, https://doi.org/10.1103/physreva.95.013807.Search in Google Scholar

[44] M. Tokman, Z. Long, S. Almutairi, et al., “Purcell enhancement of the parametric down-conversion in two-dimensional nonlinear materials,” APL Photonics, vol. 4, p. 034403, 2019, https://doi.org/10.1063/1.5044539.Search in Google Scholar

[45] K. H. Madsen and P. Lodahl, “Quantitative analysis of quantum dot dynamics and emission spectra in cavity quantum electrodynamics,” New J. Phys., vol. 15, p. 025013, 2013, https://doi.org/10.1088/1367-2630/15/2/025013.Search in Google Scholar

[46] S. Haroche and J.-M. Raymond, Exploring the Quantum. Atoms, Cavities, and Photons, Oxford, UK, Oxford University Press, 2006.10.1093/acprof:oso/9780198509141.001.0001Search in Google Scholar

[47] K. Blum, Density Matrix Theory and Applications, Heidelberg, Springer, 2012.10.1007/978-3-642-20561-3Search in Google Scholar

[48] V. M. Fain and Y. I. Khanin, Quantum Electronics. Basic Theory, Cambridge, MA, MIT, 1969.Search in Google Scholar

Received: 2020-06-29
Accepted: 2020-08-22
Published Online: 2020-09-15

© 2020 Mikhail Tokman et al., published by De Gruyter, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.

Downloaded on 28.3.2024 from https://www.degruyter.com/document/doi/10.1515/nanoph-2020-0353/html
Scroll to top button