Theory of non-Markovian dynamics in resonance fluorescence spectrum

TL;DR

This paper provides a comprehensive theoretical analysis of non-Markovian effects on the resonance fluorescence spectrum of a driven quantum dot, highlighting how phonon interactions influence spectral features under various conditions.

Contribution

It introduces an exact solution for dot-phonon coupling in a non-Markovian regime and compares it with Markovian approximations, revealing the impact of phonon bath temperature and detuning.

Findings

Phonon coupling effects vanish at resonant pump.

Sideband narrowing and asymmetry increase with phonon temperature.

Exact and Born approximation results coincide for dot-phonon interactions.

Abstract

We present a detailed theoretical study of non-Markovian dynamics in the fluorescence spectrum of a driven semiconductor quantum dot (QD), embedded in a cavity and coupled to a three-dimensional (3D) acoustic phonon reservoir. In particular, we investigate the effect of pure dephasing on one of the side-peaks of the Mollow-triplet spectrum, expressed in terms of the off-diagonal element of the reduced system operator. The QD is modeled as a two-level system with an excited state representing a single exciton, and ground state represents the absence of an exciton. Coupling to the radiative modes of the cavity is treated within usual Born-Markov approximation, whereas dot-phonon coupling is discussed within non-Markovian regime beyond Born approximation. Using an equation-of-motion technique, the dot-phonon coupling is solved exactly and found that the exact result coincides with that of obtained within Born approximation. Furthermore, a Markov approximation is carried out with respect to the phonon interaction and compared with the non- Markovian lineshape for different values of the phonon bath temperature. We have found that coupling to the phonons vanishes for a resonant pump laser. For a non-resonant pump, we have characterzied the effect of dot-laser detuning and temperature of the phonon bath on the lineshape. The sideband undergoes a distinct narrowing and aquires an asymmetric shape with increasing phonon bath temperature. We have explained this behavior using a dressed-state picture of the QD levels.

Equations

$$H_0$$

$$H_R$$

$$H_{SR}$$

$$H_D=\frac{{\omega }_0}{2}\left(|\alpha ⟩⟨\alpha |-|\beta ⟩⟨\beta |\right),$$

$$H_{DR}=\sum _k{g}_k^D\left(|\alpha ⟩⟨\beta |a_k+|\beta ⟩⟨\alpha |a_k^{†}\right),$$

$$H=H_0+H_D+H_{DR}.$$

$$S({\omega }_0)=|{\wp }_{\alpha \beta }{|}^2\mathrm{Re}{\int }_0^{\infty }d\tau ⟨E^{(-)}(0)E^{(+)}(\tau )⟩e^{i{\omega }_0\tau },$$

$$E^{(+)}(t)=\sum _k{\epsilon }_k{a}_k(t),$$

$$S({\omega }_0)={\bar{I}}^2\mathrm{Re}{\int }_0^{\infty }d\tau ⟨{\sigma }_{ab}(0){\sigma }_{ba}(\tau )⟩e^{i{\omega }_0\tau },$$

$$S({\omega }_0)={\bar{I}}^2\mathrm{Re}{\int }_0^{\infty }d\tau Tr\{{\sigma }_{ba}\Omega (t)\}{e}^{i{\omega }_{0}t},$$

$$\Omega (t)=e^{-i{H}_{0}t}\bar{\rho }{\sigma }_{ab}{e}^{i{H}_{0}t},$$

$$H(t)=$$

$$+\sum _k{\omega }_k{a}_k^{†}{a}_k+\sum _k{g}_k({\sigma }_{ab}+{\sigma }_{ba})(a_k+a_k^{†})$$

$$+\sum _k{\omega }_q{b}_q^{†}{b}_q+\sum _q{\lambda }_q{\sigma }_{aa}(b_q+b_q^{†}),$$

$$\tilde{H}=$$

$$+\sum _k{g}_k({\sigma }_{ab}{a}_k+{\sigma }_{ba}{a}_k^{†})$$

$$+\sum _k{\omega }_q{b}_q^{†}{b}_q+\sum _k{\lambda }_q{\sigma }_{aa}(b_q+b_q^{†}),$$

$${\sigma }_z(t)$$

$${\sigma }_{ab}(t)$$

$$|+⟩$$

$$|-⟩$$

$$\bar{H}\simeq H_0+H_V,$$

$$H_0$$

$$H_S$$

$$H_R$$

$$H_P$$

$$H_V$$

$$H_{dR}$$

$$H_{dP}$$

$$H_{SR}$$

$$S({\Delta }_0)=$$

$$\displaystyle-{\bf s}^{2}\Omega_{-+}(s)\bigg{]}_{s=-i\Delta_{0}},$$

$$S({\Delta }_0)\simeq S_{+}({\Delta }_0)={\bar{I}}^2{\mathbf{c}}^2\mathrm{Re}[{\Omega }_{+-}(s=-i{\Delta }_0)].$$

$$\rho (t_0)={\rho }_R(t_0)\otimes {\rho }_P(t_0)\otimes {\rho }_S(t_0),$$

$${\rho }_R(t_0)=\prod _k|0_k⟩⟨0_k|,$$

$${\rho }_P(t_0)=\frac{\exp (-H_P/k_{B}T)}{\mathrm{Tr}[\exp (-H_P/k_{B}T)]}.$$

$$\Omega (t)=e^{-i\bar{H}t}\bar{\rho }{\sigma }_{ab}{e}^{i\bar{H}t}$$

$$\Omega (0)=\bar{\rho }{\sigma }_{ab},$$

Full text

††preprint: APS/123-QED

I Introduction

The laws of quantum mechanics allow for quantum computersvon Neumann (1955); Benioff (1980); Feynman (1982), which are known to be significantly more powerful than classical computers. In a quantum computer, information is stored in quantum bits (qubits), rather than classical bitsDivincenzo (1997); Schumacher (1995). A single qubit represents a zero or one and is a two-level system with two energy levels which can be used to store and process the informationLoss and DiVincenzo (1998). An example of a two-level system, frequently used in quantum optics, is composed of the ground and excited states of an atom. Semiconductor QDs can be modeled as a two-level system with one exciton in the excited stateZrenner et al. (2002); Brunner et al. (2009); Strauß et al. (2016). Semiconductor QDs embedded inside a cavity has been a subject of intense research as a promising candidate for quantum computation and information tasks, as well as a source of single photonKuhlmann et al. (2015).

Recently, experiments on the cavity embedded QDs have been reported to show different spectral features of the Mollow-triplet fluorescence spectrumUlhaq et al. (2013); Cui and Raymer (2006); McCutcheon and Nazir (2010). In particular, a dot coupled to the acoustic phonon bath in the super-ohmic regime has shown modified spectral features as a function of the phonon bath temperature. More precisely, the triplet sideband is observed to show a systematic spectral sideband broadening for both resonant and off-resonant cases. This problem was studies experimentallyUlrich et al. (2011) and anlyzed theoreticallyRoy and Hughes (2012) in terms of the usual Born-Markov approximation. However, the pure dephasing process a 3D phonon bath is in the form a super-ohmic independent boson model (IBM), which is known to be highly non-MarkovianMcCutcheon and Nazir (2013); Weiler et al. (2012); Krummheuer et al. (2002); Vagov et al. (2014), and these results have been studied and discussed in terms of usual Born-Markov approximationHughes and Agarwal (2017); Kryuchkyan et al. (2017); Roy and Hughes (2012).

The system correlation function, for non-Markovian interactions (for e.g. nuclear spinsCoish and Loss (2004), phononsKaer et al. (2010)), decays with a typical time scale given by the correlation time $\tau_{c}$ which never dies-off to zero. In other words, the correlation time is non-zero and can be larger than the system decay time $\tau_{S}$, which is the signature of a strongly history-dependent non-Markovian interaction.

Furthermore, the equation of motion for correlation function has an additional term known as irrelevant part, which is non-zero for non-Markovian interactions. This additional term vanishes for a Markov approximation, when the well-known QRT can be applied to find the system correlationSwain (1981). Recent theories [for e.g. Ref. Roy and Hughes, 2012] discussing fluorescence spectrum in solid-state systems rely on a history-independent Markov process and apply the usual QRT, giving rise to an exponential decay of the system correlation. Often, physical processes in solid-state systemsThanopulos et al. (2017); Toyli et al. (2016); McCutcheon (2016) are highly non-Markovian (history-dependent), and so the resulting spectrum using the QRT can no longer be used to describe their spectral properties.

In this paper, we analyze the effect of pure dephasing due to a 3D acoustic phonon bath on the fluorescence spectrum of a cavity coupled to a semiconductor QD. The associated emission spectrum can be directly related to a correlation function for system observables, which we evaluate beyond the Markov approximation using a Nakajima-Zwanzig GMEFick and Sauermann (1990); Swain (1981). Assuming a large band-width, coupling to the radiation modes of the cavity can be treated within Born-Markov approximation, which is relevant to cavity quantum electrodynamics (cavity-QED) experimentsMuller et al. (2007). The dot is represented by a two-level system with an exciton in the excited state coupled to a phonon reservoir, and can be modeled with the usual IBMMahan (1990).

The resultant fluorescence spectrum has three components due to dressing of the levels by a pump laserMollow (1969); Cohen-Tannoudji and Grynberg (2004), and we project the system into a dressed state basis which allows us to characterize the three components of triplet separatelyHopfmann et al. (2017); Laucht et al. (2016). We have solved the dot-phonon coupling using an exact approach beyond Born approximation, and the exact result coincides with that of obtained within Born approximation. Phonon coupling gives rise to a frequency dependent frequent-shift and dephasing, which bring non-Lorentzian features in the fluorescence spectrum. Frequency-shift and dephasing due to phonon interaction are strongly temperature-dependent and vanish for the lower temperatures. We found that in the dressed-state basis levels of interest are coupled asymmetrically to the phonons and have vanishing dephasing and frequency-shift for a resonant dot-laser frequency. We have also observed that the sideband undergoes a distinct narrowing and becomes asymmetric with increasing temperature, which is explained using the dressed-states energy levels.

This paper is organized as follows: In Sec. II, we start with discussing the setup and establish the formula for resonance fluorescence spectrum of a general two-level system. In Sec. III, we introduce the model Hamiltonian for a driven cavity-QED two-level system interacting with a phonon bath. In Sec. IV, we discuss and derive the exact form of Nakajima-Zwanzig GME for the dynamics of the reduced density matrix and correlation function. We obtain the expressions for the lineshape functions in both Markovian and non-Markovian regimes. In Sec. V, we present our results and discuss the plots in different parametric regimes. In Sec. VI, we conclude with a discussion and summary of the results. Other technical details are discussed in the Appendixes.

II Fluorescence spectrum

We start with the model Hamiltonian of a general two-level system interacting with radiation modes of the electromagnetic field, which can be written as system ($H_{S}$), field ($H_{R}$), and interaction ($H_{SR}$) in terms of the standard Jaynes-Cummings model within a rotating-wave approximation (RWA):

[TABLE]

where $\sigma_{ab}=|a\rangle\langle b|$ and $\sigma_{ba}=|b\rangle\langle a|$ are the raising and lowering operators between excited state $|a\rangle$ and ground state $|b\rangle$ in the Hilbert space of the system and $a_{k}$, $a_{k}^{\dagger}$ are the annihilation and creation operators in the Hilbert space of a set of electromagnetic modes coupled to the system. Coupling to the radiation modes is given by a coupling constant $g_{k}$ to the mode of frequency $\omega_{k}$. For simplicity of the units, we have also used $\hbar=1$. Here, we adopt the well-known theory of gedanken spectrum analyzer given in Ref. [Scully and Zubairy, 1997], and assume that the radiation field emitted by the system is detected by a two-level atom (detector) with a transition frequency $\omega_{\alpha}-\omega_{\beta}=\omega_{0}$, which is prepared in its ground state $\ket{\beta}$ initially, see Fig. 1.

The Hamiltonian for the detector is given by

[TABLE]

and coupling between detector atom $H_{D}$ and the radiation field, included in the Hamiltonian $H$, is given by the Hamiltonian

[TABLE]

where $g_{k}^{D}$ is the coupling of detector to a field mode $k$. Therefore, the Hamiltonian for entire system including system and detector, as well as coupling to the radiation modes can be written as:

[TABLE]

According to Wiener-Khintchine theorem, fluorescence spectrum, in the stationary regime and in the interaction picture with respect to the detector, is given by Fourier transform of the correlation functionScully and Zubairy (1997)

[TABLE]

where $\wp_{\alpha\beta}$ is the detector dipole matrix element with positive-frequency part of the electric field is defined by

[TABLE]

and the negative-frequency part of the electric field is $E^{(-)}(t)=[E^{(+)}(t)]^{{\dagger}}$. The quantity $\varepsilon_{k}=\sqrt{\hbar\omega_{k}/(2\epsilon_{0}V)}$ is the electric field per photon and $V$ is the effective volume of a cubic cavity resonator. We rewrite the correlation function, in Eq. (7), in terms of the system operators using the methods described in Refs. Scully and Zubairy, 1997; Cohen-Tannoudji and Grynberg, 2004:

[TABLE]

here $\bar{I}$ is the detector response function [discussed in Appendix A]. The average $\langle\dots\rangle=Tr\{\dots\bar{\rho}\}$ , in Eq. (9), is given with respect to the stationary density matrix $\bar{\rho}$, where $\bar{\rho}=\lim\limits_{T\to\infty}\frac{1}{T}\int_{0}^{T}dt\,\rho(t)$. Using the cyclic property of trace, fluorescence spectrum can be written as

[TABLE]

here operator $\Omega(t)$ is defined as

[TABLE]

and $H_{0}=H_{S}+H_{R}+H_{SR}$. Since $\sigma_{ab}$ and $\sigma_{ba}$ are operators in the system Hilbert space and $[H_{D},H_{0}]=0$, the evolution of $\Omega(t)$ is determined by the Hamiltonian of emitting system and radiation field, $H_{0}$, in the absence of detector and can be computed without using the well-known QRTSwain (1981).

III Model

III.1 Hamiltonian

We consider a driven two-level cavity-QED system with an excited state $\ket{a}$ representing a single exciton, and a ground state $\ket{b}$ with no exciton. The QD interacts with the cavity photons and a phonon reservoir which is coupled to the excited state $\ket{a}$, as shown in Fig. 1. The Hamiltonian for the total system reads,

[TABLE]

where $\omega_{ab}$ is the transition frequency of the two-level system and $\omega$ is the frequency of laser field. The photon (phonon) modes are represented by bosonic fields with frequencies $\omega_{k}$ ($\omega_{q}$) with creation and annihilation operators $a_{k}^{\dagger}$ ($b_{q}^{\dagger}$) and $a_{k}$ ($b_{q}$), respectively. The system-photon (system-phonon) coupling strength is given by $g_{k}$ ($\lambda_{q}$), and $\Omega$ (Rabi frequency) is the coupling between the two-level system and laser field. The system operators are denoted by $\sigma_{ij}=|i\rangle\langle j|$ where $i,j\in\{a,b\}$ and $\sigma_{z}=\sigma_{aa}-\sigma_{bb}$.

The explicit time dependence in $H(t)$ [Eq. (12)] can be removed by going to a rotating frame and applying a RWA. We perform the RWA on both the driving term and the system-photon coupling term, in which we neglect the rapidly-oscillating terms and keeping only the time-independent part. The resulting RWA Hamiltonian in the rotated frame is then

[TABLE]

where $\Delta=\omega_{ab}-\omega$ and $\Delta_{k}=\omega_{k}-\omega$ are the detunings of atomic and cavity frequencies from the laser pump frequency $\omega$. In general, operators in the rest frame are transformed to the rotating frame according to the following relations:

[TABLE]

Due to presence of the intense laser field the two bare states are strongly coupled to each other and give rise to two dressed states. The dressed states can be written, in terms of the bare states, as:

[TABLE]

where ${\bf c}=\cos\theta$ and ${\bf s}=\sin\theta$ with the mixing angle $\theta$ is given by $\tan\theta=\sqrt{(\Omega_{R}-\Delta)/(\Omega_{R}+\Delta)}$, where $\Omega_{R}=\sqrt{\Omega^{2}+\Delta^{2}}$ is the dressed Rabi frequency.

Furthermore, following closely the discussing on IBM in Ref. [Mahan, 1990], we apply another transformation by using a canonical transformation $H^{\prime}=e^{B}\tilde{H}e^{-B}$, where $B=\frac{1}{2}\mathbb{1}\otimes\sum_{q}\frac{\lambda_{q}}{\omega_{q}}(b_{q}^{\dagger}-b_{q})$ is an anti-hermitian operator, to write the total Hamiltonian in terms of the free and perturbed parts as, after transforming it to the dressed-state basis

[TABLE]

where the free part is

[TABLE]

and the perturbed part is given by

[TABLE]

where $\sigma_{ij}=|i\rangle\langle j|$, $i,j\in\{+,-\}$ and $\sigma_{3}=|+\rangle\langle+|-|-\rangle\langle-|$. The polaron transformation introduces a frequency shift as: $\Omega_{R}^{\prime}=\Omega_{R}-\Delta_{P}$, where $\Delta_{P}=({\bf c}^{2}-{\bf s}^{2})\sum_{q}\lambda_{q}^{2}/\omega_{q}$ is the polaron shift. We have also performed a Schrieffer-Wolff transformation on the above Hamiltonian to get rid of the energy exchange term ($T_{1}$ lifetime process) due to phonons, see Appendix B. The separation into pure dephasing and transition terms is determined by the form of system coupling, i.e. dephasing terms contain the diagonal coupling and transition terms contain the off-diagonal couplings. Accordingly, fluorescence spectrum, in Eq. (10), can be written in Laplace domain using the dressed-state representation as a three-peak spectrum [Ref. Cohen-Tannoudji and Grynberg, 2004]:

[TABLE]

where $\Delta_{0}=\omega_{0}-\omega$ is the probe detuning, we have used the relations $\Omega_{ij}=\langle j|\Omega|i\rangle$, $i,j\in\{+,-\}$ and $\Omega_{z}=\Omega_{++}-\Omega_{--}$. The polaron transformation does not affect the fluorescence spectrum in above expression as $B$ acts on the phonon mode Hilbert space and commutes with the system operators. It can be seen that the first term in Eq. (27) gives the central peak whereas last two terms give rise to satellite peaks of the Mollow-triplet, shown in Fig. 2. We are interested in the influence of pure dephasing due to phonon interaction, which only affects the off-diagonal elements of the system operator. From here and what follows, we will study the effect of pure dephasing on the fluorescence spectrum, with a particular focus on one of the Mollow-triplet sidebands (Stokes line). Alternatively, when the width of each side peak $\sim\Gamma$, [see Eq. (63), below] is small compared to the peak separation $\sim\Omega_{R}$, we approximate the spectrum near the side peak centered at $\Delta_{0}\simeq\Omega_{R}$ by, see Fig. 2

[TABLE]

In order to compute the spectrum, we will evaluate the dynamics of matrix element $\Omega_{+-}(t)={\rm Tr}[\sigma_{-+}\Omega(t)]$ in the dressed-state representation with Hamiltonian in the polaron frame, given by Eq. (18).

III.2 Initial conditions

The radiation and phonon modes are decoupled from the system for times $t_{0}<0$ (where $t_{0}$ is a time in the distant past), and prepared independently in the states described by density matrices $\rho_{R}(t_{0})$, $\rho_{P}(t_{0})$ and $\rho_{S}(t_{0})$, respectively. The interactions (photons and phonons) are switched on at this time $t=t_{0}$, and the state of entire system is described by the full density matrix $\rho(t_{0})$:

[TABLE]

where initial density matrix for photon is described by vacuum of the cavity modes

[TABLE]

and phonon modes are described by a canonical ensemble at temperature $T$:

[TABLE]

We recall the operator $\Omega(t)$, where

[TABLE]

is analogous to the density matrix operator with a modified initial condition, given by

[TABLE]

where stationary density matrix $\bar{\rho}$ and hence $\Omega(0)$ account for conditions that accumulate between the system and interactions in the time interval $t\in[t_{0},0]$. We choose an initial condition when exciton is in excited state $\ket{a}$ given by $\rho_{S}(t_{0})=\ket{a}\bra{a}$, and evolves in the presence of pump laser. We switch on the detector at $t=0$ and subsequently calculate the dynamics of $\Omega_{+-}(t)$ for $t>0$, with initial condition given by the steady state density matrix accumulated between the time interval $[t_{0},0]$.

IV Generalized Master Equation

We are interested in the dynamics of reduced system operator, after tracing over variables of photon and phonon modes; $\Omega_{S}(t)={\rm Tr_{R}Tr_{P}}\Omega(t)$. To study the dynamics of reduced system operator $\Omega_{S}(t)$, we introduce a projection superoperator $P$, defined by its action on an operator: $P\mathcal{O}(t)=\rho_{R}(t_{0})\rho_{P}(t_{0}){\rm Tr_{R}Tr_{P}}\mathcal{O}(t)$.

Both operators, $\rho(t)$ and $\Omega(t)$, follow the same von-Neumann equation, and can be written in form of exact Nakajima-Zwanzig GMESwain (1981), using $Q\rho(t_{0})=0$

[TABLE]

where $\Sigma(t)$ is the self-energy superoperator

[TABLE]

and $L$ is the full Liouvillian superoperator, defined as $L_{\alpha}\mathcal{O}=[H_{\alpha},\mathcal{O}]$ and $\alpha=0(S,R,P),V(dR,dP,SR)$. We have used the properties of projection operator: $P^{2}=P$ and

[TABLE]

also introducing its complement $Q=\mathbb{1}-P$. We can derive an equation of motion for $\Omega(t)$ analogous to the equation for $\rho(t)$ [Eq. (34)]. However, an additional term appears because $Q\Omega(0)=Q\bar{\rho}\sigma_{ab}\neq 0$, and we have assumed that the full density matrix operator is not separable for all times i.e. $\rho(t)\neq\rho_{R}(t_{0})\otimes\rho_{P}(t_{0})\otimes\rho_{S}(t)$. The resulting motion equation for $P\Omega(t)$ is then,

[TABLE]

where $\Sigma(t)$ is defined in Eq. (35) and the last term in the above equation contains $Q\Omega(0)$, i.e. the irrelevant part of $\Omega(0)$ is non-zero. When the radiation and phonon modes are described by Eqs. (30) and (31), the projection operator $P$ follows some useful identities

[TABLE]

We apply above identities (38)-(41), and perform the partial traces on Eqs. (34) and (37) over variables of radiation and phonon modes to obtain an equation of reduced system operators,

[TABLE]

where $\Sigma_{S}(t)$ and $\Phi_{S}(t)$ are reduced self-energy and irrelevant part superoperators, respectively, and $\mathcal{O}_{S}(t)=\mathrm{Tr}_{R}\mathrm{Tr}_{P}\mathcal{O}(t)=\sum_{\alpha\beta\in\{+,-\}}\mathcal{O}_{\alpha\beta}(t)\ket{\alpha}\bra{\beta}$ is the reduced system operator.

Comparing Eqs. (42) and (43), we found that the first two terms are identical but there is an additional irrelevant part, $\Phi_{S}(t)$, is present in the equation for $\Omega_{S}(t)$. This term vanishes in a Markov approximation and usual QRT can be applied to find the system correlationSwain (1981). In addition to this, we have also assumed that the full density matrix is not separable for all times which is another reason that QRT is no longer valid in the present case. Furthermore, for a non-Markovian equation the irrelevant term is non-zero and usual QRT can not be used to compute the system correlation in this casede Vega and Alonso (2008); Lax (2000); Jin et al. (2016). The equation for off-diagonal matrix element, $\Omega_{+-}$, is coupled to both diagonal $\Omega_{++},\,\Omega_{--}$ and off-diagonal $\Omega_{-+}$ elements of the reduced system operator, and can be written in the form

[TABLE]

with a non-zero irrelevant part matrix element expressed as (see Appendix C)

[TABLE]

The self-energy superoperator can be decomposed into three terms as

[TABLE]

where the first and second terms in above expression are pure dephasing ($T_{2}^{*}$ pure dephasing time) processes due to photon and phonon couplings, respectively, and the third term gives rise to transition ($T_{1}$ lifetime) due to radiative modes coupling. We treat the dot-cavity coupling self-energy within Born-Markov approximation to second-order in perturbation Liouvillain $L_{V}$, see Appendix D.1. Applying the continuum of modes for cavity given by Lorentzian density of state [Eq. (127)], we find the self-energy matrix elements in their Laplace domian, defined as $f(s)=\int_{0}^{\infty}e^{-st}f(t)dt$,

[TABLE]

[TABLE]

where $\Delta_{c}=\omega_{c}-\omega$ is the detuning of cavity from laser pump frequency and $\Gamma_{c}$ is cavity bandwidth. Here, Born approximation is justified by finding that the higher order terms in reduced self-energy are suppressed by a small parameter $\sim g^{2}/(\Omega^{\prime}_{R}\Gamma_{c})$. Similarly for phonon modes, applying a continuum of modes [see Appendix D.2] for the deformation potential coupling mechanismKrummheuer et al. (2002) $N(\epsilon)|\lambda(\epsilon)|^{2}=\alpha_{P}|\epsilon|^{3}e^{-|\epsilon|/\epsilon_{c}}$, where $\alpha_{P}$ is the phonon coupling parameter in the units of ${\emph{\rm}frequency}^{-2}$ and $\epsilon_{c}$ is the phonon cut-off frequency, we obtain

[TABLE]

here $n_{B}(\epsilon)$ is Bose function. In above expression, the dot-phonon interaction self-energy is evaluated using exact approach within non-Markovian limit and it is found that Born approximation is exact in this case, see Appendix D.3. Furthermore, the reduced self-energy matrix elements couple to the populations can as well be written within usual Born approximation

[TABLE]

[TABLE]

and similarly for the coherence $\Omega_{-+}$,

[TABLE]

Equation for the coherence in Eq. (46) contains both diagonal and off-diagonal elements of the self-energy superoperator, some of them are fast moving compare to others. In next section, we will perform a secular approximationCohen-Tannoudji and Grynberg (2004) to get rid of the fast oscillating terms.

IV.1 Secular approximation

The secular approximation consists in neglecting the fast oscillating terms in Markov equation-of-motion, and the equation for $\Omega_{+-}$ [Eq. (46)] is decoupled from populations and coherence within a secular approximationCohen-Tannoudji and Grynberg (2004). Here, we consider a general equation of motion for the operator $\Omega(t)$, without irrelevant part matrix elements:

[TABLE]

and would like to perform the secular approximation in order to get rid of the fast oscillating terms. To this end, introducing a rotating frame

[TABLE]

where $\Delta\omega$ is the total frequency shift given implicitly by the expression

[TABLE]

In the weak coupling regime, $\Omega_{R}^{\prime}\gg\Delta\omega$, the frequency shift to the leading order in $\Delta\omega$

[TABLE]

Here, the purpose of introducing a rotating frame is to get rid of all oscillating parts from $\Omega^{\prime}(t)$ and hence obtain an equation for $\Omega^{\prime}_{+-}(t)$

[TABLE]

where $\tilde{\Sigma}_{+-,+-}(t)=e^{i(\Omega_{R}^{\prime}+\Delta\omega)t}\Sigma_{+-,+-}(t)$. Due to presence of the oscillatory exponential in the last two terms in RHS of the above equation, we decompose it into slow-varying and fast-oscillating parts as:

[TABLE]

where S and F stand for slow and fast, respectively. Carrying out Markov approximation on slow term by replacing: $t^{\prime}\rightarrow t-t^{\prime}$ , $\Omega^{\prime}_{+-}(t-t^{\prime})\rightarrow\Omega^{\prime}_{+-}(t)$, and finally extending the upper limit of integration to infinity, we obtain

[TABLE]

where

[TABLE]

Condition for validity of Markov approximation: $e^{i(\Omega^{\prime}_{R}+\Delta\omega)t^{\prime}}\Sigma_{+-,+-}(t^{\prime})$ decays on a time scale $\tau_{c}\ll T_{2}$ which is the decay time of $\Omega_{+-}^{{}^{\prime}S}(t)$ and given by the relationFick and Sauermann (1990)

[TABLE]

On substituting Eqs. (49) and (50) in the above inequality and for $\Gamma_{c}\gg\Omega_{R}^{\prime},\Delta\omega$ and $\Delta_{c}$, it leads to the condition $g/\Gamma_{c}\ll 1$ and similarly for the phonon coupling $\alpha_{P}\Delta\omega^{3}e^{\Delta\omega/\epsilon_{c}}(2n_{B}(\Delta\omega)+1)/\epsilon_{c}\ll 1$. Similarly, for the fast term

[TABLE]

and the condition for validity of Markov approximation:

[TABLE]

which leads to a similar condition $g/\Gamma_{c}\ll 1$. For large $\Omega_{R}^{\prime}+\Delta\omega$ and due to presence of the highly oscillatory exponential $e^{i(\Omega_{R}^{\prime}+\Delta\omega)t}$, the effects of terms $\Omega^{\prime}_{-+}$ will eventually average out to the smaller values compared to $\Omega^{\prime}_{+-}$. Therefore, in secular approximation when

[TABLE]

we neglect the fast oscillating terms in Markovian approximation since this term will oscillate fast and average out to a smaller value compared to the slow term. On substituting for $\Sigma_{+-,-+}^{SR}(s)$ from Eq. (54) and for $\Gamma_{c}\gg\Delta_{c}$, $\Omega_{R}^{\prime}\gg\Delta\omega$, we get an explicit condition for the validity of secular approximation as: $g^{2}/(\Omega_{R}^{\prime}\Gamma_{c})\ll 1$. In the similar manner, we can also neglect the contribution form $\Omega_{++}(t)$ and $\Omega_{--}(t)$ from Eq. (46). Going back to lab frame and within secular approximation, we obtain the expression in Laplace transform:

[TABLE]

with initial condition expressed in terms of reduced self-energy matrix elements given by Eqs. (52) and (53), see Appendix E:

[TABLE]

The irrelevant part matrix elements, due to photon $G_{+-,+-}^{R}$ and phonon $G_{+-,+-}^{P}$, both are identically zero under a Markov approximation. On substituting for $\Omega_{+-}(s=-i\Delta_{0})$ from Eq. (68) in the expression (28), we obtain an expression for one-peak Markovian spectrum

[TABLE]

which is a Lorentzian line centered at $\Delta_{0}=\Omega_{R}^{\prime}+\Delta\omega$ with a width given by $\Gamma$, frequency-shift $\Delta\omega$ and decay $\Gamma$ are given by Eqs. (58) and (63), respectively, and the pre-factor is given by

[TABLE]

Here, we have substituted the expression for the self-energy matrix elements given by Eqs. (52) and (53). Expression for the lineshape in Eq. (68) is Markovian with respect to both photon and phonon interactions. However, we are interested in the non-Markovian regime with respect to phonon coupling and an equation for $\Omega_{+-}(t)$ always contains an extra small term $G_{+-,+-}^{P}(t)$ due to non-Markovian interaction, known as irrelevant part matrix element. Assuming that the irrelevant part is associated with a smallness in the present problem and in order to get the further insight, we will estimate the typical size of its contribution and find a regime where non-Markovian correction is dominant compared to its irrelevant part contribution. Equation for $\Omega_{+-}(t)$ in Laplace domain with its irrelevant part matrix element can be written as,

[TABLE]

The smallness of irrelevant part compared to non-Markovian self-energy due to phonon interaction can be justified by the inequality given by

[TABLE]

and we only keep the contribution from self-energy matrix element. Furthermore, expanding Eq. (72) in the powers of self-energy and ignoring the higher order terms, we have

[TABLE]

and assuming that irrelevant part gives rise to a small contribution and comparing it with self-energy contribution leads to the following inequality,

[TABLE]

The one-peak spectrum is centered around $\Delta_{0}\sim\Omega_{R}^{\prime}$, with a width mostly dominated by Markovian decay rate $\Gamma$, estimating the size of above inequality around $\Delta_{0}\sim\Omega_{R}^{\prime}+\Gamma$, and it is found that irrelevant part matrix element is always suppressed by a small parameter, $\Gamma/\Omega_{R}^{\prime}\ll 1$, compared to self-energy matrix element contribution, also see Appendix C. Following the above discussion, we neglect the irrelevant part from the equation for $\Omega_{+-}(t)$ and after going back to lab frame we have an equation written in Laplace domain,

[TABLE]

Above expression is Markovian with respect to photon coupling but non-Markovian in terms of phonon interaction, where Markovian frequency shift ($\Delta\omega_{R}$) and decay rate ($\Gamma_{R}$) are given by

[TABLE]

Similarly, non-Markovian frequency-dependent shift ($\Delta\omega_{P}(\Delta_{0})$) is expressed as

[TABLE]

and dephasing ($\Gamma_{P}(\Delta_{0})$) is given by

[TABLE]

On substituting for $\Omega_{+-}(s=-i\Delta_{0})$ from Eq. (76) in the expression (28), we obtain an expression for the one-peak non-Markovian spectrum

[TABLE]

where pre-factor $X$ is given by Eq. (71). It should be noted here that we have not performed Born-Markov approximation in terms of phonon interaction. In Markovian regime frequency shift $\Delta\omega_{P}(\Delta_{0})$ and dephasing $\Gamma_{P}(\Delta_{0})$ are replaced by their $\Delta_{0}=\Omega_{R}^{\prime}+\Delta\omega$ frequency parts and give rise to an exponential decay and hence to a Lorentzian line centered at $\Delta_{0}=\Omega_{R}^{\prime}+\Delta\omega$ with a width given by $\Gamma$. Whereas in the non-Markovian regime frequency-shift and dephasing due to phonon interaction are frequency dependent and lead to a non-Markovian (non-exponential) decay giving rise to non-Lorentzian features in the lineshape. We apply above obtained theoretical results to InAs/GaAs QDs, and use the cavity and phonon parameters given in Refs. M.E. and Rumyantsev, 1996; A. and Kundrotas, 1994; Ulrich et al., 2011; Weiler et al., 2012; Brunner et al., 2009; Krummheuer et al., 2002.

V Results and discussion

In this section, we plot and analyze the results obtained in the previous sections for different parameter regimes. Typical phonon parameters for GaAs are obtained from Refs. Weiler et al., 2012; M.E. and Rumyantsev, 1996; A. and Kundrotas, 1994: phonon cut-off $\omega_{c}=1\,{\rm meV}$ and coupling $\alpha_{P}=2.08\times 10^{-7}\,\mu{\rm eV}^{-2}$. For the laser and cavity, we choose following parametersRoy and Hughes (2012); Ulrich et al. (2011) $\Omega=500\,\mu{\rm eV}$, $g=50\,\mu{\rm eV}$, $\Gamma_{c}=2\,{\rm meV}$. We have varied the other parameters in the plots and explained along with the figures.

V.1 Temperature dependent frequency-shift and dephasing

In Fig. 3, we plot the frequency-shift and dephasing as a function of probe-detuning for different values of phonon bath temperatures. We observe that the effect is strongly temperature-dependent. It increases linearly in the high temperature limit and vanishes for small temperatures. In the present parameter regime for a typical GaAs QD, the dominant contribution from phonon interaction is mainly due to the frequency-shift, see Fig. 3. Here, we have assumed that the dot and cavity are resonant and set $\Delta=\Delta_{c}=500\,\mu{\rm eV}$.

V.2 Temperature dependent one-peak fluorescence spectra

In Fig. 4, we also plot the associated one-peak spectra for different phonon bath temperatures for the fixed detunings and analyze the effect of frequency-shift and dephasing on the lineshape. We observe a distinct narrowing and asymmetry in the side-peaks with increasing temperature mainly due to the frequency-shift; as frequency-shift changes with increasing temperature giving rise to non-Lorentzian features in the lineshape. This behavior is not observed in the Markovian lineshape as both frequency-shift and dephasing are constants in this case.

In the dressed-state basis, levels involved in the transitions of interest (Stokes line) are coupled asymmetrically with the phonon modes, and frequency-shift due to phonons pulls these energy levels away from the resonance bringing an additional shift. This extra shift, due to phonons, increases the level separation and reduces the number of channels for the radiative decay. In other words, some of the photons are used to compensate for this additional shift which leads in less photons coming out to the outer world or seen by the detector, eventually causing a narrowing in the sideband. Moreover, the frequency-shift appears to be an odd function with respect to probe detuning (Fig. 3) which leads to different probability of emitting and absorbing the phonons on the either side of probe detuning.

We also observe that the Lamb shift and additional broadening, due to phonon interaction, vanish for the case of resonant pump laser. This can be explained in the dressed-state representation, where both energy-levels of interest, correspond to one of the transitions (Stokes line), are coupled asymmetrically to the phonons. For the zero detuning they become degenerate leading to the vanishing frequency-shift and dephasing. We have also observed a little shifts in the peak positions due to strong dependence of frequency-shift on the probe-detuning.

VI Conclusion

We have discussed the dynamics of a driven cavity-QED system coupled to a 3d acoustic phonon reservoir with non-Markovian pure dephasing mechanism. We have observed highly modified non-Lorentzian features in the associated spectrum due to phonon interaction, and for solid-state systems it can only be described within non-Markovian regime. We have shown that quantum regression theorem is not valid in this case because of a non-zero irrelevant part, and also the full density matrix is not separable for all times. The system correlation is solved without using the quantum regression theorem beyond usual Born-Markov approximation. We have obtained the analytical formulas for both Markovian and non-Markovian one-peak lineshapes in terms of the model parameters. Both expressions look like Lorentzian lines centered around $\Delta_{0}\simeq\Omega_{R}^{\prime}$, but this analogy is not correct as shift and dephasing are also probe-dependent in case of non-Markovian lineshape. We have investigated the one-peak spectrum of the Mollow-triplet for different values of phonon bath temperature for the resonant dot and cavity. We have found vanishing shift and dephasing when laser is resonant with the dot, because levels are equally and asymmetrically coupled to phonon modes rendering them degenerate for a resonant pump laser. We have derived an exact form of Najakjima-Zwanzig generalized master equation and showed that Markov and non-Markov solutions give significantly different results. We show that non-Markovian contribution is significantly large and can be clearly seen in the spectrum. A distinct narrowing of the sideband has been reported which is contrary to the recent results leading to broadening in the presence of phonons. We have shown that the frequency-dependent shift has strong temperature dependence which causes different features like narrowing and asymmetry in the lineshape. This procedure can also be used to systematically account for features in optical spectra of a general multi-level system due to genuine non-Markovian dynamics.

Acknowledgments

AK acknowledge financial support from the Icelandic Research Fund RANNIS, CIFAR, Canada, the National Key Research and Development Program of China (Grant No. 2016YFA0301200), and the NFSC grants (No 11574025 and No. U1530401). AK thanks Bill Coish, Sigurdur I. Erlingsson, Stefano Chesi, Li-jing Jin, and Tilen Cadez for the useful discussions and feedback.

Appendix A Detector response function

The detector response function in Eq. (9)

[TABLE]

where $g_{k}^{D}=|\mathcal{\wp}_{\alpha\beta}|\varepsilon_{k}$. Applying the continuum of modes after replacing the sum via integral $\sum_{k}g_{k}^{2}\to\int_{0}^{\infty}D(\epsilon)|g(\epsilon)|^{2}d\epsilon$ and using the well-known formula

[TABLE]

($\mathcal{P}$ indicates the principal part) one can write the square of detector response function in terms of principal part and delta function as

[TABLE]

where $D_{c}(\epsilon)$ and $D(\nu)$ are photonic density of states of cavity and open space, respectively. The photonic density of states of the cavity is described by a Lorentzian density of states given by Eq. (127).

Appendix B Schreiffer-Wolff Transformation

Hamiltonian $H^{\prime}$ resulting from the polaron transformation, is

[TABLE]

Last two terms in above Hamiltonian commute with rest of the Hamiltonian and can be igonred within a secular approximation for large $\Omega_{R}^{\prime}$. We get rid of the energy exchange process due to phonon interaction using leading-order Schrieffer-Wolff transformationSchrieffer and Wolff (1966) and start from the full Hamiltonian:

[TABLE]

where individual terms are defined as

[TABLE]

We apply a transformation, $\bar{H}=e^{A}He^{-A}$, generated by an anti-hermitian operator $A=-A^{{\dagger}}$ to eliminate the transition terms to the first order. Using Baker-Campbell-Hausdorff formula and expanding $\bar{H}$ in the powers of $A$, we obtain

[TABLE]

In order to get rid of transition term $V_{2}$ to leading order, we set $V_{2}=-[A,H_{1}]$, where $A$ can be written as

[TABLE]

and $L_{1}\mathcal{O}=[H_{1},\mathcal{O}]$. Here $A$ is of order of transition term $V_{2}$. Substituting for $A$, we obtain Hamiltonian up to the second or higher order in $V_{2}$

[TABLE]

Using the definitions of $H_{1}$ and $V_{2}$, we obtain the expression for A as:

[TABLE]

Therefore, the transformed Hamiltonian to the first order in the transition terms due to phonon can be well approximated and written as free and perturbed parts as

[TABLE]

where the free part is

[TABLE]

and perturbed parts is given by

[TABLE]

which are Eqs. (18)-(26) in the main text.

Appendix C Irrelevant part matrix element

The irrelevant part is given by Eqn. (45) can be written within Born approximation after transforming to Laplace domain

[TABLE]

In particular, we want the matrix element, $\Phi_{+-,+-}(s)$ due to phonon interaction, which can be simplified and written as

[TABLE]

where irrelevant part of the stationary density matrix can be foundSwain (1981) using GME discussed in Sec. IV

[TABLE]

The propagators in Eqn. (110) has no poles at $s=0^{+}$, where $0^{+}$ is a positive infinitesimal. After performing the limit in the above expression, we substitute for $[Q\Omega(0)]_{+-}$ in the expression for $\Phi_{+-,+-}(s)$, to obtain and expression for the irrelevant term as

[TABLE]

where $\Omega_{+-}(0)=[\bar{\rho}_{S}\sigma_{ab}]_{+-}$, see Eqs. (178) and (179). Above expression can be written in terms of irrelevant part matrix element in the problem as

[TABLE]

Solving above expressions and applying the continuum of modes, we found that irrelevant part matrix element will be suppressed by $1/\Omega_{R}^{\prime}$ compared to the self-energy matrix element due to phonon interaction. In this limit, contribution from irrelevant part can be neglected compared to the contribution from the non-Markovian self-energy.

Appendix D Self-energy calculations

The reduced self-energy superoperator $\Sigma_{S}(t)$ in Eq. (44) can be transformed in to Laplace domain and using $L_{V}=L_{dR}+L_{SR}+L_{dP}$, we obtain

[TABLE]

dropping $Q$ from the exponential in Eq. (44) will not affect the final expression DiVincenzo and Loss (2005). Cross terms in above expression will vanish because $L_{dR}$ and $L_{SR}$ act on an operator in the radiation mode Hilbert space whereas $L_{dP}$ act on phonon Hilbert space and will not contribute in the final trace. Moreover, the cross terms between $L_{dR}$ and $L_{SR}$ will give rise to off-block diagonal matrix elements in self-energy matrix and will be neglected within secular approximation, see Sec. IV.1 in main text. The self-energy superoperator in Eq. (114) can be decomposed into three different parts as:

[TABLE]

where first and second terms give rise to pure dephasing ($T_{2}^{*}$ process) due to radiation and phonon modes, respectively, whereas last term in above expression leads to transition ($T_{1}$ process) due to radiation modes coupling. Free propagator in reduced self-energy expression can be expanded in the powers of interacting Liouvillian $L_{V}$ as Coish and Loss (2004)

[TABLE]

because of the form of couplings in present model only even powers $2k$ in above expression will survive in the final trace. In order to find the matrix elements of the self-energy superoperator, we write the superoperators in matrix form in dressed-state basis as

[TABLE]

where $[L_{S}]_{\alpha\beta,\gamma\delta}=Tr\{\ket{\beta}\bra{\alpha}S\ket{\gamma}\bra{\delta}\}$ and $\{\alpha,\beta\}\in\{+,-\}$. In the dressed state basis, non-interacting Liouvillian is diagonal and can be inverted to write as $2\times 2$ blocks

[TABLE]

where parallel block is

[TABLE]

and perpendicular block is given by

[TABLE]

In the similar fashion, we can find other matrices as well

[TABLE]

here we have defined new Liouvillian for commutation and anti-commutation relations: $L_{X(Y)}^{\pm}\mathcal{O}=[X_{R(P)},\mathcal{O}]_{\pm}$, where operators $X_{R}$ and $X_{P}$ are given as

[TABLE]

Furthermore, we also find the matrix for superoperator $L_{SR}$ in the dressed-state basis

[TABLE]

where we have defined the operators for left and right multiplications as:

[TABLE]

and $\mathcal{O}_{R}$ is an operator in the radiation mode Hilbert space. Reduced self-energy matrix elements of interest can be calculated according to Eq. (46) in the main text.

D.1 Self-energy for photon interaction

We apply a continuum of modes for the cavity density of states given by a Lorentzian spectrum Ma et al. (2012),

[TABLE]

where $\Delta_{c}=\omega_{c}-\omega$ is the detuning of cavity from laser pump frequency, and $\Gamma_{c}$ is cavity bandwidth. Reduced self-energy matrix elements due to radiation mode coupling giving rise to transition are given as

[TABLE]

and self-energy matrix element responsible for pure dephasing due to cavity coupling can be found as

[TABLE]

For a large band-width cavity, coupling to its radiative modes to system is treated under Markov approximation and above self-energies are replaced by their $s=-i(\Omega_{R}^{\prime}+\Delta\omega)$ frequency parts, refer main text for details.

D.2 Self-energy for phonon interaction

Similarly, we apply a continuum of modes for 3-d acoustic phononsKrummheuer et al. (2002) with an exponential cut-off at $\epsilon=\epsilon_{c}$

[TABLE]

we obtain the expression for self-energy in Laplace transform as

[TABLE]

where $\alpha_{P}=$ is the phonon coupling paramter in the units of $\mathrm{freq.}^{-2}$ and $n_{B}(\epsilon)$ is Bose function. On further simplification, above self-energy matrix element can be decomposed into real and imaginary parts after setting $s=-i\Delta_{0}$, where $\Delta_{0}$ is detuning of the probe from the pump laser frequency, as:

[TABLE]

where $\Delta\omega_{P}(\Delta_{0})={\rm Re}[\Sigma_{+-,+-}^{dP}(\Delta_{0})]$ is the frequency-shift and $\Gamma_{P}(\Delta_{0})=-{\rm Im}[\Sigma_{+-,+-}^{dP}(\Delta_{0})]$ is the dephasing due to phonon interaction.

D.3 Self-energy matrix element calculated exactly

In the previous section, we have computed the self-energy matrix element for phonon interaction to the second order in Born approximation. In this section, we will discuss an equation-of-motion method to find the phonon interaction self-energy for all orders in perturbed Liouvillain due to phonons $L_{Y}^{+}$ beyond Born approximation, and show that exact approach recovers the result obtained within Born approximation. Using the general form of superoperators matrices, the expression for self-energy matrix element due to phonons can be written in Laplace domain as:

[TABLE]

also in the time domain, we have

[TABLE]

On further simplification, one obtains

[TABLE]

where

[TABLE]

Above expression gives rise to a differential equation

[TABLE]

which can be written as

[TABLE]

Introducing: $\tilde{X}_{P}(t)=e^{-iL_{P}t}X_{P}(t)$, we obtain an equation of motion for $\tilde{X}_{P}(t)$

[TABLE]

where

[TABLE]

Solving for one $q$, the solution for $\tilde{X}_{P}(t)$ takes the form

[TABLE]

where

[TABLE]

Using Eqns. (143) and (146), we have (for one $q$)

[TABLE]

where

[TABLE]

Introducing

[TABLE]

and taking the time derivative, one can find an expression

[TABLE]

Using the shifting operator $S_{q}=e^{\frac{\lambda_{q}^{{}^{\prime}}}{\omega_{q}}(b_{q}-b_{q}^{{\dagger}})}$ above expression can be written as

[TABLE]

this implies

[TABLE]

Multiplying by $S_{q}^{{\dagger}}$ on both sides

[TABLE]

solving for $S_{q}^{{\dagger}}\tilde{U}_{q}(t)$

[TABLE]

using $U_{q}(t)=e^{-iH_{q}t}\tilde{U}_{q}(t)$ and $U_{q}(0)=1$, we obtain

[TABLE]

Similarly, for $W_{q}(t)$

[TABLE]

Substituting for $U_{q}(t)$ and $W_{q}(t)$, also using $U_{q}(0)=W_{q}(0)=1$, we obtain an expression for $X_{P,q}(t)$

[TABLE]

Recalling Eq. (142):

[TABLE]

and substituting for $X_{P,q}(t)$, we have

[TABLE]

Above expression can be solved for $q=q^{\prime}=q^{\prime\prime}=q^{\prime\prime\prime}$, because other modes do not contribution in the final trace and do not conserve the particle number for the different modes. On further simplification and doing some algebraic manipulation, one can obtain a simplified expression as

[TABLE]

and similarly for $\mathcal{C}_{2}(t)$

[TABLE]

Considering the factor in parenthesis which is common in both expressions, we have

[TABLE]

where we have defined:

[TABLE]

Using Baker-Campbell-Hausdorff formula:

[TABLE]

similarly for $Y_{q}(t)$,

[TABLE]

Substituting for $Y_{q}(t)$ and $Z_{q}$ in the expression for $X_{P}(t)$, and then substituting for $X_{P}(t)$ in the expressions for $\mathcal{C}_{1}(t)$ and $\mathcal{C}_{2}(t)$, we perform the final trace to obtain the following expressions

[TABLE]

After substituting the expressions for $\mathcal{C}_{1}(t)$ and $\mathcal{C}_{2}(t)$ in the expression for reduced self-energy matrix element and performing the final trace, we obtain

[TABLE]

Applying the continuum of modes and going back to Laplace domain,

[TABLE]

above expression recovers the result obtained for self-energy matrix element within Born approximation, which is Eq. (51) in the main text.

Appendix E Initial Condition

In this section, we will discuss the stationary density matrix and find the initial condition for the operator $\Omega(t)$ given in Eq. (33). The stationary density matrix $\bar{\rho}$ accounts for conditions that accumulate between the system and interactions in the time interval $t\in[t_{0},0]$. In this time interval, $[t_{0},0]$, the stationary density matrix is given by its value at $t=0$ and we can replace $\rho_{S}(t^{\prime})\rightarrow\rho_{S}(t)$ in Eq. (42), to obtain

[TABLE]

Changing integration variable: $\tau=t-t^{\prime}$ to obtain

[TABLE]

one can extend the upper limit of above integration to infinity by setting $t_{0}\to-\infty$, and solving the differential equation to obtain

[TABLE]

and taking Laplace transform on both sides, we obtain

[TABLE]

here Laplace transform is defined as $f(s)=\int_{0}^{\infty}e^{-st}f(t)dt$. We choose an initial condition when exciton is in excited state $\ket{a}$ given by $\rho_{S}(t_{0})=\ket{a}\bra{a}$, and evolves in the presence of pump laser. After performing a secular approximation and using the definition of stationary limit

[TABLE]

one can find all the elements of stationary density matrix operator as:

[TABLE]

where $\Sigma_{z}=\Sigma_{++,++}-\Sigma_{++,--}$. Similarly,

[TABLE]

Recalling, the initial condition for operator $\Omega(t)$ given by Eq. (33):

[TABLE]

and matrix element of interest can be extracted after performing the trace over phonon and photon modes, and substituting for stationary density matrix element $\bar{\rho}_{++}$,

[TABLE]

Dynamics of the operator $\Omega(t)$ is evaluated using Hamiltonian given by Eq. (18), with an initial condition given by above expression which is Eq. (69) in the main text.