Papers in Physics, vol. 8, art. 080008 (2016) Received: 6 September 2016, Accepted: 5 November 2016 Edited by: J. P. Paz Licence: Creative Commons Attribution 3.0 DOI: http://dx.doi.org/10.4279/PIP.080008 www.papersinphysics.org ISSN 1852-4249 Green’s functions technique for calculating the emission spectrum in a quantum dot-cavity system Edgar A. Gómez,1∗ J. D. Hernández-Rivero,2 Herbert Vinck-Posada3 We introduce the Green’s functions technique as an alternative theory to the quantum regression theorem formalism for calculating the two-time correlation functions in open quantum systems at the steady state. In order to investigate the potential of this theoret- ical approach, we consider a dissipative system composed of a single quantum dot inside a semiconductor cavity and the emission spectrum is computed due to the quantum dot as well as the cavity. We propose an algorithm based on the Green’s functions technique for computing the emission spectrum that can easily be adapted to more complex open quantum systems. We found that the numerical results based on the Green’s functions technique are in perfect agreement with the quantum regression theorem formalism. More- over, it allows overcoming the inherent theoretical difficulties associated with the direct application of the quantum regression theorem in open quantum systems. I. Introduction The measurement and control of light produced by quantum systems have been the focus of interest of the cavity quantum electrodynamics [1, 2]. Spe- cially, the emission of light powered by solid-state devices coupled to nanocavities is an extensive area of research due to its promising technological appli- cations, such as infrared and low-threshold lasers [3, 4], single and entangled photon sources [5, 6], as well as various applications in quantum cryptogra- phy [7] and quantum information theory [8]. Ex- periments with semiconductor quantum dots (QDs) ∗E-mail: eagomez@uniquindio.edu.co 1 Universidad del Quind́ıo, A.A. 2639, 630004 Armenia, Colombia. 2 Departamento de F́ısica, Universidade Federal de Minas Gerais, A.A. 486, 31270-901 Pampulha, Belo Horizonte, Minas Gerais, Brazil. 3 Universidad Nacional de Colombia, A.A. 055051, 111321 Bogotá, Colombia. embedded in microcavities have revealed a plethora of quantum effects and offer desirable properties for harnessing coherent quantum phenomena at the single photon level. For example, the Purcell en- hancement [9], photon antibunching [10], vacuum Rabi splitting [11] and strong light matter coupling [12]. These and many other quantum phenom- ena are being confirmed experimentally by observ- ing the power spectral density of the light (PSD) emitted by the quantum-dot cavity systems (QD- cavity). Thus, the PSD, or the so-called emission spectrum of the system, becomes the only relevant information that allows to study the properties of the light via measurements of correlation functions, as it is stated by the Wiener-Khintchine theorem [13]. In order to compute the emission spectrum of the QD-cavity systems in the framework of open quantum systems, different approaches have been elaborated from the theoretical point of view. For example, the method of thermodynamic Green’s functions has been applied to the determination of the susceptibilities and absorption spectrum of 080008-1 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. atomic systems embedded in nanocavities [14], the time-resolved photo-luminescence approach whose application allows to determine the emission spec- trum when an additional subsystem is considered, the so-called the photon reservoir [15]. These theo- retical approaches are based on several approxima- tions and therefore, they have their own limitations when they are considered in more general scenar- ios. In consequence, these methods are not used extensively. Frequently, the emission spectrum of QD-cavity systems is computed through the Quantum Regres- sion Theorem (QRT) [16–18], since it relates the evolution of mean values of observables and the two-time correlation functions. It is worth men- tioning that the QRT approach can be difficult to implement in a computer program because com- putational complexity increases significantly as the number of QDs or modes inside the cavity are being considered; more precisely, the dimensionality asso- ciated with the Hilbert space is large. In general, the QRT approach is time-consuming because it is required to solve a large system of coupled differen- tial equations and numerical instabilities that can arise. Moreover, theoretical complications related to dynamics of the operators involved can appear, as we will point out in the next section. In spite of this, the QRT approach is widely used in theoret- ical works, for example, in studies about photolu- minescence spectra of coupled light-matter systems in microcavities in the presence of a continuous and incoherent pumping [19, 20]. Also, in studies considering the relation between dynamical regimes and entanglement in QD-cavity systems [21, 22]. In the past, the Green’s functions technique (GFT) was successfully applied for calculating the emis- sion spectrum for a very simple quantum system, e.g., the micro-maser [23]. Nevertheless, this ap- proach has not been widely noticed in many sig- nificant situations in open quantum systems. The purpose of this work is to provide a simple and ef- ficient numerical method based on the GFT in or- der to overcome the inherent difficulties associated with the direct application of the QRT approach by solving the dynamics of the system in the fre- quency domain directly. This paper is structured as follows: the theoretical background of the QRT as well as the GFT are presented in section II. A concrete application of our proposed methodology for calculating the emission spectrum of the QD- cavity system is considered in section III. Moreover, for comparison purposes with the GFT, we discuss in some detail the methodology of the QRT for cal- culating the emission spectrum of the cavity. The numerical results for the emission spectrum of the quantum dot, as well as of the cavity obtained from both the GFT and the QRT, are shown in section IV. A discussion about our findings is summarized in section V. II. Theoretical background i. Quantum regression theorem One of the most important measurements when the light excites resonantly a QD-cavity system is the emission spectrum of the system. From a theo- retical point of view, it is assumed that it corre- sponds to a stationary and ergodic process which can be calculated as a PSD of light by using the well-known Wiener-Khintchine theorem [13]. This theorem states that the emission spectrum is given by the Fourier Transform of the two-time correla- tion function of the operator field â; explicitly, it is S(ω) = Re lim t→∞ ∫ ∞ 0 〈â†(t + τ)â(t)〉eiωτdτ. (1) In order to calculate the two-time correlation function involved in Eq. (1), a theoretical ap- proach based on the QRT is frequently consid- ered. It states that if a set of operators {Ôi(t + τ)} satisfy the dynamical equations d dτ 〈Ôi(t + τ)〉 = ∑ j Lij〈Ôj(t + τ)〉 then d dτ 〈Ôi(t + τ)Ô(t)〉 =∑ j Lij〈Ôj(t+τ)Ô(t)〉 is valid for any operator Ô(t) at an arbitrary time t. Here, Lij represents the matrix of coefficients associated with the coupled linear equations of motion. It is worth mention- ing that validity of this theorem holds whenever a closed set of operators is associated with the dy- namics. In general, to obtain the closed set of op- erators can be difficult or an impossible task, since there must be added as many operators as neces- sary in order to close the dynamics of the system. For example, to calculate the emission spectrum in a model of QD-cavity system [20, 21], two new op- erators are required because the field operators in the interaction picture do not lead to a complete set. 080008-2 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. ii. Green’s functions technique Let us consider a QD-cavity system and an operator  which does not operate on the reservoirs, then its single-time expectation value in the Heisenberg representation is given by 〈 ˆ̃A(t + τ)〉 = TrS⊗R[ ˆ̃ A(t + τ) ˆ̃ρS⊗R(t)]. (2) The density operator system-reservoir can be evolved from an initial state at time 0 to an arbi- trary time t via ˆ̃ρS⊗R(t) = Û †(t, 0)ρ̂S⊗R(0)Û(t, 0), with Û(t, 0) being a unitary time-evolution opera- tor involving the Hamiltonian terms of the system and reservoirs. Moreover, the operator ˆ̃ρS⊗R(t) = ˆ̃ρS(t) ⊗ ˆ̃ρR(t) depicts the composite density opera- tor of the system and reservoir. It is worth point- ing out that tilde means that the operator has been transformed to the Heisenberg representation and that the dynamics of the system depend directly on ˆ̃ρS⊗R(t) for all times. The validity of the Marko- vian approximation requires that the state of the system is sufficiently well described when it is con- sidered that ˆ̃ρS(t) = TrR( ˆ̃ρS⊗R(t)). Therefore, it is sufficient to write ˆ̃ρS⊗R(t) = ˆ̃ρS(t) ⊗ ˆ̃ρR(t) for all times. If we assume that at t = 0 the ini- tial state of the system is the steady state, then ˆ̃ρS⊗R(0) = ρ̂ (ss) S⊗R. Here, the superscript ”(ss)” should be understood to be the steady state of the system-reservoir. After tracing over degrees of free- dom of the reservoirs, we have that the Eq. (2) takes the form 〈 ˆ̃A(τ)〉 = TrS[ ˆ̃ A(0) ˆ̃ρS(τ)], (3) where the reduced density operator for the system is given by ˆ̃ρS(τ) = TrR[Û(τ, 0) ˆ̃ρS⊗R(0)Û †(τ, 0)] and ˆ̃ A(0) = Â. If ˆ̃ρS(τ) satisfies the Lindblad mas- ter equation dˆ̃ρS(τ)/dτ = Lˆ̃ρS(τ) with L the su- peroperator defined as Lˆ̃ρS(τ) = −i[ĤS, ˆ̃ρS(τ)] +∑ j Γj 2 (2X̂j ˆ̃ρS(τ)X̂ † j−X̂ † jX̂j ˆ̃ρS(τ)−ˆ̃ρS(τ)X̂ † jX̂j) for an operator X̂j, then the expectation value 〈 ˆ̃ A(τ)〉 can be computed by solving the dynamics associ- ated to the Lindblad master equation. It is worth mentioning that the Hamiltonian operator ĤS de- scribes the QD-cavity system and Γj corresponds to the damping (pumping) rate associated to the operator X̂j. In order to calculate the two-time correlation func- tion 〈 ˆ̃A(t + τ) ˆ̃B(t)〉 where ˆ̃A(t + τ) = Û†(t + τ,t) ˆ̃ A(t)Û(t + τ,t) and ˆ̃ B(t) = Û†(t, 0)B̂Û(t, 0) are arbitrary Heisenberg operators which do not oper- ate on the reservoirs, we proceed similarly to the case of the single-time expectation value, it is 〈 ˆ̃A(τ) ˆ̃B(0)〉 = TrS⊗R[ ˆ̃ A(τ) ˆ̃ B(0) ˆ̃ρS⊗R(0)], = TrS[ ˆ̃ A(0) ˆ̃ G(τ)], (4) where we have used the well-known properties of the unitary time-evolution operator and the fact that the system at time t = 0 is in the steady state. We have defined the operator ˆ̃ G(τ) = TrR[Û(τ, 0) ˆ̃ B(0) ˆ̃ρS⊗R(0)Û †(τ, 0)] (5) where the trace operation is performed on the reser- voirs only. By performing the time derivation of Eq. (4), we have that d dτ 〈 ˆ̃A(τ) ˆ̃B(0)〉 = TrS[ ˆ̃ A(0) d ˆ̃ G(τ) dτ ]. (6) where d ˆ̃ G(τ) dτ = d dτ TrR[Û(τ, 0) ˆ̃ B(0) ˆ̃ρS⊗R(0)Û †(τ, 0)], = d dτ TrR[ ˆ̃ B(τ) ˆ̃ρS⊗R(τ)], = TrR[ ˆ̃ B(τ) dˆ̃ρS⊗R(τ) dτ + ˆ̃ρS⊗R(τ) d ˆ̃ B(τ) dτ ], = TrR[ ˆ̃ B(τ) dˆ̃ρS⊗R(τ) dτ ], + TrR[ d ˆ̃ B(τ) dτ ˆ̃ρS⊗R(τ)]. (7) Notice that the last term vanishes since TrR[ d ˆ̃ B(τ) dτ ˆ̃ρS⊗R(τ)] = d dτ TrR[ ˆ̃ B(0) ˆ̃ρS⊗R(0)] is independent of time τ. Thus, the Eq. (7) can be reduced to the form d ˆ̃ G(τ) dτ = TrR[ dˆ̃ρS⊗R(τ) dτ ˆ̃ B(τ)], = TrR[Lˆ̃ρS⊗R(τ) ˆ̃ B(τ)], = LTrR[ ˆ̃ρS⊗R(τ) ˆ̃ B(τ)], = LTrR[Û(τ, 0) ˆ̃ B(0) ˆ̃ρS⊗R(0)Û †(τ, 0)], = L ˆ̃G(τ) (8) 080008-3 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. where we have taken into account that the super- operator L acts only on the system Hilbert space and not on the reservoir. It is straightforward to conclude that ˆ̃ G(τ) is an operator that obeys the same dynamical equations as ˆ̃ρS(τ). More pre- cisely, d ˆ̃ G(τ)/dτ = L ˆ̃G(τ) with the boundary con- dition ˆ̃ G(0) = ˆ̃ B(0) ˆ̃ρS(0) at the steady state of the system. We also conclude that the two-time correlation function in the long-time limit can be written as lim t→∞ 〈 ˆ̃A(t + τ) ˆ̃B(t)〉 = TrS[ ˆ̃ G(τ)], (9) where ˆ̃ G(τ) = TrR[Û(τ)B̂ρ̂ (ss) S⊗RÛ †(τ)] is defined as the Green’s functions operator and the operators Â, B̂ and ρ̂ (ss) S⊗R are in the Schrödinger representation.  = ˆ̃ A(0) and B̂ = ˆ̃ B(0) are operators considered at the steady state of the system. In the remainder of the paper, we assume that Û(τ) ≡ Û(τ, 0). Par- ticularly, the Eq. (9) takes the form of the Eq. (1) after performing the integral transformation. More precisely, by taking the real part of the Laplace transform to the Eq. (9), we obtain an expression in terms of the Green’s functions operator in the frequency domain as follows S(ω) = Re TrS[ ˆ̃ G(iω)]. (10) Notice that the operators  and B̂ should be defined appropriately for describing the emission spectrum due to the cavity or the quantum dot. Moreover, the wide tilde is used to indicate that the Laplace transform was taken. The superscript ”(ss)” should be understood to be the steady state of the reduced density operator of the system. Af- ter taking the Laplace transform of the Eq. (9), we obtain an expression for the emission spectrum of the system in terms of the Green’s functions oper- ator in the frequency domain as follows S(ω) = 1 πnc ReTrS[ ˆ̃ G(iω)]. (11) Prior to leaving this section, we mention that this result will be the starting point for calculating the emission spectrum due to the cavity as well as the quantum dot by considering the photon and fermionic operators in a separated way. III. Application to the QD-cavity system i. Model In order to illustrate the potential of the Green’s function technique for calculating the emission spectrum in a QD-cavity system, we will consider an open quantum system composed of a quantum dot interacting with a confined mode of the electro- magnetic field inside a semiconductor cavity. This quantum system is well described by the Jaynes- Cummings Hamiltonian [24] ĤS = ωXσ̂ †σ̂ + (ωX −∆)â†â + g(σ̂↠+ âσ̂†), (12) where the quantum dot is described as a fermionic system with only two possible states, e.g., |G〉 and |X〉 are the ground and excited state. σ̂ = |G〉〈X| and â (σ̂† = |X〉〈G| and â†) are the annihilation (creation) operators for the fermionic system and the cavity mode, respectively. The parameter g is the light-matter coupling constant. Moreover, note that we have set h̄ = 1. We also define the detun- ing between frequencies of the quantum dot and the cavity mode as ∆ = ωX − ωa, where ωX and ωa are the energies associated to an exciton and the photons inside the cavity, respectively. This Hamiltonian system is far from describing any real physical situation since it is completely integrable [25] and no measurements could be done since the light remains always inside the cavity. In order to incorporate the effects of the environ- ment on the dynamics of the system, we consider the usual approach to model an open quantum sys- tem by considering a whole system-reservoir Hamil- tonian which is frequently split in three parts. Namely, Ĥ = ĤS + ĤSR + ĤR, where ĤS de- fines the Hamiltonian term of the QD-cavity sys- tem as it is defined in the Eq. (12). The Hamilto- nian terms ĤSR and ĤR corresponding to a bilin- ear coupling between the system-reservoir and its respective reservoirs ĤR have been discussed in de- tail by Perea et al. in [26]. The reader can find a detailed discussion of the Markovian master equa- tion in [27, 28]. In the framework of open quan- tum systems, different reservoirs have been pro- posed in order to describe the dissipation, deco- herence or decays. Particularly, for QD-cavity sys- tems, a reservoir is considered for describing the 080008-4 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. physical situation where the photons are absorbed in a semiconductor and electron-hole pairs (exci- tons) can be produced which can be associated to either electrical injection or the capture of excitons optically created at frequencies larger than the typ- ical ones of our system. This process corresponds to the so-called continuous and incoherent pumping of the QD. Also, when the excitons are coupled to the leaky modes of the cavity with energy different than the cavity mode, there is a residual density of states inside the cavity and this process is responsi- ble for the spontaneous emission (radiative recom- bination) to an independent reservoir of photons. Another physical process is known as the coherent emission and it is due to the direct dissipation of the cavity mode, more precisely, the cavity mode is cou- pled to the continuum of photonic modes out of the cavity. For obtaining the master equation for the QD-cavity system, it is convenient to consider the interaction picture with respect to ĤS + ĤSR and assume the validity of the Born-Markov approxi- mation. After tracing out the degrees of freedom of all the reservoirs, one arrives to the Lindblad master equation for the reduced density matrix of the system dˆ̃ρS dτ = −i [ ĤS, ˆ̃ρS ] + κ 2 (2âˆ̃ρSâ † − â†âˆ̃ρS − ˆ̃ρSâ†â) + γ 2 (2σ̂ ˆ̃ρSσ̂ † − σ̂†σ̂ ˆ̃ρS − ˆ̃ρSσ̂†σ̂) + P 2 (2σ̂† ˆ̃ρSσ̂ − σ̂σ̂† ˆ̃ρS − ˆ̃ρSσ̂σ̂†). (13) The parameter γ is the decay rate due to the spon- taneous emission, κ is the decay rate of the cavity photons across the cavity mirrors, and P is the rate at which the excitons are being pumped. Figure 1 shows a scheme of the simplified model of the QD- cavity system showing the processes of continuous pumping P and cavity loses κ. The physical pro- cess begins when the light from the pumping laser enters into the cavity and excites one of the quan- tum dots in the QD layer. Thus, light from this source couples to the cavity and a fraction of pho- tons escapes through the partly transparent mirror from the cavity and goes to the spectrometer for measurements of the emission spectrum. A general approach for solving the dynamics of the system consists in writing the corresponding Bloch equations for the reduced density matrix of the sys- tem in the bared basis. It is an extended Hilbert space formed by taking the tensor product of the state vectors for each of the system components, {|G〉, |X〉}⊗ {|n〉}∞n=0. In this basis, the reduced density matrix ρ̂S can be written in terms of its ma- trix elements as ρ̃Sαn,βm ≡〈αn|ˆ̃ρS(τ)|βm〉. Hence, the Eq. (13) explicitly reads dρ̃Sαn,βm dτ = i [ (ωX − ∆)(m−n)ρ̃Sαn,βm + ωX(δβXρ̃Sαn,Xm − δαXρ̃SXn,βm) ] + ig [(√ m + 1δβXρ̃Sαn,Gm+1 + √ mδβGρ̃Sαn,Xm−1 ) − (√ nδαGρ̃SXn−1,βm + √ n + 1δαXρ̃SGn+1,βm )] + κ 2 ( 2 √ (m + 1)(n + 1)ρ̃Sαn+1,βm+1 − (n + m)ρ̃Sαn,βm ) − γ 2 ( δαXρ̃SXn,βm − 2δαGδβGρ̃SXn,Xm + δβXρ̃Sαn,Xm ) + P 2 ( 2δαXδβXρ̃SGn,Gm − δαGρ̃SGn,βm − δβGρ̃Sαn,Gm ) . (14) Note that we use the convention that all indices written in Greek alphabet are used for the fermionic states and take only two possible values |G〉, |X〉. The indices written in Latin alphabet are used for the Fock states which take the possible val- ues 0, 1, 2, 3 . . . Additionally, it is worth mentioning that our proposed method does not require to solve a system of coupled differential equations, instead of it, we solve a reduced set of algebraic equations that speed up the numerical solution. Prior to leaving this section, we point out that the number of excitations of the system is defined by the operator N̂ = â†â + σ̂†σ̂. The closed system and the number of excitations of the system is con- served, i.e., [ĤS,N̂] = 0. It allows us to orga- nize the states of the system through the number of excitations criterion such that the density ma- trix elements ρ̃Gn,Gn, ρ̃Xn−1,Xn−1, ρ̃Gn,Xn−1 and ρ̃Xn−1,Gn are related by having the same number of 080008-5 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. Figure 1: The picture represents a QD-cavity sys- tem showing the processes of continuous pumping P and cavity loses κ. Figure 2: Ladder of bared states for a two-level quantum dot coupled to a single cavity mode. The double headed green arrow depicts the matter cou- pling constant g, dashed red lines the emission of the cavity mode κ, solid black lines the exciton pumping rate P and solid blue lines the sponta- neous emission rate γ. quanta; sub-spaces of a fixed number of excitation evolve independently from each other. The Fig. 2 shows a schematic representation of the action of the dissipative processes involved in the dynamics of the system according to the excitation number (Nexc). ii. Emission spectrum of the cavity based on the GFT In order to compute the emission spectrum of the cavity, we will consider the two-time correlation function accordingly with the Eq. (9) for the field operator as follows lim t→∞ 〈ˆ̃a†(t + τ)ˆ̃a(t)〉 = TrS[↠ˆ̃ G(τ)]. (15) Where we have considered that the field operator is given by ˆ̃a(0) = â at the steady state. After per- forming the partial trace over the degrees of free- dom of the system, we have that TrS[â † ˆ̃G(τ)] = ∑ α,β,γ,l,m,n √ (l + 1)(m + 1) × TrR[Uαl,βm(τ) × 〈βm + 1|ρ̂(ss)S⊗R|γn〉U † γn,αl+1(τ)], (16) where the matrix elements for the time evolution operator are given by Uαl,βm(τ) = 〈αl|Û(τ)|βm〉 and U † γn,αl+1(τ) = 〈γn|Û †(τ)|αl + 1〉. In what follows, we assume the validity of the Markovian approximation, it means that the correlations be- tween the system and the reservoir must be unim- portant even at the steady state. Thus, the den- sity operator system-reservoir can be written as ρ̂ (ss) S⊗R = ρ̂ (ss) S ⊗ ρ̂ (ss) R which implies that 〈βm + 1|ρ̂(ss)S⊗R|γn〉 = ρ̂ (ss) R 〈βm + 1|ρ̂ (ss) S |γn〉. (17) Replacing the previous expression in Eq. (16), it is straightforward to show that the two-time correla- tion function reads TrS[â † ˆ̃G(τ)] = ∑ αl √ l + 1〈αl| ˆ̃G(τ)|αl + 1〉, (18) where the Green’s functions operator ˆ̃ G(τ) is given by ˆ̃ G(τ) = TrR [ Û(τ)ρ̂ (ss) R ∑ βγmn (√ m + 1 |βm〉〈γn| × 〈βm + 1| ρ̂(ss)S |γn〉 ) Û†(τ) ] . (19) 080008-6 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. As we pointed out in section II, this operator must obey the same master equation as the reduced den- sity operator of the system. In fact, the terms that only contribute in the Eq. (18) are given by the ma- trix elements G̃βm,γn(τ) ≡ 〈βm| ˆ̃ G(τ) |γn〉 of the Green’s functions operator. This is due to the fact that the projection operator |βm〉〈γn| enters into ˆ̃ G(τ) in the same way as into the reduced density operator of the system. In order to identify these matrix elements, it should be considered that for the QD-cavity system, the dynamics of all coherences asymptotically van- ish and there only remains the reduced density matrix elements which are ruled by the number of excitations criterion, i.e., ρGn,Gn, ρXn−1,Xn−1, ρGn,Xn−1, ρXn−1,Gn. Thus, the Eq. (17) can be written as follows 〈βm + 1|ρ̂(ss)S⊗R|γn〉 = ρ̂ (ss) R ( δβGδγGδm+1,n + δβXδγXδm,n−1 + δβGδγXδm,n + δβXδγGδm+1,n−1 ) × ρ(ss)Sβm+1,γn. (20) By replacing the Eq. (20) into Eq. (19), we find that the Green’s functions operator explicitly reads ˆ̃ G(τ) = TrR [ Û(τ)ρ̂ (ss) R ∑ m √ m + 1 × ( |Gm〉〈Gm + 1|ρ(ss)SGm+1,Gm+1 + |Xm〉〈Xm + 1|ρ(ss)SXm+1,Xm+1 + |Gm〉〈Xm|ρ(ss)SGm+1,Xm + |Xm〉〈Gm + 2| × ρ(ss)SXm+1,Gm+2 ) Û†(τ) ] . (21) Note that from this expression, it is easy to identify the nonzero matrix elements of the Green’s func- tions operator that contribute to the emission spec- trum. Finally, after performing the Laplace trans- form to the Eq. (18), we have that the emission spectrum of the cavity is given by S(ω) = Re ∑ l √ l + 1 ( G̃Gl,Gl+1(iω) + G̃Xl,Xl+1(iω) + G̃Gl,Xl(iω) + G̃Xl,Gl+2(iω) ) . (22) It is worth mentioning that the initial conditions may be obtained by evaluating the Green’s func- tion operator at τ = 0. Moreover, by using the fact that the time evolution operators become the iden- tity and TrR[ρ̂ (ss) R ] = 1. We obtain a set of initial conditions given by G̃Gl,Gl+1(0) = √ l + 1ρ (ss) SGl+1,Gl+1, G̃Xl,Xl+1(0) = √ l + 1ρ (ss) SXl+1,Xl+1, G̃Gl,Xl(0) = √ l + 1ρ (ss) SGl+1,Xl, G̃Xl,Gl+2(0) = √ l + 1ρ (ss) SXl+1,Gl+2. (23) Note that this set of initial conditions corresponds to the steady state of the reduced density matrix of the system. A general algorithm based on the GFT for computing the emission spectrum is presented in the appendix. We mention that this approach can be adapted easily for calculating the emission spectrum due to the cavity as well as the quantum dot. iii. Emission spectrum of the quantum dot based on the GFT In order to compute the emission spectrum of the quantum dot, we will consider the two-time corre- lation function given by Eq. (9), but for the case of the matter operator lim t→∞ 〈ˆ̃σ†(t + τ)ˆ̃σ(t)〉 = TrS[σ̂† ˆ̃ G(τ)] (24) where we have considered that the matter opera- tor is given by ˆ̃σ(0) = σ̂ at the steady state. It is straightforward to show, after performing the par- tial trace over the degrees of freedom of the system, that the two-time correlation function reads TrS[σ̂ † ˆ̃G(τ)] = ∑ αl δαX〈Gl| ˆ̃ G(τ)|αl〉, (25) 080008-7 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. where the Green’s functions operator ˆ̃ G(τ) is given by ˆ̃ G(τ) = TrR [ Û(τ) ∑ βγmn ( δβX |Gm〉〈γm| × 〈βm| ρ̂(ss)S⊗R |γn〉 ) Û†(τ) ] . (26) Assuming again the validity of the Markovian ap- proximation and taking into account the number of excitations criterion, we have that the density operator system-reservoir can be written as 〈βm|ρ̂(ss)S⊗R|γn〉 = ρ̂ (ss) R ( δβGδγGδm,n + δβXδγXδm,n + δβGδγXδm,n+1 + δβXδγGδm,n−1 ) ρ (ss) Sβm,γn. (27) By inserting the Eq. (27) into Eq. (26), we find that the Green’s functions operator explicitly reads ˆ̃ G(τ) = TrR [ Û(τ)ρ̂ (ss) R × ∑ m ( |Gm〉〈Xm|rho(ss)SXm,Xm + |Gm〉〈Gm + 1|ρ(ss)SXm,Gm+1 ) × Û†(τ) ] . (28) Analogously as in section ii, we identify the nonzero matrix elements of the Green’s functions operator that contribute to the emission spectrum. After performing the Laplace transform to Eq. (25), the emission spectrum of the quantum dot is given by S(ω) = Re ∑ l ( G̃Gl,Xl(iω)+G̃Gl,Gl+1(iω) ) . (29) Taking into account that the initial conditions are obtained by evaluating the Green’s function opera- tor at τ = 0, we have the time evolution operators become the identity and TrR[ρ̂ (ss) R ] = 1, thus, we obtain a set of initial conditions given by G̃Gl,Xl(0) = ρ (ss) SXl,Xl, G̃Gl,Gl+1(0) = ρ (ss) SXl,Gl+1, G̃Xl,Xl+1(0) = 0, G̃Gl+2,Xl(0) = 0. (30) at the steady-state. iv. Emission spectrum of the cavity based on the QRT In what follows, we apply the QRT approach for the model of QD-cavity system described in sec- tion III. In order to compute the emission spectrum of the cavity, the knowledge of the two-time corre- lation function for the field operator is required, it is 〈â†(τ)â(0)〉 in concordance with the Eq. (1). Moreover, by following the approach presented in Ref. [26], it is straightforward to show that the two- time correlation function is given by 〈â†(τ)â(0)〉 = ∑ n √ n + 1 ( 〈â†Gn(τ)â(0)〉 + 〈â†Xn(τ)â(0)〉 ) , (31) where the following definitions have been used â † Gn = |Gn + 1〉〈Gn| , â † Xn = |Xn + 1〉〈Xn| , σ̂†n = |Xn〉〈Gn| , ζ̂n = |Gn + 1〉〈Xn− 1| . (32) It is worth mentioning that the last two operators should be added in order to close the dynamics of the system accordingly to the QRT as pointed out in section II. More precisely, we are interested in solving the dynamical equations associated to the expectation values 〈â†Gn(τ)â(0)〉 and 〈â † Xn(τ)â(0)〉 as a function of τ. Therefore, we need to solve a set of coupled differential equations given by d dτ 〈â†Gn(τ)â(0)〉 = ∑ j Lij〈â † Gn(τ)â(0)〉, d dτ 〈â†Xn(τ)â(0)〉 = ∑ j Lij〈â † Xn(τ)â(0)〉, d dτ 〈σ̂†n(τ)â(0)〉 = ∑ j Lij〈σ̂†n(τ)â(0)〉, d dτ 〈ζ̂n(τ)â(0)〉 = ∑ j Lij〈ζ̂n(τ)â(0)〉. (33) In order to find explicitly the set of dynamical equations and its corresponding initial conditions, 080008-8 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. we obtain first the set of differential equations for the single-time expectation values for the operators given by Eq. (32), it is explicitly d dτ 〈â†Gn(τ)〉 = ( −P − i∆ −nκ− κ 2 + iωX ) × 〈â†Gn(τ)〉 + κ √ (n + 1)(n + 2) × 〈â†Gn+1(τ)〉 + γ〈â † Xn(τ)〉 − ig √ n〈ζ̂n(τ)〉 + ig √ n + 1 × 〈σ̂†n(τ)〉, d dτ 〈σ̂†n(τ)〉 = ig √ n + 1〈â†Gn(τ)〉 − ig √ n〈â†Xn−1(τ)〉 + 1 2 (−P −γ − 2nκ + 2iωX) × 〈σ̂†n(τ)〉 + (n + 1)κ〈σ̂†n+1(τ)〉, d dτ 〈â†Xn−1(τ)〉 = P〈â † Gn−1(τ)〉 + κ √ n(n + 1) × 〈â†Xn(τ)〉 + (−γ − i∆ −nκ + κ 2 + iωX) × 〈â†Xn−1(τ)〉 + ig √ n + 1〈ζ̂n(τ)〉 − ig √ n〈σ̂†n(τ)〉, d dτ 〈ζ̂n(τ)〉 = −ig √ n〈â†Gn(τ)〉 + ig √ n + 1 × 〈â†Xn−1(τ)〉 + (− P 2 − γ 2 − 2i∆ −nκ + iωX) × 〈ζ̂n(τ)〉 + √ n(n + 2)κ〈ζ̂n+1(τ)〉. (34) The QRT approach implies that the follow- ing two-time correlation functions 〈â†Gn(τ)â(0)〉, 〈â†Xn(τ)â(0)〉, 〈σ̂ † n(τ)â(0)〉 and 〈ζ̂n(τ)â(0)〉 satisfy the same dynamical equations given by Eq. (34), subject to the initial conditions 〈â†Gn(0)â(0)〉 = √ n + 1ρGn+1,Gn+1(0), 〈â†Xn(0)â(0)〉 = √ n + 1ρXn+1,Xn+1(0), 〈σ̂†n(0)â(0)〉 = √ n + 1ρGn+1,Xn(0), 〈ζ̂n(0)â(0)〉 = √ nρXn,Gn+1(0). (35) 995 997 999 1001 1003 1005 ω(meV ) 0 0.4 0.8 1.2 1.6 S (ω ) a Figure 3: Emission spectrum of the cavity based on the GFT as a solid blue line and the correspond- ing numerical calculation based on the QRT as a dashed red line. The parameters values are g = 1 meV, γ = 0.005 meV, κ = 0.2 meV, P = 0.3 meV, ∆ = 2 meV, ωa = 1000 meV. More precisely, it is done explicitly by perform- ing the following replacements: 〈â†Gn(τ)〉 → 〈â†Gn(τ)â(0)〉, 〈â † Xn(τ)〉 → 〈â † Xn(τ)â(0)〉, 〈σ̂†n(τ)〉 → 〈σ̂†n(τ)â(0)〉 and 〈ζ̂n(τ)〉 → 〈ζ̂n(τ)â(0)〉. The parameters of the system ωX, ∆,g,κ,γ de- termine the dynamics of the two-time correlation function, as well as setting the initial conditions that will be propagated according to the dynam- ical equations given by Eq. (34). Since we are interested in the light that the quantum system emits, we have considered the steady state of the system as the initial state into equation Eq. (35). IV. Results and discussion In this section, we compare the numerical calcula- tions based on the GFT and the QRT approach for the emission spectrum of the cavity as well as the quantum dot. In particular, the QD-cavity sys- tem can display two different dynamical regimes by changing the parameters of the system and two regimes can be achieved when the loss and 080008-9 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. pump rates are modified. In fact, the relation g > |κ−γ|/4 holds for the strong coupling regime and the relation g < |κ−γ|/4 remains valid for the weak coupling regime. Figure 3 shows the numeri- cal results for the emission spectrum associated to the cavity in the strong coupling regime, where the emission spectrum of the cavity based on the GFT is shown as a solid blue line and the emission spec- trum based on the QRT approach as a dashed red line. The parameters of the system are g = 1 meV, γ = 0.005 meV, κ = 0.2 meV, P = 0.3 meV, ∆ = 2 meV and ωa = 1000 meV. Particularly, for this set of parameters values, we can identify two different peaks which are associated to the energy of the cav- ity and the quantum dot, they are ωa ≈ 998.3 meV and ωX ≈ 1000.3 meV. We have considered the rel- ative error as a quantitative measure of the discrep- ancy between the GFT and the QRT approaches. More precisely, by monitoring the numerical com- putations of the emission spectrum, we have esti- mated that the maximum relative error is on the order of 10−3 in all numerical calculations that we have performed. Similarly, the emission spectrum of the cavity based on the GFT (solid blue line) and the QRT approach (dashed red line) for the strong coupling regime are shown in Fig. 4. Here, the parameters values of the system are given by g = 1 meV, γ = 0.005 meV, κ = 2 meV, P = 0.005 meV, ∆ = 0.0 meV and ωa = 1000 meV. Note that we have considered the resonant case, more precisely, the same energy values for the cavity and the quantum dot. Here, the emission spectrums do not match but repel each other, resulting in a struc- ture of two separate peaks for a distance of approxi- mately two times the coupling constant, i.e., 2g ≈ 2 meV. It is worth mentioning that this quantum ef- fect is well-known as Rabi splitting in QD-cavity systems. The emission spectrum of the quantum dot in the weak coupling regime is shown in Fig. 5. The numerical result for the emission spectrum of the quantum dot based on the GFT is shown as a solid blue line and the corresponding numerical result for QRT approach is shown as a dashed red line. We set the weak coupling regime by consider- ing high values of the decay and pump rates κ = 5 meV and P = 1 meV, respectively. The rest of the parameters values are g = 1 meV, γ = 0.1 meV, ∆ = 5 meV and ωa = 1000 meV. We conclude that the method based on the GFT is in perfect agreement with the QRT approach and reproduces 995 997 999 1001 1003 1005 ω(meV ) 0 0.15 0.30 0.45 S (ω ) a Figure 4: Emission spectrum of the cavity based on the GFT as a solid blue line and the corresponding numerical calculation based on the QRT approach as a dashed red line. The parameters values are g = 1 meV, γ = 0.005 meV, κ = 2 meV, P = 0.005 meV, ∆ = 0 meV, ωa = 1000 meV. very well the emission spectrum of the QD-cavity system. For comparison purposes with our GFT ap- proach, we have also implemented the numerical method based on the QRT for the QD-cavity system (see details in section iv). In the con- ducted simulations, we have considered the same truncation level in the bare-state basis, e.g., Nexc = 10. Moreover, we have solved numerically the dynamical equations of the system given by the Eq. (34) until time tmax = 2 17 ps for obtaining an acceptable resolution in frequency domain, it is ∆ω ≈ 0.048 meV. In order to test the performance of the GFT approach in terms of efficiency, we have compared the computational time involved on the numerical calculation of the emission spectrum of the cavity at four different excitation numbers Nexc. Table 1 shows in first column the excitation number. Second and third columns show the elapsed time (CPU time) in seconds during the simulations for the GFT and the QRT approach, respectively. It is worth mentioning that we have considered, for this comparison, exactly the same 080008-10 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. Table 1: Comparison of computational time (CPU time) between the Green’s Functions Technique (GFT) and the Quantum Regression Theorem (QRT) approach for the numerical calculation of the emission spectrum of the cavity. The computa- tional times were made using a commercial Intel(R) Core(TM) i7 − 4770 processor of 3.4 GHz ×8, and 12 GB RAM. Excitation CPU time(s) CPU time(s) number (Nexc) for the GFT for the QRT 5 0.4 92.5 10 2.0 273.4 20 14.0 390.2 40 100.2 673.8 resolution in the frequency domain and the numer- ical calculations were carried out with the same parameters values as in Fig. 4 for both methods. It is straightforward to observe that the QRT ap- proach is time-consuming compared with the GFT approach when the excitation number is increased. From the computational point of view, it is due to the fact that the QRT approach requires solving a large number of coupled differential equations in contrast to the GFT approach which requires a relatively small system of algebraic equations. V. Conclusions We have presented the GFT as an alternative methodology to the QRT approach for calculating the two-time correlation functions in open quantum systems. We have applied the GFT and the QRT approach for calculating the emission spectrum in a QD-cavity system. In particular, the performance of the GFT in terms of accuracy and efficiency by comparison of the emission spectrum of the cavity and the quantum dot is demonstrated, as well as by comparison of the computational times involved during the numerical simulations. In fact, we have shown that the GFT offers a computational advan- tage, namely, the speeding up numerical calcula- tions. We conclude that the GFT allows to over- come the inherent theoretical difficulties presented 995 997 999 1001 1003 1005 ω(meV ) 0 4 8 12 16 S (ω ) a Figure 5: Emission spectrum of the quantum dot based on the GFT as a solid blue line and the corre- sponding numerical calculation based on the QRT approach as a dashed red line. The parameters val- ues are κ = 5 meV, P = 1 meV, g = 1 meV, γ = 0.1 meV, ∆ = 5 meV and ωa = 1000 meV in the QRT approach, i.e., to find a closure condi- tion on the set of operators involved in the dynam- ical equations. We mention that our methodology based on the GFT can be extended for calculat- ing the emission spectrum in significant situations where the quantum dots are in biexcitonic regime or when the quantum dots are coupled to photonic cavities. Acknowledgements - EAG acknowledges the fi- nacial support from Vicerrectoŕıa de Investiga- ciones at Universidad del Quind́ıo through research grant No. 752. HVP acknowledges the financial support from Colciencias, within the project with code 11017249692, and HERMES code 31361. Appendix The emission spectrum in a QD-cavity system can easily be computed by taking into account that the dynamics of the operators Ĝ(τ) and ρ̂S(τ) are gov- 080008-11 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. erned by the same Linblad master equation, i.e., dĜ(τ)/dτ = LĜ(τ) with L the Liouvillian super- operator discussed in section II. Moreover, it has effectively a larger tensor rank than the reduced density operator of the system. Thus, we can write the dynamical equations for the Green’s functions operator in a component form dGα̃(τ) dτ = ∑ β̃ Lα̃β̃Gβ̃(τ), (36) together with the initial condition Gβ̃(0). Here, the symbol α̃ corresponds to a composite index for la- beling the states of the reduced density operator of the system, e.g., for indexing both matter and photon states in the QD-cavity system, see section III. for details. It is worth mentioning that Gβ̃ and Lα̃β̃ act as a column vector and a matrix in this notation. In order to obtain the solution to the Eq. (36) in frequency domain, we perform a Laplace Transform and it explicitly takes the form −G̃α̃(0) = ∑ β̃(Lα̃β̃ − iωδα̃β̃)G̃β̃(iω). It is straight- forward to obtain the solution by performing the matrix inversion to Mα̃β̃ = (iωδα̃β̃ −Lα̃β̃) and the emission spectrum is computed easily in terms of the initial conditions given by G̃β̃(iω) = ∑ α̃ M−1 β̃α̃ G̃α̃(0). (37) The initial conditions are obtained by evaluating the Green’s function operator at τ = 0. [1] H Walther et al., Cavity quantum electrody- namics, Rep. Prog. Phys. 69, 395603 (2006). [2] A Kavokin, J J Baumberg, G Malpuech, F P Laussy, Microcavities, Oxford University Press (2007). [3] H Altug, D Englund, J Vuckovic, Ultrafast photonic crystal nanocavity laser, Nat. Phys. 2, 484 (2006). [4] Y Mu, C M Savage, One-atom lasers, Phys. Rev. A 46, 5944 (1992). [5] R M Stevenson, R J Young, P Atkinson, K Cooper, D A Ritchie, A J Shields, A semi- conductor source of triggered entangled photon pairs, Nature 439, 179 (2006). [6] T M Stace, G J Milburn, C H W Barnes, Entangled two-photon source using biexciton emission of an asymmetric quantum dot in a cavity, Phys. Rev. B 67, 085317 (2003). [7] N Gisin, G Ribordy, W Tittel, H Zbinden, Quantum cryptography, Rev. Mod. Phys. 74, 145 (2002). [8] C Monroe, Quantum information processing with atoms and photons, Nature 416, 238 (2002). [9] Y Todorov, I Sagnes, I Abram, C Minot, Purcell enhancement of spontaneous emission from quantum cascades inside mirror-grating metal cavities at THz frequencies, Phys. Rev. Lett. 99, 223603 (2007). [10] J Wiersig et al., Direct observation of cor- relations between individual photon emission events of a microcavity laser, Nature 460, 245 (2009). [11] G Khitrova et al., Vacuum Rabi splitting in semiconductors, Nat. Phys. 2, 81 (2006). [12] J P Reithmaier et al., Strong coupling in a single quantum dot-semiconductor microcavity system, Nature 432, 197 (2004). [13] L Mandel, E Wolf, Optical Coherence and Quantum Optics, Cambridge University Press, (1997). [14] O Jedrkiewicz, R Loudon, Atomic dynamics in microcavities: absorption spectra by Green function method, J. Opt. B. - Quantum S. O. 2, R47 (2000). [15] N V Hieu, N B Ha, Time-resolved lumi- nescence of the coupled quantum dotmicrocav- ity system: general theory, Adv. Nat. Sci.: Nanosci. Nanotechnol. 1, 045001 (2010). [16] D F Walls, G J Milburn, Quantum Optics, Springer-Verlag, Berlin (1994). [17] M Lax, Quantum Noise. IV. Quantum Theory of Noise Sources, Phys. Rev. 145, 110 (1966). [18] S Swain, Master equation derivation of quan- tum regression theorem, J. Phys. A. - Math. Gen. 14, 2577 (1981). 080008-12 Papers in Physics, vol. 8, art. 080008 (2016) / E. A. Gómez et al. [19] E del Valle, F P Laussy, C Tejedor, Lumines- cence spectra of quantum dots in microcavities. II. Fermions, Phys. Rev. B 79, 235326 (2009). [20] N Quesada, H Vinck-Posada, B A Rodŕıguez, Density operator of a system pumped with po- laritons: a Jaynes-Cummings-like approach, J. Phys. Condens. Matter 23, 025301 (2011). [21] C A Vera, N Quesada, H Vinck-Posada, B A Rodŕıguez, Characterization of dynamical regimes and entanglement sudden death in a microcavity quantum dot system, J. Phys. Con- dens. Matter 21, 395603 (2009). [22] N Ishida, T Byrnes, F Nori, Y Yamamoto, Photoluminescence of a microcavity quantum dot system in the quantum strong-coupling regimes, Sci. Rep. 3, 1180 (2013). [23] T Quang, G S Agarwal, J Bergou, M O Scully, H Walther, K Vogel, W P Schleich, Calcula- tion of the micromaser spectrum. I. Greens- function approach and approximate analytical techniques, Phys. Rev. A 48, 803 (1993). [24] E T Jaynes, F W Cummings, Comparison of quantum and semiclassical radiation theo- ries with application to the beam maser, Proc. IEEE 51, 89 (1963). [25] M O Scully, M S Zubairy, Quantum Op- tics, Cambridge, Cambridge University Press (1996). [26] J I Perea, D Porras, C Tejedor, Dynamics of the excitations of a quantum dot in a micro- cavity, Phys. Rev. B 70, 115304 (2004). [27] H P Breuer, F Petruccione, The theory of open quantum systems, Oxford, University Press (2002). [28] A Rivas, S F Huelga, Open Quantum Systems, Springer (2012). 080008-13