Generation and dynamics of entangled fermion-photon-phonon states in nanocavities
Mikhail Tokman, Maria Erukhimova, Yongrui Wang, Qianfan Chen, Alexey Belyanin
GGeneration and dynamics of entangled fermion-photon-phononstates in nanocavities
Mikhail Tokman, Maria Erukhimova, YongruiWang, Qianfan Chen, and Alexey Belyanin Institute of Applied Physics, Russian Academyof Sciences, Nizhny Novgorod, 603950, Russia Department of Physics and Astronomy,Texas A&M University, College Station, TX, 77843 USA (Dated: July 7, 2020)
Abstract
We develop the analytic theory describing the formation and evolution of entangled quantumstates for a fermionic quantum emitter coupled to a quantized electromagnetic field in a nanocav-ity and quantized phonon or mechanical vibrational modes. The theory is applicable to a broadrange of cavity quantum optomechanics problems and emerging research on plasmonic nanocavitiescoupled to single molecules and other quantum emitters. The optimal conditions for a tri-stateentanglement are realized near the parametric resonances in a coupled system. The model in-cludes decoherence effects due to coupling of the fermion, photon, and phonon subsystems totheir dissipative reservoirs within the stochastic evolution approach, which is derived from theHeisenberg-Langevin formalism. Our theory provides analytic expressions for the time evolutionof the quantum state and observables, and the emission spectra. The limit of a classical acousticpumping and the interplay between parametric and standard one-photon resonances are analyzed. a r X i v : . [ qu a n t - ph ] J u l . INTRODUCTION There is a lot of recent interest in the quantum dynamics of fermion systems coupled toan electromagnetic (EM) mode in a cavity and quantum or classical mechanical/acousticoscillations or phonon vibrations. This problem is related to the burgeoning fields of cavityoptomechanics [1–3] and quantum acoustics [4–6]. Another example where this situation canbe realized is a molecule placed in a plasmonic nanocavity [7, 8]. In this case the fermionsystem may comprise two or more electron states forming an optical transition, whereas thephonon field is simply a vibrational mode of a molecule. One can also imagine a situationwhere a quantum emitter such as a quantum dot or an optically active defect in a solid matrixis coupled to the quantized phonon modes of a crystal lattice, which would be an extensionof 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 quan-tum information applications. Indeed, the presence of a classical or quantized acoustic modeprovides an extra handle to control the quantum state of a coupled fermion-boson quantumsystem. In the extreme quantum limit in which the fermionic degree of freedom and allbosonic degrees of freedom (both photons and phonons) are quantized, a strong enough cou-pling between them leads to an entangled fermion-photon-phonon state, which is a complexenough system to implement basic gates for quantum computation or other applications.Such a system has not been realized experimentally. However, many ingredients have beenalready demonstrated, such as strong coupling between a nanocavity mode and a singlemolecule [11], numerous examples of strong coupling between nanocavity modes and a sin-gle fermionic quantum emitter such as a color center [12] or a quantum dot (QD) (see e.g.Refs. [13, 14] for semiconductor cavity-QD systems and Refs. [15–17] for plasmonic cavities),strong coupling and entanglement of acoustic phonons [18, 19], resolving the energy levelsof a nanomechanical oscillator [6], or cooling a macroscopic system into its motional groundstate [20].Interaction of three or more modes of oscillations, whether they are classical or quan-tized, is strongly enhanced close to the parametric resonance, which is therefore the mostinteresting region to study. Fortunately for theorists, the analysis near the parametric res-onance is greatly simplified, because some form of a slowly varying amplitude method for2lassical 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 lowerthan the characteristic energies in the system, such as the optical transition or vibrationalenergy. The emerging studies of the so-called ultra-strong coupling regime [24] have to gobeyond the RWA. Nevertheless, for the vast majority of experiments, including nonpertur-bative strong coupling dynamics and entanglement, the RWA is adequate and provides somecrucial simplifications that allow one to obtain analytic solutions.In particular, within Schr¨odinger’s description, the equations of motion for the compo-nents of an infinitely dimensional state vector | Ψ (cid:105) that describes a coupled fermion-bosonsystem can be split into the blocks of low dimensions if the RWA is applied. This is true evenif the dynamics of the fermion subsystem is nonperturbative, e.g. the effects of saturationare important. Note that there is no such simplification in the Heisenberg representation,i.e. when solving the equations of motion ddt ˆ g = i (cid:126) (cid:104) ˆ H, ˆ g (cid:105) , (1)where (cid:98) g is the Heisenberg operator of a certain physical observable g and ˆ H is the Hamilto-nian of the system. Operator-valued Eqs. (1) are generally impossible to split into smallerblocks, even within the RWA. This happens because some matrix elements g AB ( t ) of theHeisenberg operator are determined by states | A (cid:105) , | B (cid:105) which belong to different blocks thatevolve independently in the Schr¨odinger picture. The simplification could only be possiblefor specially selected initial conditions in which the Heisenberg operator is determined ona “truncated” basis belonging to only one of the independent blocks. The Schr¨odinger’sapproach also leads to fewer equations for the state vector components than the approachbased on the von Neumann master equation for the elements of the density matrix.Obviously, the Schr¨odinger equation in its standard form cannot be applied to describeopen systems coupled to a dissipative reservoir. In this case the stochastic versions of theequation of evolution for the state vector have been developed, e.g. the method of quantumjumps [23, 25]. This method is optimal for numerical analysis in the Monte-Carlo typeschemes. Here we formulate the stochastic equation for the state vector derived from theHeisenberg-Langevin approach which is more conducive to the analytic treatment. Its keyelement is an assumption that there exists the operator of evolution ˆ U , which is determinedunambiguously not only by the parameters of the dynamical system but also by the statistical3roperties of a dissipative reservoir.The paper is structured as follows. Section II formulates the model and the Hamiltonianfor coupled quantized fermion, photon, and phonon fields in a nanocavity. Section IIIderives the solution for the quantum states of a closed system in the vicinity of a parametricresonance and analyze its properties. In Section IV we provide the stochastic equationdescribing the evolution of quantum states of an open system in contact with a dissipativereservoir and describe the observables. In Section V we consider the case of a classicalacoustic pumping. Section VI describes the interplay of parametric and standard one-photonresonances and provides the conditions under which these resonances can be separated.Section VII gives an example of manipulating entangled electron-photon states by an acousticpumping. Appendix contains the derivation of the stochastic equation of evolution from theHeisenberg-Langevin approach and compares with Lindblad density-matrix formalism. II. A COUPLED QUANTIZED ELECTRON-PHOTON-PHONON SYSTEM: THEMODEL
Consider a quantized electron system coupled to the quantum EM field of a nanocavityand classical or quantized vibrational (phonon) modes, see Fig. 1 which sketches two out ofmany possible scenarios.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 phononsubsystems due to couplings to their respective dissipative reservoirs are also indicated.Figure 1b implies that it is a quantum dot which experiences vibrations, but our treatmentbelow works for any mechanism of relative displacement between the electron system andthe field of an EM cavity mode, including the situations where it is the wall of a nanocavitywhich 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-amplitudeacoustic oscillations, classical vs. quantum phonon mode, etc.4 Au tipAu substrate µ ω μ Ω Ωγ (a) µ ω μ Ω Ωγ W (b) FIG. 1: (a) A sketch of a molecule in a nanocavity created by a metallic nanotip and asubstrate; (b) A sketch of a quantum dot coupled to optical and mechanical vibrationalmodes in a nanocavity.
A. The fermion subsystem
Consider the simplest version of the fermion subsystem: two electron states | (cid:105) and | (cid:105) with energies 0 and W , respectively. We will call it an “atom” for brevity, although it canbe electron states of a molecule, a quantum dot, or any other electron system. Introducecreation and annihilation operators of the excited state | (cid:105) , ˆ σ = | (cid:105) (cid:104) | , ˆ σ † = | (cid:105) (cid:104) | , whichsatisfy standard commutation relations for fermions:ˆ σ † | (cid:105) = | (cid:105) , ˆ σ | (cid:105) = | (cid:105) , ˆ σ ˆ σ = ˆ σ † ˆ σ † = 0; (cid:2) ˆ σ, ˆ σ † (cid:3) + = ˆ σ ˆ σ † + ˆ σ † ˆ σ = 1 . The Hamiltonian of an atom is ˆ H a = W ˆ σ † ˆ σ. (2)5e will also need the dipole moment operator, ˆd = d (cid:0) ˆ σ † + ˆ σ (cid:1) , (3)where d = (cid:104) | ˆd | (cid:105) is a real vector. For a finite motion we can always choose the coordinaterepresentation of stationary states in terms of real functions. B. Quantized EM modes of a cavity
We use a standard representation for the electric field operator in a cavity: ˆE = (cid:88) i (cid:104) E i ( r ) ˆ c i + E ∗ i ( r ) ˆ c † i (cid:105) , (4)where ˆ c † i , ˆ c i are creation and annihilation operators for photons at frequency ω i ; the functions E i ( r ) describe the spatial structure of the EM modes in a cavity. The relation between themodal frequency ω i and the function E i ( r ) can be found by solving the boundary-valueproblem of the classical electrodynamics [23]. The normalization conditions [26] (cid:90) V ∂ [ ω i ε ( ω i , r )] ω i ∂ω i E ∗ i ( r ) E i ( r ) d r = 4 π (cid:126) ω i (5)ensure correct bosonic commutators (cid:104) ˆ c i , ˆ c † i (cid:105) = δ ij and the field Hamiltonian in the formˆ H em = (cid:126) (cid:88) i ω i (cid:18) ˆ c † i ˆ c i + 12 (cid:19) . (6)Here V is a quantization volume and ε ( ω, r ) is the dielectric function of a dispersive mediumthat fills the cavity. C. The quantized phonon field
We assume that our two-level atom is dressed by a phonon field which can be describedby the displacement operator: ˆq = (cid:88) i ˆq i ; ˆq i = Q i ( r ) ˆ b i + Q ∗ i ( r ) ˆ b † i (7)Here ˆ b i and ˆ b † i are annihilation and creation operators of phonons and the functions Q i ( r )determine the spatial structure of oscillations at frequencies Ω i . Expression (7) can be used6hen the amplitude of oscillations is small enough. One can always choose the normalizationof functions Q i ( r ) corresponding to standard commutation relations for bosons, (cid:104) ˆ b i , ˆ b † j (cid:105) = δ ij and a standard form for the Hamiltonian of mechanical oscillations:ˆ H p = (cid:126) (cid:88) i Ω i (cid:18) ˆ b † i ˆ b i + 12 (cid:19) . (8) D. An atom coupled to quantized EM and phonon fields
Now we can combine all ingredients into a coupled quantized system. Adding the inter-action Hamiltonian with a EM cavity mode in the electric dipole approximation, − ˆd · ˆE , theHamiltonian of an atom coupled to a single mode EM fieldˆ H = ˆ H em + ˆ H a − d (cid:0) ˆ σ † + ˆ σ (cid:1) · (cid:2) E ( r ) ˆ c + E ∗ ( r ) ˆ c † (cid:3) r = r a , (9)where r = r a denotes the position of an atom inside the cavity. The effect of “dressing” ofthe coupled atom-EM field system by mechanical oscillations in its most general form canbe included by adding the Hamiltonian of phonon modes ˆ H p and substituting r a = ⇒ r a + ˆq in Eq. (9). This will work for an arbitrary relative displacement of an atom with respect tothe EM cavity mode. Keeping only one phonon mode for simplicity, in which ˆq = Q ( r ) ˆ b + Q ∗ ( r ) ˆ b † , (10)and expanding in Taylor series in the vicinity of r = r a , we obtain the total Hamiltonian,ˆ H = ˆ H em + ˆ H a + ˆ H p − (cid:0) χ ˆ σ † ˆ c + χ ∗ ˆ σ ˆ c † + χ ˆ σ ˆ c + χ ∗ ˆ σ † ˆ c † (cid:1) − (cid:16) η ˆ σ † ˆ c ˆ b + η ∗ ˆ σ ˆ c † ˆ b † + η ˆ σ † ˆ c ˆ b † + η ∗ ˆ σ ˆ c † ˆ b + η ˆ σ ˆ c ˆ b + η ∗ ˆ σ † ˆ c † ˆ b † + η ˆ σ ˆ c ˆ b † + η ∗ ˆ σ † ˆ c † ˆ b (cid:17) (11)where χ = (cid:16) ˆd · ˆE (cid:17) r = r a , η = [ d ( Q · ∇ ) E ] r = r a , η = [ d ( Q ∗ ·∇ ) E ] r = r a . Note that we can always take the functions E ( r ) and Q ( r ) to be real at the position ofan atom , but we cannot keep the derivatives real at the same time if the modal structure ∝ e i k · r . However, for ideal cavity modes the latter is possible. As we will see below, thebest conditions for electron-photon-phonon entanglement are reached in the vicinity of the parametric resonance : W (cid:126) ≈ ω ± Ω. (12)7hen the upper sign is chosen in Eq. (12), the RWA applied to the Hamiltonian (11) yieldsˆ H = ˆ H em + ˆ H a + ˆ H p − (cid:16) η ˆ σ † ˆ c ˆ b + η ∗ ˆ σ ˆ c † ˆ b † (cid:17) (13)where η ≡ η . For the lower sign in Eq. (12), the RWA Hamiltonian isˆ H = ˆ H em + ˆ H a + ˆ H p − (cid:16) η ˆ σ † ˆ c ˆ b † + η ∗ ˆ σ ˆ c † ˆ b (cid:17) (14)where η ≡ η . E. An atom coupled to the quantized EM field and dressed by a classical acousticfield
For classical acoustic oscillations the operator ˆq = Q ( r ) ˆ b + Q ∗ ( r ) ˆ b † in Eq. (10) becomesa classical function q = Q ( r ) e − iΩt + Q ∗ ( r ) e iΩt (15)where Q is a coordinate-dependent complex amplitude of classical oscillations. Near theparametric resonance (cid:0) ω + Ω ≈ W (cid:126) (cid:1) the RWA Hamiltonian takes the formˆ H = ˆ H em + ˆ H a − (cid:0) R ˆ σ † ˆ ce − iΩt + R ∗ ˆ σ ˆ c † e iΩt (cid:1) . (16)where R = [ d ( Q · ∇ ) E ] r = r a . The value of the acoustic frequency Ω in Eq. (16) can be ofeither 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 theabove expression for R .Qualitatively, Hamiltonian (13) corresponds to the decay of the fermionic excitationinto a photon and phonon; Hamiltonian (14) corresponds to the decay of a photon into aphonon and fermionic excitation, whereas Hamiltonian (16) describes parametric decay of aphoton into an atomic excitation and back, mediated by classical acoustic oscillations. III. PARAMETRIC RESONANCE IN A CLOSED SYSTEM
When the system is closed and there is no dissipation, the general analytic solution tothe dynamics of coupled fermions, photons, and phonons can be obtained in the RWA. Wewrite the state vector as 8 = ∞ (cid:88) α,n =0 ( C αn | α (cid:105) | n (cid:105) | (cid:105) + C αn | α (cid:105) | n (cid:105) | (cid:105) ) . (17)Here Greek letters denote phonon states, Latin letters denote photon states, and numbers0, 1 describe fermion states. We will keep the same sequence of symbols throughout thepaper: C phonon photon fermion | phonon (cid:105) | photon (cid:105) | f ermion (cid:105) . Next, we substitute Eq. (17) into the Schr¨odinger equation, i (cid:126) ∂∂t | Ψ (cid:105) = ˆ H | Ψ (cid:105) (18)Where ˆ H is the RWA Hamiltonian. For definiteness, we consider the vicinity of the para-metric resonance with a plus sign, ω + Ω ≈ W (cid:126) , which corresponds to the Hamiltonian (13).In this case the equations for the coefficients in Eq. (17) can be separated into the pairs ofcoupled equations ddt C αn C ( α − n − + iω α,n − iΩ ( α,n ) ∗ R − iΩ ( α,n ) R iω α,n − i∆ C αn C ( α − n − = 0 , (19)and a separate equation for the lowest-energy state: · C + iω , C = 0 , (20)where Ω ( α,n ) R = η (cid:126) √ αn, ω α,n = Ω (cid:18) α + 12 (cid:19) + ω (cid:18) n + 12 (cid:19) , ∆ = Ω + ω − W (cid:126) . Note that approximate Eqs. (19),(20) preserve the norm exactly : | C | + ∞ , ∞ (cid:88) α =1 ,n =1 (cid:16) | C αn | + (cid:12)(cid:12) C ( α − n − (cid:12)(cid:12) (cid:17) = ∞ , ∞ (cid:88) α =0 ,n =0 (cid:0) | C αn | + | C αn | (cid:1) = const . The solution to Eq. (20) is trivial: C ( t ) = C (0) exp ( − iω , t ). The solution to Eqs. (19)takes the form C αn C ( α − n − = Ae − Λ ( α,n )1 t a ( α,n )1 + Be − Λ ( α,n )2 t a ( α,n )2 , (21)9here the constants A and B are determined from initial conditions. Here the eigenvalues Λ ( α,n )1 , and eigenvectors a ( α,n )1 , of the matrix of coefficients in Eq. (19) are given by Λ ( α,n )1 , = iω α,n − iδ ( α,n )1 , , a ( α,n )1 , = δ ( α,n )1 , Ω ( α,n ) ∗ R , (22)where δ ( α,n )1 , = ∆ ± (cid:114) ∆ (cid:12)(cid:12)(cid:12) Ω ( α,n ) R (cid:12)(cid:12)(cid:12) . (23) Λ ( ) Λ ( ) - - - Δ ( Ω R ( ) ) E i gen f r equen c i e s ( Ω R ( , ) ) FIG. 2: Frequency eigenvalues of the coupled electron-photon-phonon quantum system asa function of detuning from the parametric resonance W (cid:126) = Ω + ω . All frequencies are inunits of the generalized Rabi frequency Ω (1 , R . The values of eigenfrequencies are shiftedvertically by ω , | ∆ =0 .Fig. 2 shows the eigenfrequencies of the system given by Eq. (22) with α = n = 1,shifted by ω , | ∆ =0 . One can see the anticrossing with splitting by 2 Ω (1 , R at the parametricresonance.As an example, consider an exact parametric resonance, W (cid:126) = Ω + ω and the simplestinitial state Ψ = | (cid:105) | (cid:105) | (cid:105) corresponding to the initially excited atom in a cavity. In thiscase the only nonzero amplitudes are C and C : C C = 12 e − i (cid:16) ω , − (cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12)(cid:17) t e − iθ + 12 e − i (cid:16) ω , + (cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12)(cid:17) t − e − iθ , (24)where ω , = Ω (cid:18) (cid:19) + ω (cid:18) (cid:19) , Ω (1 , R = η (cid:126) = (cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12) e iθ . Ψ = e − iω , t (cid:104) ie − iθ sin (cid:16)(cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17) | (cid:105) | (cid:105) | (cid:105) + cos (cid:16)(cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17) | (cid:105) | (cid:105) | (cid:105) (cid:105) . (25)This is clearly an entangled electron-photon-phonon state, which is not surprising. In theabsence of dissipation, any coupling between these subsystems leads to entanglement.The dynamics of the corresponding physical observables, such as the energy of the fieldand the atom, is Rabi oscillations at the frequency which generalizes a standard Rabi fre-quency to the case of a parametric photon-phonon-atom resonance and which depends onboth the spatial structure of the photon and phonon fields and their occupation numbers: (cid:104) Ψ | ˆE | Ψ (cid:105) = | E ( r ) | (cid:104) − cos (cid:16) (cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17)(cid:105) (26) (cid:104) Ψ | ˆ H a | Ψ (cid:105) = W (cid:16) (cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17) IV. DYNAMICS OF AN OPEN ELECTRON-PHOTON-PHONON SYSTEMA. 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 evolutionfor the state vector, which is derived in Appendix. This is basically the Schr¨odinger equationmodified by adding a linear relaxation operator and the noise source term with appropriatecorrelation properties. The latter are related to the parameters of the relaxation operator,which is a manifestation of the fluctuation-dissipation theorem [31]. In Appendix we de-rived the general form of the stochastic equation of evolution from the Heisenberg-Langevinequations [23, 27, 28] and showed how physically reasonable constraints on the observablesdetermine the properties of the noise sources. We also demonstrated the relationship be-tween our approach and the Lindblad method of solving the master equation.11 Ω R ( ) t F i e l d i n t en s i t y (a) Ω R ( ) t A t o m ene r g y (b) FIG. 3: (a) Normalized field intensity, (cid:104) Ψ | ˆE | Ψ (cid:105) / | E ( r ) | , and (b) normalized atom energy (cid:104) Ψ | ˆ H a | Ψ (cid:105) /W as a function of time in units of the generalized Rabi frequency Ω (1 , R .Within our approach the system is described by a state vector which has a fluctuatingcomponent: | Ψ (cid:105) = | Ψ (cid:105) + (cid:102) | Ψ (cid:105) , where the straight bar means averaging over the statistics ofnoise and the wavy bar denotes the fluctuating component. This state vector is of coursevery different from the state vector obtained by solving a standard Schr¨odinger equationfor a closed system. In fact, coupling to a dissipative reservoir leads to the formation of amixed state, which can be described by a density matrix ˆ ρ = | Ψ (cid:105) · (cid:104) Ψ | + (cid:102) | Ψ (cid:105) (cid:102) (cid:104) Ψ | . One shouldview the stochastic equation approach as a convenient formalism for calculating physicalobservables.Following the derivation in Appendix, Eqs. (19),(20) are modified as12 dt C αn C ( α − n − + iω α,n + γ αn − iΩ ( α,n ) ∗ R − iΩ ( α,n ) R iω α,n − i∆ + γ ( α − n − C αn C ( α − n − = − i (cid:126) R αn R ( α − n − , (28)˙ C + ( iω , + γ ) C = − i (cid:126) R . (29)Coupling to a reservoir introduces two main differences to Eqs. (28),(29) as compared toEqs. (19),(20) for a closed system. First, eigenfrequencies acquire imaginary parts whichdescribe relaxation: ω α,n = ⇒ ω α,n − iγ αn , ω α,n − ∆ = ⇒ ω α,n − ∆ − iγ ( α − n − , ω , = ⇒ ω , − iγ . The relaxation constants are determined by the properties of all subsystems. They arederived 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 − i (cid:126) R αn , − i (cid:126) R ( α − n − and − i (cid:126) R . They are equal to 0 after averaging over the noise statistics: R αn = R ( α − n − = R . The averages of the quadratic combinations of noise source termsare nonzero. Including the noise sources is crucial for consistency of the formalism: it en-sures the conservation of the norm of the state vector and leads to a physically meaningfulequilibrium state. Note that the Weisskopf-Wigner theory does not enforce the conservationof the norm. B. Evolution of the state amplitudes and observables
The solution to Eq. (29) is C = e − ( iω , + γ ) t (cid:20) C (0) − i (cid:126) (cid:90) t e ( iω , + γ ) τ R ( τ ) dτ (cid:21) . (30)The solution to Eqs. (28) is determined again by the eigenvalues and eigenvectors of thematrix of coefficients, which are now modified by relaxation rates: Λ ( α,n )1 , = iω α,n − iδ ( α,n )1 , , a ( α,n )1 , = δ ( α,n )1 , − iγ αn Ω ( α,n ) ∗ R , (31)13here δ ( α,n )1 , = ∆ i γ αn + γ ( α − n − ± (cid:115) (cid:2) ∆ + i (cid:0) γ ( α − n − − γ αn (cid:1)(cid:3) (cid:12)(cid:12)(cid:12) Ω ( α,n ) R (cid:12)(cid:12)(cid:12) . (32)The solution to Eqs. (28) takes the form C αn C ( α − n − = e − Λ ( α,n )1 t a ( α,n )1 (cid:32) A − i (cid:126) (cid:90) t e Λ ( α,n )1 τ R αn ( τ ) a ( α,n )2 − R ( α − n − ( τ ) a ( α,n )2 − a ( α,n )1 dτ (cid:33) + e − Λ ( α,n )2 t a ( α,n )2 (cid:32) B − i (cid:126) (cid:90) t e Λ ( α,n )2 τ R ( α − n − ( τ ) − R αn ( τ ) a ( α,n )1 a ( α,n )2 − a ( α,n )1 dτ (cid:33) (33)Where the constants A and B are determined by initial conditions.As an example, we consider the reservoir at low temperatures, when the steady-statepopulation should go to the ground state | (cid:105) | (cid:105) | (cid:105) . In this case we can take γ = 0, asshown below. We will also assume that the only nonzero correlator of noise is delta-correlatedin time: R ( t + ξ ) R ∗ ( t ) = (cid:126) δ ( ξ ) D . (34)Then Eqs. (29) and (30) yield ddt | C | = D , (35)whereas Eqs. (28) give ddt (cid:16) | C αn | + (cid:12)(cid:12) C ( α − n − (cid:12)(cid:12) (cid:17) = − (cid:16) γ αn | C αn | + γ ( α − n − (cid:12)(cid:12) C ( α − n − (cid:12)(cid:12) (cid:17) . (36)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 beconserved. This gives D = 2 ∞ , ∞ (cid:88) α =1 ,n =1 (cid:16) γ αn | C αn | + γ ( α − n − (cid:12)(cid:12) C ( α − n − (cid:12)(cid:12) (cid:17) . (37)In Appendix we discuss in detail the dependence of the noise correlator on the averageddyadic components of the state vector. We also show how to find the correlators whichensure that the system approaches thermal distribution at a finite temperature.14he above formalism allows us to obtain analytic solutions to the state vector and observ-ables at any temperatures and detunings from the parametric resonance, while still withinthe RWA limits. However, the resulting expressions are very cumbersome and they arebetter to visualize in the plots. Let’s give an example of the solution at zero reservoir tem-perature and exactly at the parametric resonance W (cid:126) = Ω + ω , when the expressions aremore manageable. Consider the initial state Ψ = | (cid:105) | (cid:105) | (cid:105) when an atom is excited andboson modes are in the ground state. In this case the only nonzero amplitudes are C , C and C . To make the algebra a bit simpler, we assume that the dissipation is weakenough and its effect on the eigenvectors a ( α,n )1 , can be neglected. As a result, we obtain Ψ = e − ( iω , − γ γ ) t (cid:104) ie − iθ sin (cid:16)(cid:12)(cid:12)(cid:12) (cid:101) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17) | (cid:105) | (cid:105) | (cid:105) + cos (cid:16)(cid:12)(cid:12)(cid:12) (cid:101) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17) | (cid:105) | (cid:105) | (cid:105) (cid:105) + C | (cid:105) | (cid:105) | (cid:105) , (38)where | C | = 1 − e − ( γ + γ ) t , (cid:101) Ω (1 , R = (cid:115)(cid:12)(cid:12)(cid:12) Ω (1 , R (cid:12)(cid:12)(cid:12) − ( γ − γ ) , θ = Arg (cid:104) Ω (1 , R (cid:105) . As we see, dissipation leads not only to the relaxation of the entangled part of the statevector, but also to the frequency shift of the Rabi oscillations. This shift is absent if γ = γ .The resulting expressions for the observables, such as the EM field intensity and theenergy of the atomic excitation are (cid:104) Ψ | ˆE | Ψ (cid:105) = | E ( r ) | (cid:104) e − ( γ + γ ) t − cos (cid:16) (cid:12)(cid:12)(cid:12) (cid:101) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17) e − ( γ + γ ) t (cid:105) , (39) (cid:104) Ψ | ˆ H a | Ψ (cid:105) = W (cid:16) (cid:12)(cid:12)(cid:12) (cid:101) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17) e − ( γ + γ ) t (40)Fig. 4 illustrates the dynamics of observables in Eqs. (39) and (40).Note that the Weisskopf-Wigner theory would give the same expression (40) for theatomic energy, but a wrong expression for the EM field intensity: (cid:104) Ψ | ˆE | Ψ (cid:105) = | E ( r ) | (cid:104) − cos (cid:16) (cid:12)(cid:12)(cid:12) (cid:101) Ω (1 , R (cid:12)(cid:12)(cid:12) t (cid:17)(cid:105) e − ( γ + γ ) t , which does not approach the correct vacuum state.15 Ω R ( ) t F i e l d i n t en s i t y (a) Ω R ( ) t A t o m ene r g y (b) FIG. 4: (a) Normalized field intensity, (cid:104) Ψ | ˆE | Ψ (cid:105) / | E ( r ) | , and (b) normalized atomenergy (cid:104) Ψ | ˆ H a | Ψ (cid:105) /W as a function of time in units of the generalized Rabi frequency Ω (1 , R . Here γ + γ = 0 . Ω (1 , R . C. Emission spectra
According to [23], the power spectrum of the emission is S ( r , ν ) = 1 π Re (cid:90) ∞ dτ G (1) ( r , r ; τ ) e iντ , (41)where G (1) ( r , r ; τ ) is the field autocorrelation function at the position r of the detector: G (1) ( r , r ; τ ) = | E ( r ) | (cid:90) ∞ dt (cid:104) ˆ c † d ( t )ˆ c d ( t + τ ) (cid:105) . (42)where ˆ c d ( t ) , ˆ c † d ( t ) are annihilation and creation operators for the photons which interact withthe detector, and the Heisenberg picture is used. We will assume that the coupling between16he photons and the detector is weak, so the photon detection does not affect the dynamics ofthe intracavity photons. According to [32], ˆ c d ( t ) ∝ ˆ c ( t ) for a nanocavity, so we can calculatethe G (1) ( r , r ; τ ) using operators for the cavity field ˆ c ( t ) , ˆ c † ( t ), up to a constant factor in theresult. Note that the lower limit of the integral over t is set to be t = 0, which requires thatno photons exist before t = 0.In the Heisenberg-Langevin approach, an operator in the Heisenberg picture can be ex-pressed through Schr¨odinger’s operators using the effective Hamiltonian ˆ H eff , which con-tains the anti-Hermitian part; see the Appendix. At the same time, the inhomogeneousterm proportional to the noise sources should be added. Including these noise terms in thesolution for the field operators when calculating the emission spectra is equivalent to takinginto account the detection of thermal radiation which seeps into the cavity from outside andspontaneous emission resulting from thermal excitation of an atom. We assume that thereservoir temperature in energy units is much lower than W and (cid:126) ω , so that the contributionof these noise terms to the emission spectra can be neglected (although noise is still neededto preserve the norm).Then, the average correlator (cid:104) ˆ c † ( t )ˆ c ( t + τ ) (cid:105) is expressed as (cid:104) ˆ c † ( t )ˆ c ( t + τ ) (cid:105) = (cid:104) Ψ ( t = 0) | e i ˆ H † eff t/ (cid:126) ˆ c † e − i ˆ H eff t/ (cid:126) e i ˆ H † eff ( t + τ ) / (cid:126) ˆ ce − i ˆ H eff ( t + τ ) / (cid:126) | Ψ ( t = 0) (cid:105) = (cid:104) Ψ ( t ) | ˆ c † e − i ˆ H eff t/ (cid:126) e i ˆ H † eff ( t + τ ) / (cid:126) ˆ c | Ψ ( t + τ ) (cid:105) , (43)where | Ψ ( t ) (cid:105) is the state vector of the system which we found in the previous subsection.It can be written as | Ψ ( t ) (cid:105) = (cid:80) ∞ n =0 C n ( t ) | n (cid:105)| Ψ α,en ( t ) (cid:105) , where | Ψ α,en ( t ) (cid:105) is the part describingphonons and electrons. Therefore. (cid:104) ˆ c † ( t )ˆ c ( t + τ ) (cid:105) = (cid:32) ∞ (cid:88) n =0 C ∗ n ( t ) (cid:104) n |(cid:104) Ψ α,en ( t ) | (cid:33) ˆ c † e − i ˆ H eff t/ (cid:126) e i ˆ H † eff ( t + τ ) / (cid:126) ˆ c (cid:32) ∞ (cid:88) n =0 C n ( t + τ ) | n (cid:105)| Ψ α,en ( t + τ ) (cid:105) (cid:33) = (cid:32) ∞ (cid:88) n =0 √ nC ∗ n ( t ) (cid:104) n − |(cid:104) Ψ α,en ( t ) | (cid:33) e − i ˆ H eff t/ (cid:126) e i ˆ H † eff ( t + τ ) / (cid:126) (cid:32) ∞ (cid:88) n =0 √ nC n ( t + τ ) | n − (cid:105)| Ψ α,en ( t + τ ) (cid:105) (cid:33) . (44)Consider a simple example when the initial state is | (cid:105)| (cid:105)| (cid:105) . Within the RWA the systemcan only reach states | (cid:105)| (cid:105)| (cid:105) , | (cid:105)| (cid:105)| (cid:105) and | (cid:105)| (cid:105)| (cid:105) . After acting with ˆ c on a state of the17ystem, a new state | (cid:105)| (cid:105)| (cid:105) can also appear, but it cannot evolve into other states. So, inthis case we have (cid:104) ˆ c † ( t )ˆ c ( t + τ ) (cid:105) = ( C ∗ ( t ) (cid:104) |(cid:104) Ψ α,e ( t ) | ) e − i ˆ H eff t/ (cid:126) e i ˆ H † eff ( t + τ ) / (cid:126) ( C ( t + τ ) | (cid:105)| Ψ α,e ( t + τ ) (cid:105) )= ( C ∗ ( t ) (cid:104) |(cid:104) |(cid:104) | ) e − i ˆ H eff t/ (cid:126) e i ˆ H † eff ( t + τ ) / (cid:126) ( C ( t + τ ) | (cid:105)| (cid:105)| (cid:105) )= C ∗ ( t ) C ( t + τ ) exp[ iω , τ − γ (2 t + τ )] , (45)where we used Eqs. (A32) and (A33) and assumed that the noise for state | (cid:105)| (cid:105)| (cid:105) has zerocorrelator. Since C ( t ) = i sin (cid:16) | ˜ Ω (1 , R | t (cid:17) exp[ − iω , t − γ + γ t ] , (46)we obtain (cid:104) ˆ c † ( t )ˆ c ( t + τ ) (cid:105) = sin (cid:16) | ˜ Ω (1 , R | t (cid:17) sin (cid:16) | ˜ Ω (1 , R | ( t + τ ) (cid:17) exp[ − iωτ ] exp [ − γ ac (2 t + τ )] , (47)where we introduced the notation γ ac ≡ γ + γ + γ . Then the power spectrum is foundto be S ( r , ν ) ∝ π | E ( r ) | | ˜ Ω (1 , R | γ ac ( | ˜ Ω (1 , R | + γ ) Re (cid:34) γ ac − i ( ν − ω )[ γ ac − i ( ν − ω )] + | ˜ Ω (1 , R | (cid:35) . (48)The normalized power spectra are shown in Fig. 5 for various values of | ˜ Ω (1 , R | /γ ac . For | ˜ Ω (1 , R | < γ ac the spectrum has a single maximum at zero detuning ν = ω . For | ˜ Ω (1 , R | > γ ac the spectra are split and their maxima (same value for all spectra) are reached at detuningsgiven by ( ν − ω ) = | ˜ Ω (1 , R | − γ . Therefore, to reach the strong coupling regime the Rabifrequency | ˜ Ω (1 , R | has to exceed the combination of the decoherence rates denoted by γ ac . D. Relaxation rates
Finally, we give explicit expressions for the relaxation constants γ αn and γ αn . They werederived in Appendix using the Lindblad master equation approach and assuming statisticalindependence of “partial” dissipative reservoirs for the atomic, EM, and phonon subsystems.The result is γ αn = γ N T a + µ ω (cid:2) n T em ω ( n + 1) + (cid:0) n T em ω + 1 (cid:1) n (cid:3) + µ Ω (cid:104) n T p Ω ( α + 1) + (cid:16) n T p Ω + 1 (cid:17) α (cid:105) , (49)18
10 -5 0 5 1000.20.40.60.811.2 0.51.02.05.0
FIG. 5: The emission spectra for | ˜ Ω (1 , R | /γ ac equal to 0.5, 1, 2 and 5. All spectra arenormalized by the same constant. γ αn = γ N T a + µ ω (cid:2) n T em ω ( n + 1) + (cid:0) n T em ω + 1 (cid:1) n (cid:3) + µ Ω (cid:104) n T p Ω ( α + 1) + (cid:16) n T p Ω + 1 (cid:17) α (cid:105) , (50)where γ , µ ω and µ Ω are partial relaxation rates of the atomic, photon, and phonon sub-systems respectively; N T a = e − WTa , N T a = e − WTa e − WTa , n T em ω = e (cid:126) ωTem − , n T p Ω = e (cid:126) ΩTp − are theiroccupation numbers at thermal equilibrium; T a,em,p are temperatures of partial atom, pho-ton, and phonon reservoirs in energy units. As a reminder, the atom energy is equal to 0 instate | (cid:105) and W in state | (cid:105) .If all reservoirs are at zero temperature, we obtain γ αn = µ ω n + µ Ω α, γ αn = γ µ ω n + µ Ω α. (51)Eq. (51) shows that γ = 0, validating our choice earlier in this section. We also obtainphysically intuitive expressions for γ and γ : γ = µ ω + µ Ω , γ = γ . V. CLASSICAL ACOUSTIC PUMPING
In this case the RWA Hamiltonian is given by Eq. (16). It depends only on quantumoperators ˆ σ, ˆ σ † and ˆ c, ˆ c † ; therefore the state vector has to be expanded over the basis states | n (cid:105) | (cid:105) and | n (cid:105) | (cid:105) : Ψ = ∞ (cid:88) n =0 ( C n | n (cid:105) | (cid:105) + C n | n (cid:105) | (cid:105) ) . (52)19ubstituting Eq. (52) in the Schr¨odinger equation with the Hamiltonian (16), we again getseparation into a block of two equations,˙ C n = − iω n C n + i R ∗ (cid:126) e iΩt C ( n − √ n, (53)˙ C ( n − = − i (cid:18) ω n − + W (cid:126) (cid:19) C ( n − + i R (cid:126) e − iΩt C n √ n, (54)and a separate equation for the amplitude of the ground state | (cid:105) | (cid:105) of the system:˙ C = − iω C , (55)where ω n = ω (cid:0) n + (cid:1) . After making the substitution C ( n − = G ( n − e − iΩt , Eqs. (53),(54)give the equations similar in form to Eqs. (19): ddt C n G ( n − + iω n − iΩ ( n ) ∗ R − iΩ ( n ) R iω n − i∆ C n G ( n − = 0 , (56)where Ω ( n ) R = R (cid:126) √ n, ∆ = Ω + ω − W (cid:126) , ω n − ∆ = ω n − + W (cid:126) . Eqs. (55),(56) are different from Eqs. (19),(20) only in one aspect: they don’t contain theindex of the quantum state of the phonon field, whereas the Rabi frequency depends on theamplitude of classical acoustic oscillations Q ( r a ), see Sec. IIE. Obviously, the solution toEqs. (55),(56) will have the same form and the expressions (26), (27) for the observables willremain the same, after dropping the index of the quantum phonon state and redefining theRabi frequency.Dissipation due to coupling to a reservoir can be included using the stochastic equationof evolution of the state vector, see the Appendix. The corresponding equations are againsimilar to those for a fully quantum problem given by Eqs. (28),(29):˙ C + i ( ω + γ ) C = − i (cid:126) R , (57) ddt C n C ( n − + iω n + γ n − iΩ ( n ) ∗ R − iΩ ( n ) R iω n − i∆ + γ ( n − C n C ( n − = − i (cid:126) R n R ( n − . (58)Since the acoustic field is now a given external pumping, the relaxation constants shouldnot depend on the parameters of a phonon reservoir. They can be obtained after obvioussimplification of Eqs. (49),(50): 20 n = γ N T a + µ ω (cid:2) n T em ω ( n + 1) + (cid:0) n T em ω + 1 (cid:1) n (cid:3) , (59) γ n = γ N T a + µ ω (cid:2) n T em ω ( n + 1) + (cid:0) n T em ω + 1 (cid:1) n (cid:3) , (60)All expressions for the state vector and observables can be obtained from the correspondingexpressions in Sec. IV after dropping the index α of the quantum state of the phonon fieldand redefining the frequency of Rabi oscillations. VI. 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, thephonon frequency Ω can be much lower than the optical frequency. In this case the overlap ofthe parametric (three-wave) resonance ω ± Ω ≈ W (cid:126) and the one-photon (two-wave) resonance ω ≈ W (cid:126) can be an issue.First of all, it is clear that the resonances can be separated only if the value of Ω exceedsthe sum of the spectral widths of the EM cavity mode and the electron transition.Second, the separation criterion imposes certain restrictions on the Rabi frequencies ofthe two resonances. To derive these restrictions, we neglect dissipation and retain in theHamiltonian (11) both the RWA terms near the parametric resonance ω ± Ω ≈ W (cid:126) , and theterms near a one-photon resonance ω ≈ W (cid:126) . Since the result will be almost the same whetherthe phonon field is quantized or classical, we will consider the classical phonon field to keepthe expressions a bit shorter. The resulting Hamiltonian isˆ H = (cid:126) ω (cid:18) ˆ c † ˆ c + 12 (cid:19) + W ˆ σ † ˆ σ − (cid:0) χ + R e − iΩt (cid:1) ˆ σ † ˆ c − (cid:0) χ ∗ + R ∗ e iΩt (cid:1) ˆ σ ˆ c † , (61)where χ = ( d · E ) r = r a , R = [ d ( Q · ∇ ) E ] r = r a ; Q is now a complex-valued amplitude ofclassical 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 (cid:126) . The change of sign in Ω corresponds to replacing Q with Q ∗ in the expressionfor R .The state vector should be sought in the form of Eq. (52). After substituting it into theSchr¨odinger equation we obtain coupled equations for the amplitudes of basis states | n (cid:105) | (cid:105) ,21 n − (cid:105) | (cid:105) : ˙ C n + iω n C n − i (cid:126) (cid:0) χ ∗ + R ∗ e iΩt (cid:1) √ nC ( n − = 0 , (62)˙ C ( n − + i (cid:18) ω n − + W (cid:126) (cid:19) C ( n − − i (cid:126) (cid:0) χ + R e − iΩt (cid:1) √ nC n = 0 , (63)and ˙ C + iω C = 0 , (64)where ω n = ω (cid:0) n + (cid:1) . To compare these equations with Eqs. (53) and (54), it is convenientto assume that the system is exactly at one of the resonances and study the behaviorof the solution with increasing the detuning from another resonance. For example, weassume an exact parametric resonance ω + Ω = W (cid:126) . In this case the detuning from thetwo-wave resonance is W (cid:126) − ω = Ω . After the substitution C n = G n e − iω n t and C ( n − = G ( n − e − i ( ω n − + W (cid:126) ) t , we obtain from Eqs. (62) and (63) that˙ G n − i (cid:126) (cid:0) χ ∗ + R ∗ e iΩt (cid:1) √ nG ( n − e − i ( W (cid:126) − ω ) t = 0 , (65)˙ G ( n − − i (cid:126) (cid:0) χ + R e − iΩt (cid:1) √ nG n e i ( W (cid:126) − ω ) t = 0 . (66)If we neglect at first the perturbation of the system in the vicinity of the two-waveresonance, the solution to Eqs. (65) and (66) at χ = 0 is G n G ( n − = Ae iΩ (3) R t + Be − iΩ (3) R t − , (67)where Ω (3) R = (cid:126) R √ n is the Rabi frequency of the parametric resonance, A , and B arearbitrary constants. The state described by Eq. (67) is obviously entangled.To write the formal solution to Eqs. (65) and (66), we make another substitution ofvariables: G n ± G ( n − = G ± . The result is˙ G ± ∓ iΩ (3) R G ± = iΩ (2) ∗ R e − iΩt G n ± iΩ (2) R e iΩt G ( n − , where Ω (2) R = (cid:126) χ √ n is the Rabi frequency corresponding to the one-photon (two-wave)resonance. The solution to the last equation is G ± = 2( A, B ) e ± iΩ (3) R t + ie ± iΩ (3) R t (cid:90) t e ∓ iΩ (3) R τ (cid:104) Ω (2) ∗ R e − iΩτ G n ( τ ) ± Ω (2) R e iΩτ G ( n − ( τ ) (cid:105) dτ. (68)22onsidering the terms proportional to Ω (2) R as perturbation, we seek the solution as G n G ( n − = Ae iΩ (3) R t + Be − iΩ (3) R t − + δG n δG ( n − . To estimate the magnitude of the perturbation, we substitute Eq. (67) into Eq. (68). Aftersome algebra we obtain that under the condition Ω (3) R (cid:28) Ω the magnitude of the perturba-tion is δG n , ( n − ∼ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Ω (2) R Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) G n , ( n − , whereas if Ω (3) R ∼ Ω the magnitude of the perturbation is δG n , ( n − ∼ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Ω (2) R Ω (3) R (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) G n , ( n − . To summarize this part, if both Rabi frequencies Ω (3) R , Ω (2) R (cid:28) Ω , the two resonances canbe treated independently for any relationship between the magnitudes of Ω (3) R and Ω (2) R . If theabove inequality is violated, one can neglect one of the resonances only if its associated Rabifrequency is much lower than the Rabi frequency of another resonance. These restrictionsare obvious from qualitative physical reasoning: either the magnitudes of the Rabi splittingsare much smaller than the distance between resonances, or one of the splittings is muchweaker than another one.When the effect of the neighboring resonance is non-negligible, it can still be taken intoaccount in the solution. Indeed, consider the solution to Eqs. (62) and (63), taking intoaccount only the two-wave resonance, i.e. taking R = 0. After obvious substitutions, wearrive at C n C ( n − = Ae − i (cid:32) ω n − Ω + (cid:114) Ω + (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) (cid:33) t × − Ω + (cid:114) Ω + (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) Ω (2) R + Be − i (cid:32) ω n − + W (cid:126) + Ω − (cid:114) Ω + (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) (cid:33) t × Ω (2) R − Ω − (cid:114) Ω + (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) ; (69)In the limit Ω (cid:29) Ω (2) R we obtain C n C ( n − ≈ Ae − i ω n + | Ω (2) R | Ω t (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) Ω + Be − i ω n − + W (cid:126) − | Ω (2) R | Ω t − (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) Ω (70)23t is clear from Eq. (70) that the entanglement of states described by C n and C ( n − isdetermined by a small parameter (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) Ω , whereas at exact resonance the entanglement isalways stronger; see Eq. (67). Therefore, when Ω (2) R (cid:28) Ω , we can neglect the contribution ofthe two-photon resonance to the entanglement of states | n (cid:105) | (cid:105) and | n − (cid:105) | (cid:105) . However, itfollows 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) butreplacing the eigenfrequencies ω n and ω n − according to Eq. (70): ω n = ⇒ ω n + (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) Ω , ω n − + W (cid:126) = ⇒ ω n − + W (cid:126) − (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) Ω . (71)If Ω (3) R (cid:28) (cid:12)(cid:12)(cid:12) Ω (2) R (cid:12)(cid:12)(cid:12) Ω these shifts can be significant in order to interpret the spectra near thethree-wave parametric resonance.The same reasoning can be carried out to analyze the effect of a detuned three-waveresonance on the solution near the two-wave resonance.These results can be verified by an exact numerical solution of Eqs. (62) and (63) for giveninitial conditions. After that, we can obtain the spectra of C n and C ( n − . Since they areoscillating functions, their spectra form discrete lines at frequencies which we denote as ω osc .As an example, we select the case of n = 1, set | Ω (2) R | = | Ω (3) R | = 0 . Ω , and choose theinitial condition as C n (0) = 0 and C ( n − (0) = 1. The frequencies ω osc of the spectral linesfor C n and C ( n − are shown in Fig. 6. Their values are shifted by ω osc, = ω | ω = W/ (cid:126) . Thearea of the dot for each spectral line is proportional to the square of its amplitude. If amarker 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.As an illustration of the violation of the condition for resonance separation, we show theoscillation frequencies for | Ω (2) R | = | Ω (3) R | = 0 . Ω in Fig. 7. Here the anticrossing picture ofisolated resonances is smeared and cannot be observed. VII. CONTROL OF ENTANGLED STATES
In order to control the quantum state of the system, turn the entanglement on/off, reador write information into a qubit, or implement a logic gate, one has to vary the parametersof a system, for example the detuning from resonance, the field amplitude of the EM mode24
FIG. 6: The frequencies ω osc of the spectral lines for C n (left panel), and C ( n − (rightpanel), with n = 1, as functions of the photon frequency ω . The photon frequencies areshifted by W/ (cid:126) , and the positions of spectral lines ω osc are shifted by ω osc, = ω | ω = W/ (cid:126) .The area of a marker is proportional to the amplitude squared of the spectral line. Bothaxes are in units of Ω . The parameters are | Ω (2) R | = | Ω (3) R | = 0 . Ω , and the initial conditionis C n (0) = 0 and C ( n − (0) = 1. -3 -2 -1 0 1 2-5-4-3-2-10123 -3 -2 -1 0 1 2-5-4-3-2-10123 FIG. 7: The frequencies ω osc of the spectral lines for C n (left panel), and C ( n − (rightpanel), with n = 1. The notations are the same as in Fig. 6. The parameters are | Ω (2) R | = | Ω (3) R | = 0 . Ω , and the initial condition is C n (0) = 0 and C ( n − (0) = 1.at the atom position, or the intensity of a classical acoustic pumping. The analytic resultsobtained in previous sections can be readily generalized when the variation of a parameteris adiabatic, i.e. slower than the optical frequencies ω or W (cid:126) . Since the space is limited, thetime-dependent problem will be considered elsewhere. Here we consider just one example,25amely turning on/off of a classical acoustic pumping q = Q ( r ) e − iΩt + Q ∗ ( r ) e iΩt . For maximal control it is beneficial to place an atom at the point where E ( r = r a ) → Q · ∇ ) E r = r a is maximized. The equations of motion for quantum state amplitudeswere derived in Sec. V, see Eqs. (53)-(55).Consider an exact parametric resonance ω + Ω = W (cid:126) for simplicity, when ω n = ω n − + W (cid:126) . The solution to Eqs. (53)-(55) when the acoustic pumping is turned off is Ψ = C (0) e − iω t | (cid:105) | (cid:105) + ∞ (cid:88) n =1 (cid:16) C n (0) e − iω n t | n (cid:105) | (cid:105) + C ( n − (0) e − i ( ω n − + W (cid:126) ) t | n − (cid:105) | (cid:105) (cid:17) The solution when the acoustic pumping is turned on is Ψ = C (0) e − iω t | (cid:105) | (cid:105) + ∞ (cid:88) n =1 [ (cid:18) A n e − i (cid:12)(cid:12)(cid:12) Ω ( n ) R (cid:12)(cid:12)(cid:12) t + B n e i (cid:12)(cid:12)(cid:12) Ω ( n ) R (cid:12)(cid:12)(cid:12) t (cid:19) e − iω n t | n (cid:105) | (cid:105) + (cid:18) − A n e − i (cid:12)(cid:12)(cid:12) Ω ( n ) R (cid:12)(cid:12)(cid:12) t + B n e i (cid:12)(cid:12)(cid:12) Ω ( n ) R (cid:12)(cid:12)(cid:12) t (cid:19) e iθ − i ( ω n − + W (cid:126) ) t | n − (cid:105) | (cid:105) ] (72)Assume that the initial quantum state before the pumping was turned on was not entan-gled, for example, an atom was in an excited state and there were no photons: Ψ = e − i ( ω + W (cid:126) ) t | (cid:105) | (cid:105) . If the acoustic pumping is turned on at t = 0, the quantum state becomes entangled: Ψ = ie − iω t − iθ sin (cid:16)(cid:12)(cid:12)(cid:12) Ω (1) R (cid:12)(cid:12)(cid:12) t (cid:17) | (cid:105) | (cid:105) + e − i ( ω + W (cid:126) ) t cos (cid:16)(cid:12)(cid:12)(cid:12) Ω (1) R (cid:12)(cid:12)(cid:12) t (cid:17) | (cid:105) | (cid:105) , (73)Then the acoustic pumping can be turned off. Depending on the turnoff moment of time, onecan obtain various entangled photon-atom states, e.g. Bell states etc. The above reasoningis valid when the turn-on/off rate is slower than the optical frequencies and the detuningfrom the two-wave resonance ω = W (cid:126) . VIII. CONCLUSIONS
In conclusion, we showed how the entanglement in a system of a fermionic quantumemitter coupled to a quantized electromagnetic field in a nanocavity and quantized phononor mechanical vibrational modes emerges in the vicinity of a parametric resonance in the26ystem. We developed analytic theory describing the formation and evolution of entangledquantum states, which can be applied to a broad range of cavity quantum optomechanicsproblems and emerging nanocavity strong-coupling experiments. The model includes de-coherence effects due to coupling of the fermion, photon, and phonon subsystems to theirdissipative reservoirs within the stochastic evolution approach, which is derived from theHeisenberg-Langevin formalism. We showed that our approach provided the results forphysical observables equivalent to those obtained from the density matrix equations withthe relaxation operator in Lindblad form. We derived analytic expressions for the timeevolution of the quantum state and observables, and the emission spectra. The limit ofa classical acoustic pumping, the control of entangled states, and the interplay betweenparametric and standard two-wave resonances were discussed.
ACKNOWLEDGMENTS
This work has been supported in part by the Air Force Office for Scientific Researchthrough Grant No. FA9550-17-1-0341 and by NSF Award No. 1936276. M.T. acknowledgesthe support from RFBR Grant No. 20-02-00100, M.E. acknowledges the support fromFederal Research Center Institute of Applied Physics of the Russian Academy of Sciences(Project No. 0035-2019-004).
Appendix A: The stochastic equation of evolution for the state vector
The description of open quantum systems within the stochastic equation of evolution forthe state vector is usually formulated for a Monte-Carlo type numerical scheme, e.g. themethod of quantum jumps [23, 25]. We developed an approach suitable for analytic deriva-tions. Our stochastic equation of evolution is basically the Schr¨odinger equation modified byadding a linear relaxation operator and the noise source term with appropriate correlationproperties. The latter are related to the parameters of the relaxation operator in such away that the expressions for the statistically averaged quantities satisfy certain physicallymeaningful conditions.The protocol of introducing the relaxation operator with a corresponding noise sourceterm to the quantum dynamics is well known in the Heisenberg picture, where it is called the27eisenberg-Langevin method [23, 27, 28]. We develop a conceptually similar approach forthe Schr¨odinger equation. Here we derive the general form of the stochastic equation of evo-lution from the Heisenberg-Langevin equations and track how certain physically reasonableconstraints on the observables determine the correlation properties of the noise sources.
1. From Heisenberg-Langevin equations to the stochastic equation for the statevector
The Heisenberg-Langevin equation for the operator (cid:98) g of a certain observable quantitytakes the form [23, 27, 28] ddt ˆ g = i (cid:126) (cid:104) ˆ H, ˆ g (cid:105) + ˆ R (ˆ g ) + ˆ L g ( t ) , (A1)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 twooperators: [ˆ g , ˆ g ] = C , where C is a constant, correct Langevin sources should ensure theconservation of commutation relations at any moment of time, despite the presence of therelaxation operator in Eq. (A1); see [28–30].The group of terms i (cid:126) (cid:104) ˆ H, ˆ g (cid:105) + ˆ R (ˆ g ) can often be written as i (cid:126) (cid:104) ˆ H, ˆ g (cid:105) + ˆ R (ˆ g ) = i (cid:126) (cid:16) ˆ H † eff ˆ g − ˆ g ˆ H eff (cid:17) , (A2)where ˆ H eff is a non-Hermitian operator. For example, if the relaxation operator describesdissipation with relaxation constant γ , so that ˆ g ∝ e − γt , then ˆ H eff = ˆ H − i (cid:126) γ ˆ1, where ˆ1is a unit operator. Note that in the master equation for the density matrix the relaxationis often introduced in a conceptually similar way [25], (cid:104) ˆ H, ˆ ρ (cid:105) = ⇒ ˆ H eff ˆ ρ − ˆ ρ ˆ H † eff , which ishowever slightly different from the form used in Eq. (A2): (cid:104) ˆ H, ˆ g (cid:105) = ⇒ ˆ H † eff ˆ g − ˆ g ˆ H eff . Thedifference is because the commutator of an unknown operator with Hamiltonian enters withopposite sign in the master equation as compared to the Heisenberg equation.Now consider the transition from the Heisenberg-Langevin equation to the stochasticequation for the state vector. The key point is to assume that there exists the operatorof evolution ˆ U ( t ), which is determined not only by the system parameters but also by theproperties of the reservoir. This operator determines the evolution of the state vector: | Ψ ( t ) (cid:105) = ˆ U ( t ) | Ψ (cid:105) , (cid:104) Ψ ( t ) | = (cid:104) Ψ | ˆ U † ( t ) , (A3)28here Ψ = Ψ (0). Hereafter we will denote the operators in the Schr¨odinger picture with in-dex “s” to distinguish them from the Heisenberg operators. An observable can be calculatedas g ( t ) = (cid:104) Ψ ( t ) | ˆ g S | Ψ ( t ) (cid:105) = (cid:104) Ψ | ˆ g ( t ) | Ψ (cid:105) Which leads to ˆ g ( t ) = ˆ U † ( t ) ˆ g S ˆ U ( t ) , (A4)Since the substitution of Eqs. (A3) and (A4) into the standard Heisenberg equation leadsto the standard Schr¨odinger equation, it makes sense to apply the same procedure to theHeisenberg-Langevin equation in order to obtain the “stochastic variant” of the Schr¨odingerequation. The solution of the latter should yield the expression for an observable, g ( t ) = (cid:104) Ψ ( t ) | ˆ g S | Ψ ( t ) (cid:105) , Which is different from the standard expression by additional averaging over the noise statis-tics.Note that an open system interacting with a reservoir is generally in a mixed state andshould be described by the density matrix. We are describing the state of the system witha state vector which has a fluctuating component. For example, in a certain basis | α (cid:105) the state vector will be C α ( t ) = C α + (cid:102) C α , where the fluctuating component is denotedwith a wavy bar. The elements of the density matrix of the corresponding mixed state are ρ αβ = C α C ∗ β = C α · C ∗ β + (cid:102) C α · (cid:102) C β ∗ .The solution to the Heisenberg-Langevin equation can be expressed through the evolutionoperator ˆ U ( t ) using Eq. (A4). The noise source terms should be chosen to ensure theconservation of commutation relations at any moment of time, despite the presence of therelaxation operator. Since commutation relations between any two operators are conservedif 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,S ˆ U and ˆ H † eff = ˆ U † ˆ H † eff,S ˆ U into Eqs. (A1),(A2), and using ˆ U † ˆ U = ˆ1, we arrive at (cid:18) ddt ˆ U † − i (cid:126) ˆ U † ˆ H † eff,S (cid:19) ˆ g S ˆ U + ˆ U † ˆ g S (cid:18) ddt ˆ U + i (cid:126) ˆ H eff,S ˆ U (cid:19) = ˆ L g (A5)Next, we introduce the operator ˆ F , defined byˆ L g = 2 ˆ U † ˆ g S ˆ F (A6)29or the operator ˆ L † g , Eq. (A6) gives ˆ L † g = 2 ˆ F † (cid:16) ˆ U † ˆ g S (cid:17) † = 2 ˆ F † ˆ g † S ˆ U . Since ˆ g and ˆ g S areHermitian operators, ˆ L g has to be Hermitian too. (One can develop the Heisenberg-Langevinformalism 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 twoterms on the left-hand side of Eq. (A5) using the relationshipˆ L g = ˆ U † ˆ g S ˆ F + ˆ F † ˆ g S ˆ U (A7)Substituting the latter into Eq. (A5), we obtain (cid:18) ddt ˆ U † − i (cid:126) ˆ U † ˆ H † eff,S − ˆ F † (cid:19) ˆ g S ˆ U + ˆ U † ˆ g S (cid:18) ddt ˆ U + i (cid:126) ˆ H eff,S ˆ U − ˆ F (cid:19) = 0 . For simplicity we will assume operator ˆ H eff to be constant with time, i.e. we won’t differ-entiate between ˆ H eff and ˆ H eff,S . The last equation is satisfied for sure if ddt ˆ U = − i (cid:126) ˆ H eff ˆ U + ˆ F , ddt ˆ U † = i (cid:126) ˆ U † ˆ H † eff + ˆ F † . (A8)Multiplying Eqs. (A8) by the initial state vector | Ψ (cid:105) from the right and from the left, weobtain the stochastic equation for the state vector and its Hermitian conjugate: ddt | Ψ (cid:105) = − i (cid:126) ˆ H eff | Ψ (cid:105) − i (cid:126) | R ( t ) (cid:105) (A9) ddt (cid:104) Ψ | = i (cid:126) (cid:104) Ψ | ˆ H † eff + i (cid:126) (cid:104) R ( t ) | (A10)Where we introduced the notations i (cid:126) ˆ F | Ψ (cid:105) = ⇒ | R ( t ) (cid:105) , − i (cid:126) (cid:104) Ψ | ˆ F † = ⇒ (cid:104) R ( t ) | . We willalso need Eqs. (A9) and (A10) in a particular basis | α (cid:105) : ddt C α = − i (cid:126) (cid:88) ν (cid:16) ˆ H eff (cid:17) αν C ν − i (cid:126) R α , (A11) ddt C ∗ α = i (cid:126) (cid:88) ν C ∗ ν (cid:16) ˆ H † eff (cid:17) να + i (cid:126) R ∗ α , (A12)where R α = (cid:104) α | R (cid:105) , (cid:16) ˆ H eff (cid:17) αβ = (cid:104) α | ˆ H eff | β (cid:105) .Applying the same procedure to the standard Heisenberg equation (1) we obtain that inEqs. (A9),(A10): ˆ H eff ≡ ˆ H † eff = ˆ H and (cid:104) R ( t ) | ≡
0, which corresponds to the standardSchr¨odinger equation and its Hermitian conjugate.30ote 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 originalHeisenberg-Langevin equation (A1). We assume that the Langevin operators in the originalequation do not contradict this physically reasonable requirement.In general, statistical properties of noise that ensure certain physically meaningful re-quirements impose certain constraints on the noise source | R (cid:105) which enters the right-handside of the stochastic equation for the state vector. In particular, it is natural to requirethat the statistically averaged quantity | R (cid:105) = 0. We will also require that the noise source | R (cid:105) has the correlation properties that preserve the norm of the state vector averaged overthe reservoir statistics: (cid:104) Ψ ( t ) | Ψ ( t ) (cid:105) = 1 . (A13)
2. Noise correlator
The solution to Eqs. (A9) and (A10) can be formally written as | Ψ (cid:105) = e − i (cid:126) ˆ H eff t | Ψ (cid:105) − i (cid:126) (cid:90) t e i (cid:126) ˆ H eff ( τ − t ) | R ( τ ) (cid:105) dτ, (A14) (cid:104) Ψ | = (cid:104) Ψ | e i (cid:126) ˆ H † eff t + i (cid:126) (cid:90) t (cid:104) R ( τ ) | e − i (cid:126) ˆ H † eff ( τ − t ) dτ, (A15)In the basis | α (cid:105) , Eqs. (A14),(A15) can be transformed into C α = (cid:104) α | e − i (cid:126) ˆ H eff t | Ψ (cid:105) − i (cid:126) (cid:90) t (cid:104) α | e i (cid:126) ˆ H eff ( τ − t ) | R ( τ ) (cid:105) dτ, (A16) C ∗ α = (cid:104) Ψ | e i (cid:126) ˆ H † eff t | α (cid:105) + i (cid:126) (cid:90) t (cid:104) R ( τ ) | e − i (cid:126) ˆ H † eff ( τ − t ) | α (cid:105) dτ. (A17)In order to calculate the observables, we need to know the expressions for the averageddyadic combinations of the amplitudes. We can find them using Eqs. (A11) and (A12): ddt C α C ∗ β = − i (cid:126) (cid:88) ν (cid:16) H ( h ) αν C ν C ∗ β − C α C ∗ ν H ( h ) νβ (cid:17) − i (cid:126) (cid:88) ν (cid:16) H ( ah ) αν C ν C ∗ β + C α C ∗ ν H ( ah ) νβ (cid:17) + (cid:18) − i (cid:126) C ∗ β R α + i (cid:126) R ∗ β C α (cid:19) , (A18)Where we separated the Hermitian and anti-Hermitian components of the effective Hamil-tonian: (cid:104) α | ˆ H eff | β (cid:105) = H ( h ) αβ + H ( ah ) αβ . Substituting Eqs. (A16) and (A17) into the last term31n Eq. (A18), we obtain − i (cid:126) C ∗ β R α + i (cid:126) C α R ∗ β = 1 (cid:126) (cid:90) − t (cid:104) R ( t + ξ ) | e − i (cid:126) ˆ H † eff ξ | β (cid:105) (cid:104) α | R ( t ) (cid:105) dξ + 1 (cid:126) (cid:90) − t (cid:104) R ( t ) | β (cid:105) (cid:104) α | e i (cid:126) ˆ H eff ξ | R ( t + ξ ) (cid:105) dξ. To proceed further with analytical results, we need to evaluate these integrals. The simplestsituation is when the noise source terms are delta-correlated in time (Markovian). In thiscase only the point ξ = 0 contributes to the integrals. As a result, Eq. (A18)) is transformedto ddt C α C ∗ β = − i (cid:126) (cid:88) ν (cid:16) H ( h ) αν C ν C ∗ β − C α C ∗ ν H ( h ) νβ (cid:17) − i (cid:126) (cid:88) ν (cid:16) H ( ah ) αν C ν C ∗ β + C α C ∗ ν H ( ah ) νβ (cid:17) + D αβ , (A19)Where the correlator D αβ is defined by R ∗ β ( t + ξ ) R α ( t ) = R ∗ β ( t ) R α ( t + ξ ) = (cid:126) δ ( ξ ) D αβ (A20)The time derivative of the norm of the state vector is given by ddt (cid:88) α | C α | = − (cid:88) α (cid:34) i (cid:126) (cid:88) ν (cid:0) H ( ah ) αν C ν C ∗ α + C α C ∗ ν H ( ah ) να (cid:1) − D αα (cid:35) (A21)Clearly, the components D αα of the noise correlator need to compensate the decrease inthe norm due to the anti-Hermitian component of the effective Hamiltonian. Therefore theexpressions for H ( ah ) αβ and D αα have to be mutually consistent. This is the manifestation ofthe fluctuation-dissipation theorem [31].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 (cid:105) introducedabove depends on the initial state | Ψ (cid:105) and the evolution operator ˆ U , and these quantitiesform 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 solutionof the conventional Schr¨odinger equation for a closed system. In particular, we postulatedthe existence of the evolution operator ˆ U determined not only by the parameters of thedynamical system but also the properties of a dissipative reservoir, although we did notspecify any particular expression for ˆ U .As an example, consider a simple diagonal anti-Hermitian operator H ( ah ) αν :32 ( ah ) αν = − i (cid:126) γ α δ αν (A22)And introduce the following models:(i) Populations relax much slower than coherences (expected for condensed matter sys-tems). In this case we can choose D α (cid:54) = β = 0, D αα = 2 γ α | C α | ; within this model thepopulation at each state will be preserved.(ii) 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 theoccupation number of the ground state approaches 1, similar to the Weisskopf-Wigner model.The adequate choice of correlators is D α (cid:54) = β = 0, D αα ∝ δ αα down , γ α down = 0. The expressionfor the remaining nonzero correlator, D α down α down = (cid:88) α (cid:54) = α down γ α | C α | , (A23)Ensures the conservation of the norm: ddt (cid:88) α (cid:54) = α down | C α | = − (cid:88) α (cid:54) = α down γ α | C α | = − ddt | C α down | . This is an example of the correlator’s dependence on the state vector that we discussedbefore.
3. Comparison with the Lindblad method
One can choose the anti-Hermitian Hamiltonian H ( ah ) αβ and correlators D αβ in the stochas-tic equation of motion in such a way that Eq. (A19) for the dyadics C n C ∗ m correspond exactlyto the equations for the density matrix elements in the Lindblad approach. Indeed, the Lind-blad form of the master equation has the form [23, 25] ddt ˆ ρ = − i (cid:126) (cid:104) ˆ H, ˆ ρ (cid:105) + ˆ L ( ˆ ρ ) (A24)where ˆ L ( ˆ ρ ) is the Lindbladian:ˆ L ( ˆ ρ ) = − (cid:88) k γ k (cid:16) ˆ l † k ˆ l k ˆ ρ + ˆ ρ ˆ l † k ˆ l k − l k ˆ ρ ˆ l † k (cid:17) , (A25)Operators ˆ l k in Eq. (A25) and their number are determined by the model which describes thecoupling of the dynamical system to the reservoir. The form of the relaxation operator given33y 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 aproper steady state given by thermal equilibrium or supported by an incoherent pumping.Eq. (A24) is convenient to represent in a slightly different form: ddt ˆ ρ = − i (cid:126) (cid:16) ˆ H eff ˆ ρ − ˆ ρ ˆ H † eff (cid:17) + δ ˆ L ( ˆ ρ ) (A26)where ˆ H eff = ˆ H − i (cid:126) (cid:88) k γ k ˆ l † k ˆ l k , δ ˆ L ( ˆ ρ ) = (cid:88) k γ k ˆ l k ˆ ρ ˆ l † k . (A27)Writing the anti-Hermitian component of the Hamiltonian in Eqs. (A11),(A12) as H ( ah ) αβ = − i (cid:126) (cid:104) α | (cid:88) k γ k ˆ l † k ˆ l k | β (cid:105) , (A28)and defining the corresponding correlator of the noise source as R ∗ β ( t + ξ ) R α ( t ) = (cid:126) δ ( ξ ) D αβ , D αβ = (cid:104) α | δ ˆ L ( ˆ ρ ) | β (cid:105) ρ mn = C n C ∗ m , (A29)We obtain the solution in which averaged over noise statistics dyadics C n C ∗ m correspondexactly to the elements of the density matrix within the Lindblad method.Instead of deriving the stochastic equation of evolution of the state vector from theHeisenberg-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 ensuringthat they lead to the same observables as the solution of the density matrix equationswith the relaxation operator in Lindblad form [25, 33]. However, the demonstration ofdirect connection between the stochastic equation of evolution of the state vector and theHeisenberg-Langevin equation provides an important physical insight.
4. Relaxation rates for coupled subystems interacting with a reservoir
Whenever we have several coupled subsystems (such as electrons, photon modes, andphonons in this paper), each coupled to its reservoir, the determination of relaxation ratesof 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 partialLindbladians and obtain the total effective Hamiltonian.34onsider the Hamiltonian (11) of the system formed by a two-level electron system cou-pled to an EM mode field and dressed by a phonon field:ˆ H = ˆ H em + ˆ H a + ˆ H p + ˆ V . (A30)Here ˆ H em = (cid:126) ω (cid:0) ˆ c † ˆ c + ˆ c ˆ c † (cid:1) is the Hamiltonian for a single EM mode field, ˆ H a = W ˆ σ † ˆ σ + W ˆ σ ˆ σ † is the Hamiltonian for a two-level “atom” with energy levels W , , ˆ H p = (cid:126) Ω (cid:16) ˆ b † ˆ b + ˆ b ˆ b † (cid:17) is the Hamiltonian for a phonon mode, ˆ V = ˆ V + ˆ V the interaction Hamiltonian, where ˆ V , describe the atom-photon and atom-photon-phonon coupling, respectively:ˆ V = − (cid:0) χ ˆ σ † ˆ c + χ ∗ ˆ σ ˆ c † + χ ˆ σ ˆ c + χ ∗ ˆ σ † ˆ c † (cid:1) , ˆ V = − (cid:16) η ˆ σ † ˆ c ˆ b + η ∗ ˆ σ ˆ c † ˆ b † + η ˆ σ † ˆ c ˆ b † + η ∗ ˆ σ ˆ c † ˆ b + η ˆ σ ˆ c ˆ b + η ∗ ˆ σ † ˆ c † ˆ b † + η ˆ σ ˆ c ˆ b † + η ∗ ˆ σ † ˆ c † ˆ b (cid:17) , where χ , η , η are coupling constants defined before.Summing up the known (see e.g. [23, 25]) partial Lindbladians of two bosonic (infiniteamount of energy levels) and one fermionic (two-level) subsystems, we obtain L ( ˆ ρ ) = − γ N T a (cid:0) ˆ σ ˆ σ † ˆ ρ + ˆ ρ ˆ σ ˆ σ † − σ † ˆ ρ ˆ σ (cid:1) − γ N T a (cid:0) ˆ σ † ˆ σ ˆ ρ + ˆ ρ ˆ σ † ˆ σ − σ ˆ ρ ˆ σ † (cid:1) − µ ω n T em ω (cid:0) ˆ c ˆ c † ˆ ρ + ˆ ρ ˆ c † ˆ c − c † ˆ ρ ˆ c (cid:1) − µ ω (cid:0) n T em ω + 1 (cid:1) (cid:0) ˆ c † ˆ c ˆ ρ + ˆ ρ ˆ c ˆ c † − c ˆ ρ ˆ c † (cid:1) − µ Ω n T p Ω (cid:16) ˆ b ˆ b † ˆ ρ + ˆ ρ ˆ b † ˆ b − b † ˆ ρ ˆ b (cid:17) − µ Ω (cid:16) n T p Ω + 1 (cid:17) (cid:16) ˆ b † ˆ b ˆ ρ + ˆ ρ ˆ b ˆ b † − b ˆ ρ ˆ b † (cid:17) (A31)where γ , µ ω and µ Ω are partial relaxation rates of the subsystems, N T a , = (cid:16) e − W − W Ta (cid:17) − e − W , − W Ta , n T em ω = (cid:16) e (cid:126) ωTem − (cid:17) − , n T p Ω = (cid:16) e (cid:126) ΩTp − (cid:17) − ,T a,em,p are the temperatures of partial reservoirs. For the Lindblad master equation in theform Eq. (A26) we get ˆ H eff = ˆ H − i ˆ Γ , (A32)whereˆ Γ = (cid:126) (cid:110) γ (cid:0) N T a ˆ σ ˆ σ † + N T a ˆ σ † ˆ σ (cid:1) + µ ω (cid:2) n T em ω ˆ c ˆ c † + (cid:0) n T em ω + 1 (cid:1) ˆ c † ˆ c (cid:3) + µ Ω (cid:104) n T p Ω ˆ b ˆ b † + (cid:16) n T p Ω + 1 (cid:17) ˆ b † ˆ b (cid:105)(cid:111) . (A33)Using the effective Hamiltonian given by Eqs. (A32),(A33), we arrive at the stochasticequation for the state vector in the following form: ddt C αn = − i W + (cid:126) ω (cid:0) n + (cid:1) + (cid:126) Ω (cid:0) α + (cid:1) (cid:126) C αn − i (cid:126) (cid:104) α | (cid:104) n | (cid:104) | ˆ V | Ψ (cid:105) − γ αn C αn − i (cid:126) R αn , (A34)35 dt C αn = − i W + (cid:126) ω (cid:0) n + (cid:1) + (cid:126) Ω (cid:0) α + (cid:1) (cid:126) C αn − i (cid:126) (cid:104) α | (cid:104) n | (cid:104) | ˆ V | Ψ (cid:105) − γ αn C αn − i (cid:126) R αn , (A35)where γ αn = γ N T a + µ ω (cid:2) n T em ω ( n + 1) + (cid:0) n T em ω + 1 (cid:1) n (cid:3) + µ Ω (cid:104) n T p Ω ( α + 1) + (cid:16) n T p Ω + 1 (cid:17) α (cid:105) , (A36) γ αn = γ N T a + µ ω (cid:2) n T em ω ( n + 1) + (cid:0) n T em ω + 1 (cid:1) n (cid:3) + µ Ω (cid:104) n T p Ω ( α + 1) + (cid:16) n T p Ω + 1 (cid:17) α (cid:105) , (A37)Eqs. (A36),(A37) determine the rules of combining the “partial” relaxation rates for severalcoupled subsystems. [1] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).[2] P. Meystre, Ann. Phys. 525, 215 (2013).[3] J.-M. Pirkkalainen, S.U. Cho, F. Massel, J. Tuorila, T.T. Heikkila, P.J. Hakonen, and M.A.Sillanpaa, Nat. Comm. 6, 6981 (2015).[4] Y. Chu, P. Kharel, W. H. Renninger, L. D. Burkhart, L.Frunzio, P. T. Rakich, R. J. Schoelkopf,Science 358, 199 (2017).[5] S, Hong, R. Riedinger, I. Marinkovic, et al., Science 358, 203 (2017).[6] P. Arrangoiz-Arriola, E. A. Wollack, Z. Wang, M. Pechal, W. Jiang, T. P. McKenna, J. D.Witmer, R. Van Laer, and A.H. Safavi-Naeini, Nature 571, 537 (2019).[7] F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C.Carnegie, H. Ohadi, B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg, Science 354,726 (2016).[8] K.-D. Park, E. A. Muller, V. Kravtsov, P. M. Sass, J. Dreyer, J. M. Atkin, and M. B. Raschke,Nano Lett. 16, 479 (2016).[9] M. S. Tame, K. R. McEnery, S. K. Ozdemir, J. Lee, S. A. Maier, and M. S. Kim, Nat. Phys.9, 329 (2013).[10] F.C.B. Maia, B.T. OCallahan, A.R. Cadore, I.D. Barcelos, L.C. Campos, K. Watanabe, T.Taniguchi, C. Deneke, A. Belyanin, M.B. Raschke, and R.O. Freitas, Nano Lett. 19, 708(2019).
11] R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A. Scherman, E. Rosta, A. Demetriadou,P. Fox, O. Hess, and J. J. Baumberg, Nat. 535, 127 (2016).[12] A. Sipahigil, R. E. Evans, D. D. Sukachev, et al., Science 354, 847 (2016).[13] T. Yoshie, A. Scherer, J. Hendrickson, G. Khitrova, H. M. Gibbs, G. Rupper, C. Ell, O. B.Shchekin, and D. G. Deppe, Nature 432, 200 (2004).[14] J. P. Reithmaier, G. Sek A. Loffler, C. Hofmann, S. Kuhn, S. Reitzenstein, L. V. Keldysh, V.D. Kulakovskii, T. L. Reinecke, and A. Forchel, Nature 432, 197 (2004).[15] H. Leng, B. Szychowski, M.-C. Daniel, and M. Pelton, Nat Commun. 9, 4012 (2018).[16] O. Bitton, S. N. Gupta, and G. Haran, Nanophotonics 8, 559 (2019).[17] K.-D. Park, M. A. May, H. Leng, J. Wang, J. A. Kropp, T. Gougousi, M. Pelton, M. B.Raschke, Sci. Adv. 2019;5: eaav5931.[18] K. J. Satzinger, Y. P. Zhong, H. Chang, et al., Nature 563, 661665 (2018).[19] A. Bienfait, Y. P. Zhong, H.-S. Chang, M.-H. Chou, C. R. Conner, E. Dumur, J. Grebel, G.A. Peairs, R. G. Povey, K. J. Satzinger, and A. N. Cleland, Phys. Rev. X 10, 021055 (2020).[20] U. Delic et al., Science 367, 892 (2020).[21] N. Bloembergen, Nonlinear Optics (World Scientific, 1996).[22] Y. A. Bogoliubov and N. N. Mitropolsky, Asymptotic Methods in the Theory of NonlinearOscillations (Gordon and Breach, 1961).[23] M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, Cambridge,1997)[24] P. Forn-Diaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).[25] M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).[26] M.Tokman, Y. Wang, I. Oladyshkin, A. R. Kutayiah, and A. Belyanin, Phys. Rev. B 93,235422 (2016).[27] C. Gardiner and P. Zoller, Quantum Noise (Springer-Verlag, Berlin, Heidelberg, 2004).[28] M. Tokman, X. Yao, and A. Belyanin, Phys. Rev. Lett. 110, 077404 (2013).[29] M. Erukhimova and M. Tokman, Phys. Rev A 95, 013807 (2017).[30] M. Tokman, Z. Long, S. Almutairi, Y. Wang, V.Vdovin, M. Belkin, and A. Belyanin, APLPhotonics 4, 034403 (2019).[31] L.D. Landau, E.M. Lifshitz, Statistical Physics, Part 1 (Pergamon, Oxford, 1965).[32] K. H. Madsen and P. Lodahl, New J. of Phys. 15, 025013 (2013).
33] K. Blum, Density Matrix Theory and Applications (Springer, Heidelberg, 2012).33] K. Blum, Density Matrix Theory and Applications (Springer, Heidelberg, 2012).