Green's functions technique for calculating the emission spectrum in a quantum dot-cavity system

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. In particular, we investigate the potential of this theoretical approach by its application to compute the emission spectrum of a dissipative system composed by a single quantum dot inside of a semiconductor cavity. We also describe a simple algorithm based on the Green's functions technique for calculating the emission spectrum of the quantum dot as well as of the cavity which can easily be implemented in any numerical linear algebra package. We find that the Green's functions technique demonstrates a better accuracy and efficiency in the calculation of the emission spectrum and it allows to overcome the inherent theoretical difficulties associated to the direct application of the quantum regression theorem approach.


Introduction
The measurement and control of light produced by quantum systems have been the focus of interest of cavity quantum electrodynamics [1,2]. Specially, the emission of light powered by solid-state devices coupled to nanocavities is an extensive area of research due to its promising technological applications, such as infrared and low-threshold lasers [3,4], single and entangled photon sources [5,6], as well as various applications in quantum cryptography [7], and quantum information [8]. Experiments with semiconductor quantum dots (QDs) 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 enhancement [9], photon anti-bunching [10], vacuum Rabi splitting [11] and strong light matter coupling [12]. These and many others quantum phenomena are being confirmed experimentally by observing the power spectral density of the light (PSD) emitted by quantumdot-cavity systems (QD-Cavity). Thus, the PSD or so-called emission spectrum is the only relevant information about the system which allows to study the properties of light via measurements on correlations functions as stated by the Wiener-Khintchine theorem [13]. In order to compute the absorption or emission spectrum in open quantum systems, more precisely, in QD-Cavity systems different approaches have been developed from theoretical point of view. For example, the method of the thermodynamic Green functions which is applied to the determination of the susceptibilities and absorption spectra of atomic systems embbeded in nanocavities [14], and the time-resolved photoluminescence approach whose application allows to determine the emission spectrum by consideration of an additional subsystem called the photon reservoir [15]. However, these methods have their own approximations and restrictions and therefore are not widely used. Frequently, the emission spectrum in QD-Cavity systems is computed through the Quantum Regression Theorem (QRT) [16,17,18], since it relates the evolution of mean values of observables and the two-time correlation functions. It is worth mentioning that this approach can be difficult to implement in a computer program, it due to that computational complexity of QRT approach increases significantly as the number of QDs or modes inside the cavity, and the dimensionalities of the Hilbert spaces are large. In general, this approach is time-consuming due to that it requires to solve a large system of coupled differential equations, and numerical instabilities can arise. Moreover, theoretical complications can appears related to dynamics of the operators involved, as we will point out in the next section. In spite of this, the QRT approach is widely used for theoretical works, for example, in studies of the luminescence spectra of coupled light-matter systems in microcavities in the presence of a continuous and incoherent pumping [19,20], and 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 calculation of the micromaser spectrum [23], as a methodology in which the two-time correlation function is treated as a Green's function that decays as the off-diagonal elements of the reduced density matrix of the system for a very specific initial conditions. Nevertheless, this approach has not been widely noticed in many significant situations in open quantum systems. Possibly, it is due to their work having a limitation of implementation. The purpose of this work is to present a simple, but efficient numerical method based on QRT formalism which overcome the inherent difficulties associated to the direct application of the QRT, by solving the dynamics of the system in the frequency domain directly. This paper is organized as follows: Section 2 review the theoretical background of quantum regression theorem and its relationship with the Green's functions technique. Section 3 deal with a concrete application of our proposed method for computing the emission spectra of a dissipative QD-Cavity system. In section 4 we show the numerical calculations of the emission of spectrum for the cavity and the quantum dot from both GFT and QRT approaches. Finally, we conclude in the last section.

Theoretical background
One of the most important measurements when the light excites resonantly a QD-Cavity system is the emission spectrum of the system. From theoretical point of view, it is assumed that corresponds to a stationary and ergodic process which can be calculated as a PSD of light using the well-known Wiener-Khintchine theorem [13]. It states that the emission spectrum is given by the Fourier Transform of the correlation function (two-time expectation value) of the operator fieldĉ, where K(τ) = ĉ † (t + τ)ĉ(t) and the normalizing factor is the population n c at the steady-state. In order to calculate the two-time expectation value is frequently used the QRT which states that if a set of operators {Ô j (t + τ)} satisfy the dynamical equations d dτ is valid for any operatorÔ(t) at arbitrary time t. It is worth mentioning that vality of this theorem holds whenever a closed set of operators are involved in the dynamics. In general, to obtain the closed set of operators can be difficult or an impossible task, since it must be added as many operators as necessary in order to close the dynamics of the system. For example, in order to compute the emission of spectrum in a simple model of QD-Cavity system [20,21] two new operadors are required due to that the field operators in the interaction picture does not lead to a complete set. Before we consider the Green's functions technique, we will briefly describe the calculation of the QRT in an alternate form which will be the starting point in the following section. Lets consider a system operadorÂ which does not operate on the reservoir, then its singletime expectation value in the Heisenberg picture is given by The operatorρ S ⊗R (t) =ρ S (t) ⊗ρ R (t) depics the composite density operador of the system and reservoir.
It is worth pointing out that the dynamics of the system depends directly onρ S ⊗R (t) for all times, but the validity of the Markovian approximation requires that the state of the system is sufficiently well described bŷ ρ S (t) = T r R (ρ S ⊗R (t)), therefore it is sufficient to writê In what follows, we change to the Schrödinger representation usingÂ(t + τ) = U † (t +τ, t)Â(t)Û(t +τ, t) withÛ(t +τ, t) being the unitary time-evolution operator, and after tracing over degrees of freedom of the reservoirs, we have where the reduced density operador for the system is given byρ Then, if theρ S (t + τ) satisfies the Markovian master equation dρ S (t + τ)/dτ = Lρ S (t + τ) with L the Liouvillian superoperator, the evolution of Â (t + τ) can be computed by solving the dynamics of the master equation. To calculate the two-time correlation function Â (t + τ)B(t) whereÂ(t + τ) andB(t) are arbitrary Heisenberg operators, we proceed in a similar manner, it is, where we have used the well-known composition and inversion properties of the evolution operator. Then, the two-time operator is given bŷ By comparison of the Eq. (3) and Eq. (4), we find that G(t + τ) is an operator that obeys the same dynamical equations asρ S (t + τ), but as function of τ. It is, dĜ(t + τ)/dτ = LĜ(t + τ) with the boundary conditionĜ(t) = B(t)ρ S (t) at arbitrary time t. Hence, in the long-time limit the QRT reads, S ⊗RÛ † (τ)] is the Green's functions operator, and the operatorsÂ,B andρ (ss) S ⊗R are written in the Schrödinger representation. The superscript "(ss)" refers to the steady state of the reduced density operator of the system. After taking the Laplace transform on Eq. (6), we obtain an expression for the emission of spectrum in terms of the Green's functions operator, it is, Prior to leaving this section, we mention that this equation will be considered for computing the emission spectrum due to the cavity as well as the quantum dot, e.g. by considering the photon and fermionic operators in a separated way. Therefore, we will describe in the next subsetion a general approach that can be applied for both cases.

Algorithm for the Green's functions technique
Before to describe a simple algorithm for calculating the emission spectrum, we take into account that the dynamics for both opertorsĜ(τ) andρ S (τ) are governed by the same Master equation, i.e., dĜ(τ)/dτ = LĜ(τ) with L the Liouvillian superoperator, that efectivelly has a larger tensor rank than the reduced density operator of the system. So, we can write the dynamical equations for the Green's functions operator in a component form: together with the initial condition Gβ(0). The symbol α is a composite index for labeling 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 3 for example. Hence, Gβ and Lαβ acts as a column vector and a matrix in this notation. (ii) To obtain the solution to the Eq.
We perform the invertion of the matrix Mαβ = (iωδαβ − Lαβ) and finally, the spectrum of emission is computed in terms of the initial conditions given by,Gβ These initial conditions are easily obtained by evaluating the Green's function operator at τ = 0.
3. Application to the quantum dot-cavity system

Model
In order to apply our proposed method for calculating the emission spectrum in QD-Cavity system, we will consider a simple but illustrative system composed of a quantum dot interacting with a confined mode of the electromagnetic field inside a semiconductor cavity. This quantum system is well described by the Jaynes-Cummings Hamiltonian [24] where the quantum dot is described as a fermionic system with only two possible states, |G and |X are the ground and excited state, respectively.σ = |G X| and a (σ † = |X G| andâ † ) are the annihilation (creation) operators for the fermionic system and the cavity mode. g is the light-matter coupling constant, and we have set = 1. We also define the detuning between frequencies of the quantum dot and the cavity mode as ∆ = ω X − ω a , moreover ω X is the energy to create an exciton and ω a is the energy associated to the photons inside de cavity, respectively. This Hamiltonian system is far away for 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 include the effects of environment on the dynamics of the system, we consider the usual approach to model an open quantum system by considering a whole system-reservoir hamiltonian which is frequently splitted in three parts. The first part corresponds to the system of quantum dot-microcavity. The second part is the hamiltonian of the reservoirs and finally, the third part which is a bilinear coupling between the system and the reservoirs [26]. After tracing out the degrees of freedom of all the reservoirs and assuming the validity of the Born-Markov approximation, one arrives to a master equation for the reduced density matrix of the system, Where γ is the decay rate due to the spontaneous 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. Fig. 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 process begin when the light from the pumping laser enters into the cavity and excites one of the quantum dots in the QD layer. Thus, light from this source couples to the cavity and a fraction of photons escapes through the partly transparent mirror from the cavity and goes to the spectrometer for measurements of the emission of spectrum. A general approach for solving the dynamics of the coupled system, consist of writting the Bloch equations for the reduced density matrix of the system 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 matrix elements as ρ S αn,βm = αn|ρ S |βm . Hence, the Eq. (11) explicitly reads, − 2δ αG δ βG ρ S Xn,Xm + δ βX ρ S αn,Xm + P 2 2δ αX δ βX ρ S Gn,Gm − δ αG ρ S Gn,βm − δ βG ρ S αn,Gm .
Note that we use the convention that all indices written in greek alphabet are used for matter states and take values |G , |X , and the indices written in latin alphabet are used for Fock states and take values 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 opera-torN =â †â +σ †σ . The closed system and the number of excitations of the system is conserved, i.e., [Ĥ S ,N] = 0. It allows us to organize the states of the system through the number of excitations criterion such that the density matrix elements ρ Gn,Gn , ρ Xn−1,Xn−1 ρ Gn,Xn−1 and ρ Xn−1,Gn are related by having the same number of quanta. It is, subspaces 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 (N exc ).

Emission spectrum of the cavity
In order to compute the emission spectrum of the cavity, we will consider the two-time correlation function according to the Eq. (6) for the photon operator as follows: After performing the partial trace over the degrees of freedom of the system we have that 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 between the system and the reservoir must be unimportant even at the steady state. Thus, the density operator system-reservoir can written asρ (ss) S ⊗R =ρ (ss) S ⊗ρ (ss) R which implies that βm + 1|ρ (ss) S ⊗R |γn =ρ (ss) Replacing the previous expression in Eq. (14), it is straightforward to shows that the two-time correlation function reads where the Green's functions operatorĜ(τ) is given bŷ As we pointed out in section 2, this operator must obey the same master equation as the reduced density operator of the system. In fact, the terms that only contribute in the Eq. (16) are given by the matrix elements G βm,γn (τ) ≡ βm|Ĝ(τ)|γn of the Green's functions operator. This is due to the fact that the projection operator |βm γn| enter intoĜ(τ) in the same way as into the reduced density operator of the system.
By replacing the Eq. (18) into Eq. (17) we find that the Green's functions operator explicitly readŝ Note that from this expression is easy to identify the nonzero matrix elements of the Green's functions operator that contribute to the emission spectrum. Finally, after performing the Laplace transform we have that the emission spectrum of the cavity is given by It is worth mentioning that the initial conditions may be obtained by evaluating the Green's function operator at τ = 0, then using the fact that the time evolution operators become the identity and T r R [ρ (ss) R ] = 1, we obtain a set of initial conditions given bỹ Note that this set of initial conditions corresponds to the asymptotic solution of the Bloch equations for the reduced density matrix of the system.

Emission spectrum of the quantum dot
In order to compute the emission spectrum of the quantum dot, we will consider the two-time correlation function given by Eq. (6), but for the case of the matter operator: It is straightforward to show after performing the partial trace over the degrees of freedom of the system that the two-time correlation function reads where the Green's functions operatorĜ(τ) is given bŷ Assuming again the validity of the Markovian approximation and taking into account the number of excitations criterion, we have that the density operator systemreservoir 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 .
By inserting the Eq. (24) into Eq. (23) we find that the Green's functions operator explicitly readŝ Analogously as in section 3.2, we identify the nonzero matrix elements of the Green's functions operator that contribute to the emission spectrum and after performing the Laplace transform the emission spectrum of the quantum dot is given by where n σ = σ †σ is the normalizing factor at the steady-state. Taking into account that the initial conditions are obtained by evaluating the Green's function operator at τ = 0, we have the time evolution operators become the identity and T r R [ρ (ss) R ] = 1, thus, we obtain a set of initial conditions given bỹ

Results and Discussion
In this section, we compare the numerical calculations based on GFT and QRT approach for the emission spectrum of the cavity as well as the quantum dot. Due to that the QD-Cavity system can display two different dynamical regimes by changing the values of the free parameters of the system and transitions between these two regimes can be achieved when the loss and pump rates are modified. Particularly, in the strong coupling regime the relation P/κ g holds and the relation P/κ g remains valid for the weak coupling regime. Fig. 3 shows the numerical calculations of the emission spectrum due to the cavity in the weak coupling regime, the parameters values are g = 1 meV, γ = 0.005 meV, κ = 0.2 meV, P = 0.3 meV, ∆ = 2 meV, ω a = 1000 meV. Panel (a) shows the emission spectrum for the GFT compared to the QRT approach. Panel (b) shows the quantity |S (ω) GFT − S (ω) QRT | as a measurement of the error between the numerical calculations in the emission spectrum. For this set of parameters values, we can easily to identify two peaks associated to the modes of the cavity and the quantum dot, it is ω a ≈ 998.3meV and ω X ≈ 1000.3meV, respectively. Fig. 4 shows the same calculations as in Fig. 3, but in the strong coupling regime and the parameters values are g = 1 meV, γ = 0.005 meV, κ = 2 meV, P = 0.005 meV, ∆ = 0.0 meV, ω a = 1000 meV. In the case of resonance, the modes associated to the cavity and the quantum dot do not match, but repel each other, resulting in a structure of two separate peaks a distance 2g ≈ 2meV. Fig. 5 shows the numerical calculations of the emission of spectrum due to the quantum dot with a high value of the rate κ = 5 meV and a smaller, although non negligible, pumping P = 1 meV. The rest of parameters values are g = 1 meV, γ = 0.1 meV, ∆ = 5 meV, ω a = 1000 meV. We observed that our numerical method based on GFT is in full agreement with the QRT approach and reproduces very well the spectrum of emission associated with this system. The quantity |S (ω) GFT − S (ω) QRT | shows the discrepancy between both methods which is the order of 10 −3 − 10 −2 as it is seen in Fig. 4 and Fig. 3. Mainly, the discrepancy between both methods is due to the numerical errors accumulated in the numerical integration of the Bloch equations in the QRT approach, it causes some differences in the spectrum with respect to the results computed by the GFT. Note that, there is no integration of any equations in the GFT, therefore, we expect a more accurate spectrum of emission. We mention that for the numerical calculations based on QRT approach, we have followed the Ref. [21]. In order to test the performance of the numerical method, we regard four calculation times for computing the emission spectrum of the cavity based on GFT and QRT approach at different excitation numbers. Table 1 shows in first column the excitation number, i.e. the truncation level in the bare-state basis for the numerical calculations involved. Second and third column show the results of elapsed time in seconds for both the GFT and QRT approach, respectively. Note that for comparison purposes all numerical calculations were performed at the same truncation level i.e. N exc = 10. Additionally, we have solved numerically the Bloch equations (see Eq. (12)) until time t max = 2 17 in order to obtain a good resolution in the frequency domain for the QRT approach, i.e. ∆ω ≈ 0.048. Hence, we have evaluated the emission spectrum for the GFT in a grid with the same resolution in the frequency domain (we emphasize that QRT approach is time consuming due to the number of coupled differential equations to be solved, rather than the number of evaluations in the grid size used). In addition, the numerical calculations were carried out with the same parameters values as in Fig. 4 for both GFT and QRT approach. We found that our numerical approach based on GFT is very efficient and accurate for calculating the emission spectrum in QD-Cavity systems. Moreover, this method can easily be implemented in the numerical linear algebra packages as well as in any programming language.

Conclusions
We have developed the Green's function technique as an alternative methodology to the QRT for calculating the two-time correlation functions in open quantum systems. In particular, we have shown the performance of the Green's function technique by calculating the emission spectrum in an open quantum system composed by a quantum dot embedded in a microcavity. This theoretical approach is rather general and allows to overcome the inherent theoretical difficulties presented in the direct application of the QRT, i.e., to find a closure condition on the set of operators involved in the dynamics equations, by considering that all coherences asymptotically vanish, and remains only the reduced density matrix elements which are ruled by the number of excitations criterion. We have shown that the Green's function technique offers several computational advantages, namely, the speeding up numerical computations via a transformation of the dynamics of the master equation in a set of linear algebraic equations, which are efficiently solvable by a numerical linear algebra routine, a faster convergence and significant reduction of compu- tational time since the emission spectrum is calculated as a sum of terms of non-diagonal matrix elements of the reduced density operator of the system. We mention that our methodology can be extended for calculating the emission spectrum in significant situations of quantum dots in biexcitonic regime or involving coupled photonic cavities.