Impact of the phonon environment on the nonlinear quantum-dot-cavity QED. I. Path-integral approach

We demonstrate a strong influence of the phonon environment on the coherent dynamics of the quantum dot (QD)-cavity system in the quantum strong coupling regime. This regime is implemented in the nonlinear QD-cavity QED and can be reliably measured by heterodyne spectral interferometry. We present a semi-analytic asymptotically exact path integral-based approach to the nonlinear optical response of this system, which includes two key ingredients: Trotter's decomposition and linked-cluster expansion. Applied to the four-wave-mixing optical polarization, this approach provides access to different excitation and measurement channels, as well as to higher-order optical nonlinearities and quantum correlators. Furthermore, it allows us to extract useful analytic approximations and analyze the nonlinear optical response in terms of quantum transitions between phonon-dressed states of the anharmonic Jaynes-Cummings (JC) ladder. Being well described by these approximations at low temperatures and small exciton-cavity coupling, the exact solution deviates from them for stronger couplings and higher temperatures, demonstrating remarkable non-Markovian effects, spectral asymmetry, and strong phonon renormalization of the JC ladder.


I. INTRODUCTION
A quantum dot (QD) embedded in an optical microcavity offers a technologically accessible platform for the study of quantum-optics phenomena in solid state.This system can be experimentally realized, for example, in self-assembled QDs inside semiconductor micropillar or photonic-crystal cavities [1,2].Some properties of the system, including generation of single photons [3][4][5], are important for development of quantum devices with applications in the field of quantum information.
The repeated mutual conversion of excitation between the QD exciton and the cavity mode, well known as strong light-matter coupling, is typically described by the Jaynes-Cummings (JC) model [6].The strong coupling itself is a linear classical effect leading to polariton formation [7] and observation of the vacuum Rabi splitting in various systems, such as quantum wells inside planar microcavities [8] and atoms inside optical cavities [9], as well as QDs coupled to cavity modes [10,11].However, owing to the fermionic nature of QD excitons, their coupling to bosonic cavity modes introduces a quantum nonlinearity, which manifests itself in anharmonic ladder-like structure of hybrid QD-cavity states.This is know in the literature as the quantum strong coupling regime recently observed QD-cavity systems [12,13] and superconducting circuits [14].This quantum nonlinearity can naturally be measured by means of nonlinear optical spectroscopy [13,15], in which higher rungs of the JC ladder are accessed through an external excitation in a form of a sequence of optical pulses and subsequent selection of different channels of optical nonlinearity [16].
In practice, due to the solid-state environment, the QD can be subject to significant non-Markovian behavior [17][18][19][20].This is mainly a result of the coupling of QD excitons to longitudinal acoustic (LA) phonons via the deformation potential [17,21,22].In specific parameter regimes, when a QD exciton is not very strongly coupled to a cavity mode, the effect of phonons may be addressed perturbatively [23], or via the polaron transformation combined with a perturbation theory [24,25].However, for a stronger QD-cavity coupling, phonons play a more significant role that requires developing non-perturbative techniques [26].This has been recently demonstrated by an exact theoretical approach to the linear optical response [27,28].The linear response is, however, harder to observe experimentally, as compared to a nonlinear optical polarization that can be reliably measured by means of a heterodyne spectral interferometry [15], and is also a suitable tool for the study of quantum nonlinearities and the quantum strong coupling regime.Alternatively, quantum nonlinearities can be studied by analyzing photon counting statistics in photoluminiscence measurements [29].
The four-wave mixing (FWM) polarization of a QDcavity system has been the focus of many experimental and theoretical works but the exciton-phonon interaction has either been neglected or treated perturbatively [13,[30][31][32].In the absence of a cavity, the effect of the phonon dephasing measured in the FWM response [33] has been addressed theoretically for ensembles of isotropic [34] and anisotropic QDs, also accounting for real phonon-assisted transitions [35].One important advantage of using a FWM scheme is a possibility to use a two-dimensional hyperspectral imaging that reveals the information about coherences in the system [36].The effect of phonons on two-dimensional spectra for a nonlinear optical response has been studied in the absence of a cavity using a path integral (PI) based approach [37][38][39] in the context of energy harvesting.In addition to the FWM, other quantities of experimental relevance, including photoluminescence and photon indistinguishability [40,41], have been recently studied in QD-cavity systems [42,43], using a PI-based approach to include a phonon contribution.With an increasing interest in understanding of non-Markovian effects and the need to study non-perturbative regimes of more complex models, there has been a rise in popularity of such treatments, including applications in the context of QDs [44][45][46].PI-based approaches offer a possibility to achieve numerically exact solutions for a QD-cavity system in a phonon environment and are capable of addressing arbitrary parameter regimes, where approximations fail.Such solutions typically require Hilbert spaces, which rapidly grow in size.However, the environment degrees of freedom can be reduced dramatically by means of introducing a memory kernel [47].The memory kernel takes into account all the necessary information introduced by temporal correlations in the system evolution as a result of its interaction with the environment [48].The complexity of the problems one can address is generally limited by the computational memory and speed.The approaches based on the early real-time PI methods developed in Ref. [47,49] typically depend on the combined use of Trotter's decomposition and Feynman-Vernon influence functional [50].Instead of the latter, Ref. [27] uses a linked-cluster, or cumulant expansion [51], that formed the basis of this work, as it is a suitable approach for treating a linear coupling betweeen the system and a continuum of bath modes.
In order to tackle the computational costs and improve efficiency, a number of further PI optimization schemes have been proposed.Early approaches use basic selective filtering methods, such as those which filter out path segments with negligible weight by Monte-Carlo importance sampling [52,53].More recently, an idea to spread the temporal correlations over path segments of increasing length, whilst their contribution is reduced, was formulated [54].Other recent works have developed an optimization scheme, which is based on rewriting the PI as a tensor network (TN), where the propagator is represented by a matrix product state (MPS).Compression of large tensors can be achieved by a selective truncation that relies on singular value decomposition [55].The propagation can then be done very efficiently using existing tensor network algorithms.This approach has seen further development [56,57] and has been recently combined with mean-field theory [58] to address systems with large Hilbert spaces.Another work has combined the approach developed in Ref. [55] with a different MPS-based algorithm, introducing a time-discrete quantum memory as a second non-Markovian reservoir, resulting in a quasitwo dimensional TN [59].
In addition to PI-based approaches, there are other numerically exact methods.These include solving hierarchical equations of motion [60][61][62] and numerical renormalization group approaches [63,64].An implementation of the density matrix (DM) renormalization group [65] is based on the optimization of matrix product states [66], which also uses singular value decomposition.One further optimization implements a mapping of the the environment onto a one-dimensional chain with effective nearest-neighbor interactions [67], which makes a subsequent application of time-adaptive DM renormalization group algorithm [68] very efficient.
In the present work, we focus on the FWM polarization of a QD-cavity system linearly coupled to LA phonons.We develop an asymptotically exact semianalytic approach which is a generalization of the method presented in Ref. [27], with key ingredients being the Trotter decomposition [69,70] and linked-cluster expansion [71,72].This method is an explicit and physically intuitive PI-based approach allowing a number of useful analytic approximations [27,28], also presented for the case of the FWM in this and the follow-up paper [73].Since the approach is non-perturbative, it allows us to explore regimes of comparable exciton-cavity and excitonphonon coupling strengths, when the similarity of the system and environment timescales opens up a possibility of phonon-assisted transitions between different polariton states [27,74,75].Although the treatment of optical nonlinearities increases the complexity of the approach, it remains physically intuitive and computationally straightforward.Furthermore, it allows us to address all possible channels of optical excitation and measurement and provides perspectives to a further application of the method to higher-order optical nonlinearities, arbitrary elements of the density matrix, and other quantum correlators and physical observables.
We demonstrate that the exciton coherent dynamics in semiconductor quantum dots can be significantly modified by the phonon environment, showing a remarkable non-Markovian behavior.Based on the developed technique, we propose two computationally simple approaches to the regime of small exciton-cavity coupling.One of these approaches has a fully analytic form while the other reduces the general method to a simple matrix multiplication.Being not limited to perturbative regimes, our microscopic approach can deal with arbitrary temperatures, as well as with situations where the interactions of the exciton with the cavity mode and the phonon environment are comparable.In the latter case, the phonon cloud around the quantum dot is unable to adiabatically adapt to a varying optical state, resulting in a non-Markovian dynamics and phonon-assisted transitions between states of the Jaynes-Cummings ladder.We observe a spectral asymmetry, a modification in anharmonicity of the ladder, and even deviation from the ladder-like structure, which becomes more pronounced with increasing temperature.By extracting complex fit parameters we analytically characterize the long-time behavior of individual transitions which contribute to the FWM signal and demonstrate the non-Markovian nature of the latter.
In this section, we describe the formalism of our asymptotically exact semi analytic approach to the FWM dynamics of the QD-cavity system surrounded by a phonon environment.We first introduce in Sec.II A the system Hamiltonian and the master equation containing Lindblad dissipators for both the QD exciton and the cavity mode.We also discuss there in detail optical excitation and measurement channels for observation of the FWM response of the system and comparing them, where appropriate, to those of the linear response.We then focus in Secs.II B and II C on two key elements of our approach: Trotter's decomposition and linked-cluster expansion.Finally, in Secs.II D, II E, and II F we show how these elements work together to accurately describe the FWM dynamics of the system.This includes the full rigorous L-neighbor approach presented in Sec.II F and a number of useful approximations following from it, such as nearest-neighbor and polaron approximation, presented in Secs.II D and II E, respectively.
A. System, its excitation, and measurement of the linear and FWM polarizations We consider a composite system, which presents a combination of two analytically solvable models, the JC model describing the exciton-photon coupling and the independent boson (IB) model [76] taking care of the exciton-phonon interaction.The composite system is, however, not solvable analytically, to the best of our knowledge.Its full Hamiltonian is given by where Here, the fermionic (bosonic) creation operator d † (a † ) corresponds to the exciton (cavity) mode with real frequency Ω x (Ω c ), and g is the exciton-cavity coupling strength.We work in the units of = 1 throughout this paper.The bosonic operator b † q creates a phonon with momentum q and frequency ω q , where q = |q|.The coupling of the exciton to the LA phonon mode q is given by the matrix element λ q , which has a general symmetry property: λ * q = λ −q .The explicit form of λ q in case of a linear phonon dispersion is given by Eq. (B1) of Appendix B.
The exciton-phonon Hamiltonian Eq. ( 3) contains a single exciton mode, so we assume that other exciton states of the QD do not couple to this state via phonons.Such a diagonal form of the exciton-phonon interaction is known in the literature as the IB model which leads in the absence of the cavity to exciton pure dephasing [17].While this model can adequately describe a measured quick relaxation of the FWM polarization [34], it does not include any mechanisms of experimentally observed [77] density or polarization relaxation at longer times.To rectify on this, we add to the model a phenomenological exponential decay rate γ x of the exciton state, which reffers to the population decay time T 2 .This can be further refined in a rigorous way within the single-mode (i.e.diagonal) exciton Hamiltonian by mapping the offdiagonal phonon coupling to other exciton states onto a diagonal quadratic interaction [71,78] that accounts for the pure dephasing T * 2 .We also introduce a decay rate γ c of the electromagnetic eigenmode of the cavity, which is strictly the negative of the imaginary part of the complex frequency of the cavity mode.Both γ x and γ c are introduced consistently via the Lindblad dissipation formalism and do not include any phenomenological pure dephasing.
The evolution of the composite system between the pulses is modeled by a master equation which has a Lindblad form: where L is the Liouvillian superoperator of the composite system.Here, the dephasing of the cavity and the exciton modes is represented by a Lindblad dissipator defined as In the dipole approximation, the excitation of the system by a sequence of ultrashort pulses is given by an interaction operator in which the external classical field due to each pulse is represented by a delta function peaked at t = t j and a pulse area µ cj E j , with the index j = I, II, . . .labeling the pulses.This corresponds to a situation in which each pulse duration is much shorter than the characteristic timescales of the system.The effective dipole moment µ cj describes the coupling strength of the mode c † j to the pulse field E j .The excitation can take place via the cavity mode (c † j = a † ) or the exciton mode (c † j = d † ), or a combination of both [28], and is further referred to as the excitation channel.The time evolution of the system due to each pulse is then given by a product of exponential operators [13,16].The DM under the action of the interaction Q(t) can then be expanded in all orders of the excitation field as a series of nested commutators.The phase of E j permits selection of a particular channel in the optical polarization, which is linear or nonlinear in terms of the powers of E j .
For the linear polarization, we need to consider only a first-order effect on the DM by a single pulse, which is Pulsed excitation scheme showing a sequence of ultrashort laser pulses with a positive delay time τ .The first pulse with the pulse area µc I EI and action Q (I) on the DM arrives at tI = −τ .After delay τ , two identical pulses with the resulting pulse area µc II EII and the total action Q (II) arrive together at tII = 0.The observation time is denoted by t.
proportional to E * I , which is due to a single pulse at t I = −τ .For a nonlinear polarization, we focus on the well-known case of degenerate FWM in which three-pulse excitation is effectively reduced to two pulses, at t I = −τ and t II = 0, with τ being the delay time between the pulses.Figure 1 schematically illustrates this pulsed excitation scheme.The action of the first pulse, denoted by Q (I) , is described by Eq. ( 7), the same way as in the linear polarization.
The second and the third pulses, both arriving at t = 0, are considered as a single pulse, with the action denoted by Q (II) .The FWM signal is proportional to E * I E 2 II in the lowest excitation order, and a third-order correction to the DM is given by describing the action Q (II) of the second pulse.We emphasize that the excitation channels for pulses I and II need not be the same, so one can have c I = c II .Having discussed the selection of the desired channels for both the linear and FWM signals, in what follows, we will ignore the constant prefactors −iµ * cI E * I and (iµ cII E II ) 2 , respectively, in Eqs.(7) and (8), consistent with [16].
The total optical polarization is given by: where the expectation value is taken with respect to the states of the exciton-cavity system, as well as phonon degrees of freedom, which is expressed by double brackets.Similar to the excitation, the response is detected through the coupling to an external field, and O represents an annihilation operator of the observation channel: either cavity O = a, or exciton O = d, or a combination (i.e. a linear superposition) of both [28].
It is convenient to expand the DM into the basis of the JC model, where ρ ηξ (t) includes the phonon degrees of freedom.In order to describe the linear and FWM polarizations, in addition to the ground state |0 , the following reduced basis [13] is used in Eq. ( 10): For the linear polarization, only the ground state |0 , the exciton state |1 and the single photon state |2 are required.In order to describe a third-order nonlinearity, two states of the second rung, are also needed.These are state |3 containing single excitations both in the exciton and cavity modes and the two-photon state |4 .As we are dealing with a linearly polarized light, the system also has in principle a biexciton state involved in the FWM dynamics.However, owing to the biexciton binding energy (which is typically much larger than g), this state is detuned out of the energy range in focus (which is of the order of g) and can be ignored.The effect of phonons on a nonlinear response of an excitonic system with a possibility of exciton-biexciton transitions was explored in [79], in the absence of a cavity.Using the expansion Eq. ( 10), the master equation ( 5) can be expressed in terms of the Liouvillian written as a matrix and the DM written as a vector [13]: When expressed as a vector in the basis of states given by Eq. ( 11), the DM takes the following form [16]: The superscripts 0 and I (II) are used to distinguish quantities describing the system, respectively, before the pulses and after the application of the first (second) pulse.They refer to a specific size of the reduced basis needed to describe the signal.The initial DM has only one component ρ 00 = ρ ph with phonons in thermal equilibrium, where k B is the Boltzmann constant and T is the temperature.Before the first pulse (at t I = −τ ) and immediately after it (or after both pulses if τ = 0), the DM is factorizable into the JC and phonon parts.In fact, the excitation pulse acts only on the JC component of the DM.Since the pulse is infinitely short, the phonon component of the DM remains unaffected at the time of excitation.At any later time these parts of the DM are entangled, and the phonon bath surrounding the QD is no longer in equilibrium.
When the DM is factorizable or any effects of phonons are neglected, the Liouvillian L JC of the exciton-cavity system can be diagonalized analytically -a detailed procedure can be found in the supplement of Ref. [16] where the FWM and higher-order responses were studied with inclusion of many rungs of the JC ladder, though in the absence of exciton-phonon interaction.Here, we exploit this procedure in order to solve the full problem, in which the exciton-phonon interaction is taken into account.

Linear response
A linear polarization can be produced by exciting the system with a single pulse.With ρ (0) being a 1×1 vector [see Eq. ( 13)], the action on the system of an ultrashort pulse at t I = −τ can be represented by a vector Q (I) .Then Eq. ( 7) takes the form: with where the index x (c) refers to the exciton (cavity) excitation channel.The effect of Q (I) is to produce a singleparticle excitation (in either the exciton or the cavity mode), converting the element |0 0| of the DM into |0 1| or |0 2| in the decomposition Eq. (10).The Liouvillian for the exciton-cavity system after the first pulse has the following form: Following the evolution of the DM, from ρ (I) (−τ ) to ρ (I) (t), a linear response is measured in a specific observation channel O.With ρ(t) and O represented as vectors, Eq. ( 9) takes the form of a dot product where, for the linear polarization, and the thermal expectation value in Eq. ( 18) is taken over the phonon bath only.

FWM response
Starting with the initial excitation, given by Eq. ( 15), followed by a delay dynamics, in which the system evolves from t = −τ to t = 0, the FWM response is produced by applying another pulse (consisting of two identical simultaneous pulses) at t II = 0.With ρ (I) given by a 2 × 1 vector in Eq. ( 13), the action of the second pulse can be represented by a matrix Q (II) , so Eq. ( 8) becomes where In terms of the decomposition Eq. (10) (13).The Liouvillian matrix then takes the form: The FWM polarization is given by and is measured (at t > 0) in an observation channel O (II) , which can be either exciton or cavity, or a linear combination of the two.

B. Trotter's decomposition
We now consider the evolution of the system between and after the pulses which is described by the master equation ( 5), or its matrix analogue Eq. ( 12).The first key element of our asymptotically exact approach to the system dynamics is the Trotter decomposition.It is based on Trotter's theorem [69] for two non-commuting operators A and B: where ∆t = t/N and N is an integer.In our case, the IB and JC models are two exactly solvable parts of the system, which are described by non-commuting operators.We therefore set A = −iL IB and B = −iL JC , the two non-commuting operators of the full composite system, described by the Liouvillian We then use Eq. ( 25), in order to separate the evolution of the full system into discrete time intervals, and assume independent evolution of the IB and JC parts within each time interval, exploiting the exact analytic solvability of each part.The time interval between the pulse time t 0 (where t 0 = −τ or t 0 = 0) and the observation time t is split into N discrete time steps, which need not be equidistant.Concentrating on the n-th time step, between t n−1 and t n , with ∆t = t n − t n−1 , and assuming that ∆t is small, the evolution of the DM over this time interval can be approximated, in line with Eq. ( 25), as To account for the effects of L IB and L JC , we introduce, respectively, a matrix operator W and a matrix M (both defined below), such that Eq. ( 27) can be written as The exciton-phonon dynamics, represented by W , follows a unitary evolution that can be fully described by the Hermitian Hamiltonians H IB and H ph : where is the interaction representation of the phonon coupling V defined in Eq. ( 4).Due to the diagonal form of the exciton-phonon interaction, W and W † are diagonal matrices, and owing to the Hermiticity of the IB model, The non-Hermitian dynamics of the exciton-cavity system is in turn described by the matrix M defined as where matrix L JC depends on the basis used and for the linear and FWM polarizations is given, respectively, by Eq. ( 17) and Eq. ( 22).
Introducing the components of the DM vector, the approximate Eq. ( 28) can be written as where the phonon operators W inin and W † inin are the matrix elements of the diagonal matrix operators W (t n , t n−1 ) and W † (t n , t n−1 ), respectively, the numbers M inin−1 are the matrix elements of M , and ρ in is the i nth component of the vector ρ(t n ) representing the DM at the time moment t n .Note that the subscript n is included in the index i n labeling matrix and vector components, in order to keep the information about this time moment and the selected time interval (between t n−1 and t n ) of the system evolution described by Eq. ( 32).Physically, i n labels the quantum state of the full system at time t n on the selected path of its evolution.In the case of non-equidistant time steps, W and M depend also on n, which is omitted for the brevity of notations.Now, we introduce a specific time-ordering operator T , which allows us to move the operator W † inin next to W inin : In doing so, the order of phonon operators should be preserved.The operator T ensures that all phonon operators in W (W † ) stand to the left (right) of ρ in normal (inverse) order.We further introduce an operator which includes the effects of both W and W † , with the help of T and the interaction This allows us to write Eq. (33) as Introduced in Eq. ( 35), operators V (±) (t ) are the same as the operator V (t ), defined by Eq. ( 30), with the only difference that under the action of T all phonon operators in V (+) (V (−) ) stand to the left (right) of ρ in the normal (inverse) order.Vectors α and β are related, respectively, to the left and right sides of the DM and have Boolean components α i and β i indicating whether or not the exciton-phonon coupling is affecting the given element ρ i of the DM from either side.Within Eq. (36), Y in has effect on M inin−1 ρ in−1 from the left (right) only if α in = 1 (β in = 1), which takes place when the corresponding component of the DM contains exciton.Otherwise α in = 0 (β in = 0), and there is no effect of Y in .Note that α and β are constant vectors, independent of the time step n, but are different for different orders of optical nonlinearity.Their explicit form for the linear and FWM polarizations is provided in this section below.Applying Eq. ( 36) to the full dynamics of the DM, from the initial time t 0 of pulsed excitation to the final time t = t N of signal observation, and using the excitation and measurement conditions, as detailed in Sec.II A, the optical polarization takes the form where in the last line, a thermal expectation value is taken over the phonon ensemble.Clearly, the phonon contribution is present only in the time-ordered product of the elements Y in consisting of V (±) operators.This expectation value will be evaluated in Sec.II C with the help of the linked-cluster theorem [76].Equation ( 37) is valid for any channel of optical nonlinearity, including linear and FWM response, provided that the system is excited by ultrashort optical pulses.Below we provide details of how to use Eq. ( 37) for calculating the linear and FWM polarization at zero delay (τ = 0).The case of a non-zero delay between the pulses (τ = 0) is considered in detail in Appendix A and also in Sec.II D.

Linear polarization
The linear polarization P Lin (t) can be extracted from Eq. (37), by using M = M (I) with L JC = L (I) JC given by Eq. ( 17), O = O (I) given by Eq. ( 19), and Q = Q (I) given by Eq. ( 16).The phonon contribution is taken care of by Y in , defined by Eq. ( 34), with α = α (I) and β = β (I) , where as follows from the DM components involved, which are given by Eq. ( 13).

C. Linked-cluster expansion
The second element of our approach to the system dynamics is the cumulant, or linked-cluster expansion.It allows us to address the exciton-phonon interaction exactly, by providing an analytic evaluation of the expectation value of the time-ordered product of operators in Eq. (37).As a result, any explicit phonon dependence is removed at a cost of introducing temporal correlations between the states of the system at different time steps.
The linked-cluster theorem [76] is valid for arbitrary operators and a wide class of time orderings, including the normal (T ) and the inverse (T inv ) time orderings.The specific form of time ordering T introduced above can be considered as a combination of T and T inv , which are distinctly separated: In fact, differently ordered operators stand on either side of the DM and do not mix.As the two classes of time ordering do not mix in T , the linkedcluster theorem holds for this time ordering, resulting in cumulants similar to those which appear in the IB model.
Within the IB model, only second-order connected diagrams contribute to the cumulant.Applying the linkedcluster theorem, the expectation value of the product of Y i operators in Eq. ( 37) takes the form Kimin (m, n) , (41) where and Ṽi (τ ) is defined by Eq. (35).Shown above is a general result for an arbitrary channel of nonlinearity.Depending on which channel is considered, vectors α and β take a particular form, see e.g.Eqs.(38) and (40).
Owing to the linked-cluster expansion, the phonon contribution to the system dynamics is thus expressed in a form of two-time correlations in system variables determining the path, indexed by i n .The possible realizations (or paths) are indicated via combinations of the basis state labels i n for all values of n on the time grid.In the optical polarization Eq. ( 37), all possible realizations are summed over.
With a general derivation given in Appendix A, Eq.( 42) becomes for m n, and is symmetric with respect to the interchange of indices: Within Eq. ( 43), each cumulant element is given by: where is the phonon propagator and is the Bose occupation number.Substituting Eqs. ( 41), (43), and (45) into Eq.( 37) allows us to find the optical polarization.
Note that we have used so far an arbitrary time grid, which is not necessary equidistant.Below, for clarity of presentation, we use an equidistant time grid, with a constant time step ∆t.As we would like the time evolution to be composed of identical memory kernels (visualized by L-shapes, see below), an equidistant grid is also more relevant because the cumulant function is a function of the difference between the two times t 1 and t 2 .In this case, all the cumulant elements depend on the difference l = |n−m| only and can be written as K nm = K mn = R l , with R l found recursively via starting from l = 2 and using R 0 = K(∆t).Here the cumulant function of the IB model is defined as [cf.Eq. ( 45)].Then for m n, Eq. ( 43) simplifies to where l = m − n 0. In the opposite case of m n, one can use the symmetry Eq. ( 44) to obtain Kimin (m, n) = K i m+l im (l) with l = n − m and the same K ij (l) defined in Eq. (50).
If the phonon bath has a memory time τ IB (introduced in Ref. [27], see also Appendix C), it is not necessary to take into account all R l on the time grid up to l = N , but it is sufficient to cover only a portion of the twodimensional time grid including a distance of about τ IB /2 from the main diagonal (see Sec. II F).In this case the number of cumulant elements R l which are taken into account is limited to some finite value l L. The number of time steps in the memory segment of the system evolution, or the number of "neighbors" L is chosen in our calculation as to satisfy a criterion ∆t τ JC , τ IB , which is imposed by the Trotter decomposition.Here τ JC is the timescale for the JC dynamics, see Appendix C. Below we first consider in Sec.II D the L = 1 case called a nearest-neighbor approximation and then proceed in Sec.II F with a general case of L 1 which we call an L-neighbor approach.

D. Nearest-neighbor approximation
The theory presented so far has focused on the case of zero delay between optical excitation pulses.For the purpose of illustration of the full method, we consider here both cases of zero and non-zero delay, illustrating and comparing them in the nearest-neighbor (NN) approximation.The latter is achieved by setting L = 1 so that only |n − m| 1 are kept in the double sum in Eq. ( 41), and the DM or the optical polarization can be expressed at any time simply as a product of matrices.This approach is valid for a relatively small exciton-cavity coupling strength g, for which τ JC τ IB , and breaks down as the two timescales become comparable.Owing to its simplicity, the NN approximation serves as a clear illustration of the Trotter decomposition method in general and provides a basic start point for the full L-neighbor (LN) solution presented in Sec.II F below.It also allows us to clearly demonstrate how one can treat an arbitrary sequence of excitation pulses, at the same time taking into account any memory effects and the phonon contribution to the evolution of the system between the pulses.
In the limit τ JC τ IB , the full cumulant in Eq.( 41) can be approximated by a reduced cumulant which includes only the diagonal (m = n) and the NN (m = n ± 1) terms.For the linear or FWM polarization with τ = 0, the cumulant for a given realization (or a path) then takes the form: where K ij (l) is defined in Eq. (50).For the linear polarization, using the explicit form of α (I) and β (I) given by Eq. (38), this reduces to where δ nm is the Kronecker delta.Note that Eq. ( 52) is the complex conjugate of the expression given by Eq. ( 19) of Ref. [27], which is obtained here in accordance with the phase selection determined by Eq. ( 7).The optical polarization can then be written as a product of matrices, where and t = N ∆t.The time step ∆t satisfies the condition ∆t ≤ τ IB , and for short times (0 < t < τ IB ) is chosen as ∆t = t/2, so that the NN approximation takes all cumulant elements into account.Eqs. ( 53), (54), and ( 55) are valid both for linear and FWM polarization at τ = 0.For the FWM (linear) polarization, one should use O, Q, α, and β with index II (I), except Q which is given by Eq. ( 39) in the case of the FWM with τ = 0.
For the linear polarization, Eq. ( 53) reproduces the result obtained in Ref. [27].Physically, M Q in Eq. ( 53) represents the state of the system following a pulsed excitation and a subsequent propagation by one time step under the influence of L JC .The full evolution of the system by a single time step is carried out through a multiplication with the matrix G, which is repeated N − 1 times.The phonon contribution at the final step of evolution is included in the matrix E, and the optical polarization is then obtained by taking the dot product with O, in this way projecting the DM onto a selected observation channel.
For the FWM polarization with τ > 0, the cumulant Eq. ( 51) modifies to include in the dynamics both the delay time τ = N I ∆τ and the observation time t = N II ∆t, with the number of discrete steps N I and N II , respectively.The time steps satisfy the condition ∆t, ∆τ ≤ τ IB .Equation ( 53) then modifies as follows where and matrices G (ζ) and K (ζ) (l) are defined as before, respectively, by Eqs. ( 54) and ( 50), with the upper index ζ (taking the values I or II) added to the matrix M and vectors α and β, see Eqs. ( 17), (22), and (31), as well as Eqs.( 38) and (40) for their definitions.Matrix is also defined as in Eq. ( 55), now with the upper index II added.The range of the matrix indices is different in different regions: G (I) and M (I) are 2 × 2 matrices; G (II) , M (II) , and E (II) are 6 × 6 matrices, while G (I−II) is a 6 × 2 matrix.
The specific cumulant elements, which are required for the NN approximation and appear in Eqs. ( 50) and (58) are obtained in the same way as R l in Eq. ( 48), and take the form Note that they are all expressed in terms of the cumulant function K(t) given by Eq. ( 49).Clearly, within the NN approximation there are only five distinct cumulant elements given by Eq. ( 59) which must be considered.
The NN approximation Eq. ( 56) for the FWM polarization with τ > 0 is illustrated in Fig. 2, indicating three different regions in the time domain, labeled with (I), (II), and (I-II).Different cumulant elements R j defined in Eq. ( 59) are distinguished by color.The correlations in memory, arising from the phonon contribution, are limited to the indicated L-shaped regions.Note in particular that the (I-II) region, where the reduced basis size changes, includes path segments connecting states of the system evolution both in the delay and the observation time regions.

E. Long-time asymptotics
A much better understanding of system's behavior can be achieved by separating the dynamics which occurs on different timescales, τ IB τ JC .We provide two different ways to describe the long-time behavior analytically: (i) applying a multi-exponential fit to our numerical solution and (ii) using the polaron approximation (PA).This section is devoted to the PA, while the multi-exponential fit approach is introduced at the end of Sec.II F below.For the review of the PA in greater depth and its behavior with the delay time, see Ref. [73].
The PA is obtained by taking the long-time limit of Eq. ( 49), K(t) → −S + iΩ p t.For the linear and τ = 0 FWM polarization the PA is given by (see Ref. [73] for the derivation) where S is a diagonal matrix, given by S (I) = diag(S, 0) for the linear and S (II) = diag(S, 0, 0, S, S, 0) for the FWM polarization, and LJC is the polaron transformed JC Liouvillian, obtained from L JC provided in Eqs. ( 17) and ( 22) by replacing the exciton transition energy Ω x and the coupling strength g with, respectively, a polaronshifted exciton energy Ωx and an effective coupling strength g, which are given by The polaron shift Ω p and Huang-Rhys [80] factor S have the following explicit form: and are determined by the phonon spectral density J(ω) = q |λ q | 2 δ(ω − ω q ).For QD systems, considered in this work, the phonons spectral density has a superohmic form J(ω) = J 0 ω 3 e −ω 2 /ω 2 0 , where J 0 and ω 0 are constants determined by the material parameters of the QD (see Appendix F of Ref. [27] for details and derivation).
For τ > τ IB , the PA for the FWM polarization can be obtained from Eq. ( 56) (see Ref. [73] for the derivation) and becomes where with and Q (II) c defined in Eq. ( 21).

F. Full numerical L-neighbor approach
Here we present a full LN approach, already mentioned at the end of Sec.II C, which is based on the scheme proposed in Ref. [27] and can be understood as a "forwardmemory" approach (see below for details).For clarity of presentation, we return back to linear and zero-delay FWM polarizations -the general case of non-zero delay is derived in Appendix A which also provides details of the formalism presented here.
The LN approach is a generalization of the NN approximation, given by Eq. ( 53), and consists of a recursive generation of tensors F (n) , starting from n = 1, and ending up with n = N , the latter determining the optical polarization: Tensors F (n) of rank L are generated with the help of a constant tensor G of rank L + 1 which is the LN analogue of the matrix G, given by Eq. ( 54), that appears in the NN approach.Calculation of the tensor G requires only L + 1 distinct cumulant elements K ij (l) defined in Eq. (50).As in Eq. ( 53), for the FWM (linear) polarization, one should use vectors, tensors, and matrices with the upper index II (I), except Q which is given by Eq. (39) in the case of the FWM with τ = 0. We shall introduce a short-hand notation, where we first specify the excitation channel and then the measurement channel (e.g.c-x refers to the situation where the system is initially excited in the cavity mode and the response is measured via the excitonic mode).In Eq. ( 68) h = p = 6 (h = p = 2) for the zero-delay c-c and c-x FWM (x-x, x-c FWM and linear) polarization.The time step is set to ∆t = τ IB /L for t τ IB and at earlier times (t < τ IB ) to ∆t = t/(L+1), in G: Diagrams showing how the system variables determining the path are linked within G (blue) and F (red) of Eq. ( 66) for L = 5.G contains path segments connecting the time of observation (indexed by j) to itself, as well as to five other nearest time intervals (indexed by in).These links are weighted by appropriate cumulant elements if the interaction is present or by zero if there is no interaction and form a memory kernel.In F, all system variables that describe L-step segment along some path, are linked together.This tensor contains a future 'history' of exciton-phonon interactions that extends over five nearest timesteps.the latter case taking all possible cumulant elements into account.This allows one to resolve the initial dynamics on finer timescales and capture the fast initial decay of the polarization.
Physically, Eq. ( 65) forms the initial DM, by applying a pulsed excitation to the ground state of the JC system and adding a propagation by one time step under the influence of L JC [see Eq. ( 31)], without contribution from the phonon dynamics yet.The subsequent evolution of the DM is described by Eq. (66), representing each single time step in the form of a tensor product, expressed as a summation over j.The final result expressed by Eq. ( 68) is obtained by taking a dot product O • ρ(t) , in accordance with Eqs. ( 18) and ( 23), projecting the DM onto a selected observation channel.
The tensor G, given by Eq. ( 67), is a memory kernel responsible for phonon-induced correlations which are schematically demonstrated in Fig. 3 for both G and F. The tensor G contains the full information required to propagate the system by a single time step and includes path segments connecting the current time interval to L nearest intervals and to itself (Fig. 3, blue).The tensor F (n) contains the information about the state of the system and its interactions with the phonon environment in the form of two-time correlations in system variables determining the path (Fig. 3, red).It includes all possible path segments within the memory kernel.In the literature, G and F, or their analogues, are sometimes referred to as the propagator tensor and the augmented density tensor, respectively (see e.g.[47,49,55]).
As already mentioned above, we follow the forward memory LN approach introduced in Ref. [27]. Figure 4 illustrates the method for L = 5, showing a time grid representing temporal correlations in system variables.The L-shaped regions (bounded by dashed lines) contain contributions within the memory kernel.A product of (II) 4. A time grid for L = 5 corresponding to the evolution of the system between the pulse time and the observation time t.The evolution along this grid consists of L-shaped regions, which are bounded by the dashed lines.The time steps within each L-shape are linked by an element of the tensor G. Grid squares are colored to distinguish different cumulant elements Eq. ( 50), which characterize the strength of non-local interaction between system variables for a pair of time steps within the memory.Cumulants (elements of the grid) not included in the path sum (as they are beyond t) are labeled by ×.
terms in each such region forms a rank-6 tensor G i5...i1j .As it is clear from direct L shapes, the memory kernel includes correlations with the future time steps, rather than the past, and any cumulant contributions after the observation time (labeled by crosses) are discarded.In Eq. ( 68), representing the physical observable extracted from the tensor F, the index p is chosen in such a way that any cumulant contributions which fall outside the range of the propagation are removed.In other words, the choice of the index p is dictated by the fact that L last L-shape regions are truncated.Choosing p = 2 or p = 6 for the FWM polarization (both corresponding to the cavity components of the DM) ensures that there is no phonon contribution beyond the observation time t, in accordance with the vanishing second and sixth elements of both vectors α (II) and β (II) , see Eq. (40).In fact, this corresponds to paths beyond time t with no change of the DM, since phonons do not interact with the cavity.
It has been demonstrated in Ref. [27] that the linear polarization shows a bi-exponential behaviour at long times.Expecting a similar property of the FWM polarization, we introduce a multi-exponential form of a fit function, which we use in Sec.III below for approximating the full numerical data at long times, t τ IB .Here h = 6 (h = 2) for the FWM (linear) polarization, corresponding to six (two) optical transitions involved in the coherence dynamics.The amplitudes A j and the corresponding frequencies ω j are complex-valued and found numerically by fitting the numerical data at long times.Note that when excited via the cavity channel, the system reaches second rung of the JC ladder, so in fact six transitions take part in the FWM dynamics, and h = 6.Excitonic excitation can only produce ground state to first rung coherences in the FWM, therefore only two possible transition frequencies contribute, so h = 2, the same as in the linear polarization.

III. RESULTS
In this work we use the exciton-phonon parameters of InGaAs QDs, same as in Ref. [27].These parameters are in turn taken from Refs.[71,72]: D c − D v = −6.5 eV, ρ m = 5.65 g cm −3 , v s = 4.6 km s −1 , l = 3.3 nm.They describe, respectively, the difference in deformation potentials of the conduction and valence bands, mass density, speed of sound in the QD material, and the electron and hole Gaussian confinement radius.The cavity parameters are extracted from the FWM experiment [13] with QDs embedded in a micropillar cavity: g = 50 µeV (corresponding to g = 39 µeV at T = 50 K), γ c = 30 µeV, and γ x = 2 µeV.We also use in all results a zero effective detuning, Re(ω x − ω c ) = 0, in which the polaron shift Ω p is compensated by the energy difference between the bare exciton and the cavity modes, see Eqs. ( 61) and (62).As for the coupling strength g, in addition to the above value of g = 50 µeV used in Sec.III A and corresponding to the experiments [13,31], we explore in Secs.III B and III C and in Appendices D and E also the regime of a much stronger coupling, with g = 0.3 meV and g = 0.8 meV.These large values of the coupling strength were shown to be achievable in modern experiments [12,75,[81][82][83].

A. Small g regime: Comparison of approaches
In this section we explore the regime of relatively small exciton-cavity coupling strength g = 0.05 meV, for which τ JC τ IB .Still, the exciton and the cavity mode are strongly coupled.In this regime, both the NN (corresponding to L = 1, see Sec.II D) and the PA (see Sec. II E) provide good approximations to the exact optical polarization calculated in the full LN approach (see Sec. II F).The latter is also analyzed in terms of the multi-exponential fit described at the end of Sec.II F.
The convergence of the LN result to the exact solution is studied in detail in Appendix E. In this section we use only L = 9 (the numerical results are referred below to as just LN).This value of L is sufficient for the coupling strengths used in this paper (with root mean square deviation from the exact result of up to 0.03-0.3%for the linear polarization), as demonstrated in Appendix E.
While the main purpose of this paper is to study the FWM polarization, we also include here some results for the linear polarization, in order to provide a broader comparison.This allows us to see changes of the phonon contribution to the polarization, induced by the third-order optical nonlinearity.For a more detailed study of the linear polarization in this system see Ref. [27].As for the FWM, we focus here on the case of zero delay τ = 0 for clarity.A generalization of the theory for arbitrary delay times is available in Appendix A and is studied in more depths in Ref. [73].
To see how linear and nonlinear effects manifest themselves in different excitation and measurement channels, we focus here on (i) the FWM polarization in the x-x channel when the optical excitation and measurement are performed through the QD exciton state, and on (ii) both the linear and FWM polarizations in the c-c channel when the excitation and measurement are done via the cavity mode.Other combinations of excitation and measurement channels in the FWM and linear polarizations are explored, respectively, in Appendix D and Ref. [28].Note that the dynamics of the x-x FWM polarization is essentially the same as the x-x linear polarization, studied in Ref. [27], since only the first rung of the JC ladder is involved, as it is clear from the form of the excitation operator Q (II) x given by Eq. ( 21).However, this simpler coherent dynamics can still be measured [15] as a nonlinear optical response of the system.As for the c-c channel, the linear and the FWM polarizations are entirely different, as they involve different quantum transitions, as discussed below in more depths.
To begin with, we consider the situation when the system is excited and the response is measured via the excitonic mode (Figs. 5 and 6).In this case the second-rung coherences are absent in the dynamics and the results reduce to the linear polarization, as mentioned above.Note that in the absence of the cavity, the zero-delay FWM and the linear response are identical up to a phase factor P FWM (t, 0) = −P * Lin (t), and the full effect of the linear exciton-phonon coupling Eq. ( 4) can be taken into account analytically, which is know as the IB model.Moreover, even for arbitrary delay times, the FWM response can be expressed entirely in terms of the linear response using an analytic relation [35] In the presence of the cavity, the linear and the zero-delay FWM polarizations in x-x and x-c channels (where the second rung is not directly excited) are also identical up to a phase factor in the full calculation.x-x FWM polarization at τ = 0, g = 0.05 meV (g = 0.039 meV), and T = 50 K, calculated in the LN (L = 9, red solid), NN (red dotted), and PA approach (black dashed line).The LN result is fitted with two complex exponentials [see Eq. ( 69)], to separate the long-time from the shorttime behavior, here shown as a relative error of the fit (green curve).The inset shows the initial dynamics.
Figures 5 and 6 show, respectively, the time dependence of the FWM signal |P FWM (t, 0)| and the corresponding spectrum | PFWM (ω)|, which is the Fourier transform of the former, PFWM (ω) = ∞ 0 P FWM (t, 0)e iωt dt.One can see a fast initial decay of the polarization (Fig. 5) followed by a multi-exponential (here, bi-exponential) behavior at later times.The former is well seen in the inset and the fit error (green line).The initial decay can be attributed to the rapid polaron formation as a result of the instantaneous excitation of the QD by a laser pulse.The polaron formation occurs on a timescale of τ IB = 3.25 ps and corresponds to a phonon broad band (BB) in the spectrum, represented by the Fourier transform of the complex absolute error of the fit (green line in Fig. 6).The Fourier transform of the fit corresponds to the zero-phonon line (ZPL) in the spectrum.The ZPL is described in this case by four complex parameters, which are illustrated in the inset.The fit error, i.e. the difference between the bi-exponential fit and the full calculation shown by the green line in Fig. 5 is in the 10 −4 − 10 −3 range for long times and gradually (exponentially) increases at earlier times, reflecting the non-Markovian dynamics of the polaron formation.
The fitted long-time behavior of the LN is captured well by the PA, although the decay rates are slightly underestimated, as can be seen from the bubble plot (the inset in Fig. 6), representing the complex frequencies and the amplitudes of both exponentials (the circles are centered at the complex frequencies, with the circle area proportional to the modulus of the amplitude).In fact, the phonon BB is not taken into account in the PA, and any phonon-related effects are limited to reduction of the exciton-cavity coupling and a change in the detuning, see Eq. (61).With this renormalization, the quantum dynamics in the PA is entirely determined by the JC model.The NN approximation is capable of capturing correctly the BB and also describes the quantum x-x FWM spectrum, i.e. the Fourier transform of the FWM polarization, calculated in the LN and PA approaches shown in Fig. 5.The full numerical LN result (red), is separated into ZPLs (the inset) and a phonon broad band (green line), corresponding, respectively, to the long-time and short-time parts of the signal.The upper (lower) part of the spectrum has a logarithmic (linear) scale.The inset shows the complex frequencies of the bi-exponential fit Eq. ( 69) (blue crosses) and PA (black dots) along with the corresponding amplitudes |Ai| given by the circle area and their phases (all shifted by π) color coded.
dynamics at longer times quite well, which is expected in this regime of small g [27,73].Being very efficient and straightforward, the NN approximation is, however, a semianalytic approach.The fully analytic PA is derived from the NN approximation in the long-time limit, so the PA and NN naturally agree at long times and differ only in that PA underestimates and NN overestimates the decay compared to the LN.
Let us now focus on the case in which the initial excitation and the measurement of the response are both done via the cavity mode.Figures 7 and 8 show, respectively, the time dependence and spectra for linear (upper panels) and FWM (lower panels) polarizations.In the linear response, the spectrum is dominated by two Lorentzian lines, corresponding to quantum transitions between the ground state and the first rung of the JC ladder, the same as in the x-x polarization, see Fig. 8(a).In the time domain, the linear c-c polarization also starts from one but does not show any initial pure dephasing, see Fig. 7(a).As a result, the BB is practically absent in the spectrum.The c-c FWM response shows a few notable differences compared to the x-x FWM and c-c linear polarizations.One difference arises from the fact that there are six transitions involved in the quantum dynamics, including four transitions between the first and the second rungs, compare the insets in Figs.8(a) and (b).The spectra in Figs.8(a) and (b) show, respectively, two and six Lorentzian lines, both in the fit and PA.As a result of destructive interference between the first-and secondrung transitions [13,16], the FWM spectrum is shrunk significantly (i.e. the spectral tails are suppressed) compared to the linear spectra.The other difference is that in the time domain, the FWM signal starts from zero, see Fig. 7(b).This is because the excitation of the system via the cavity mode does not immediately lead to an optical nonlinearity, because the cavity itself is linear (photons do not interact directly with each other).In other words, the excitation has to be converted from photon to exciton and back to photon, in order to produce a nonlinearity measured in the c-c FWM, with the rise time of the signal proportional to t 2 [16].
The absence of a visible phonon BB in the spectrum is a feature common to all c-c polarizations.Indeed, this is what we see in Fig. 8, and there is no apparent fast decay in the LN data at short times, see Fig. 7.This is because the polaron formation is adiabatically slow in the case of cavity excitation and small g.In fact, the polaron cloud forms around the QD or disappears on the timescale τ JC , and since τ JC τ IB the quick non-Markovian dynamics of the polaron formation is spread over the much slower JC evolution.As a result, there is a better overall agreement between LN, PA, and NN, both in the linear and FWM signals.Also, the fit error is reduced by an order of magnitude compared to the x-x polarization (Fig. 7), reaching 10 −5 for the FWM polarization.Still, there is a few times smaller (compared to x-x) leftover after the multi-exponential fit at short times and, as a consequence, a small BB is present in the spectra (green lines in Fig. 8, magnified by a factor of 10 3 and 10 5 for the linear and the FWM response, respectively).Surprisingly, the BB in the FWM corresponds to a much longer non-Markovian time compared to the linear and even the x-x FWM polarizations (10-20 ps versus 3 ps).This is due to the quantum nonlinearity of the JC system, leading to a quicker polariton dynamics when the second-rung transitions are involved.These non-Markovian features enhance dramatically for higher-order nonlinearities, involving higher rungs of the JC ladder, and for larger values of g.The latter is demonstrated and discussed in more depths in Sec.III B below.

B. Non-Markovian dynamics with increasing g
While the two approximate solutions, NN and PA, provide simple and convenient results, they show strong deviation from the correct behavior as the exciton-cavity coupling strength g and the temperature T are increased.For the linear polarization, increasing g and T are known to cause the breakdown of the PA and approaches based on the polaron master equation [27,28].Here, we present results for the c-c FWM and linear response, focusing on the role of the coupling strength g in the quantum dynamics.The temperature dependence is studied in Sec.III C below.
As the timescales τ JC and τ IB become comparable, the LN shows a slower convergence with L, forcing us to use a finer time grid, with an increased memory kernel (which contains 6 L+1 elements for the c-c FWM).In fact, with g increasing, the time step must be reduced also to resolve a faster JC system dynamics.To understand the effect of g on phonon induced dephasing, it is useful to first consider a simpler case of the linear polarization, where the kernel memory size is much smaller (contains 2 L+1 elements).In this case we can store more time steps in memory and use a finer grid.To make sure that the convergence was reached even in a regime of comparable timescales, L = 15 was used in Ref. [27].In this work, we use L = 9, which is also sufficient, with RMSD of 0.03-0.3%,as demonstrated in Appendix E. For the FWM response, we do not expect to see any significant change in convergence with L (see Appendix E) and all relevant features of the physical behavior can be captured with L = 9.
Results for the linear and FWM c-c polarizations are shown in Figs.9-12 for the exciton-cavity coupling strengths g = 0.3 meV and 0.8 meV.The separation of the time dependence into the long-time and short-time dynamics, using the multi-exponential fit Eq. ( 69), is attempted for these larger coupling strengths, and the results of such fits are also included in Figs.9-12 and compared with the multi-exponential behaviour predicted by the PA, now demonstrating a clear deviation from the latter.In the LN, the actual Rabi splittings of both rungs are larger than those predicted by the PA, but for g = 0.3 meV, the PA offers a reasonably good approximation.For g = 0.8 meV instead, a better estimate of the Rabi splittings is given by the bare JC model, i.e. without the polaron renormalization, but both approximations, with or without the renormalization fail to capture the full dynamics.In particular, the scaling of the splitting with the square root of rung number n is no longer observed in the LN: The average scaling of the Rabi splitting with √ n is reduced (for the second rung) by around 13% for g = 0.3 meV and 5% for g = 0.8 meV, as compared to only 0.4% reduction for g = 0.05 meV.One can infer from this comparison a phonon-induced modification of the JC energy level structure.
As demonstrated by Fig. 7 in Sec.III A, the presence of first to second rung transitions in the c-c FWM dynamics increases dramatically the timescale of the initial non-Markovian behavior as compared to the linear or x-x FWM polarization, involving only ground state to first rung transitions.This effect is getting stronger as g increases.Moreover, this non-Markovian timescale grows with g both in the linear and FWM polarizations, as it is clear from Figs. 9 and 11 (green lines).Consequently, the corresponding non-Markovian contribution to the spectrum is getting narrower and higher, reducing the ZPL wights, see Figs. 10 and 12.
For g = 0.3 meV, the linear and the FWM response spectra in Fig. 10 are very different, which is not seen as clearly in the time dependence in Fig. 9.However, it is clear from Fig. 9 that the initial deviation from a multiexponential behavior is more prominent in the FWM dynamics, same as in Fig. 7 but with longer non-Markovian times both for the linear and FWM polarizations.We note that the overall spectrum in Fig. 10 (b) appears as a triplet.A tendency of nonlinearities of the JC model to manifest as triplet structures in the presence of pure dephasing, has been previously reported in Ref [84], but in the photoluminiscence spectra of the system subject to continuous incoherent pumping.The spectra in Fig. 12 for g = 0.8 meV show an increased asymmetry (which  is discussed in Sec.III C) and higher deviation from the PA, while the temporal dynamics demonstrates a further increase of the non-Markovian times for both linear and FWM polarizations.
To understand the enhancement with g of the non-Markovian component in the linear and FWM c-c polarizations, let us recall that within the IB model, the exciton-phonon interaction results in a non-Markovian pure dephasing on the timescale τ IB , which is the memory time of the phonon bath and the time of the polaron formation around the QD in the absence of the cavity or for a small coupling strength g between the cavity and the QD.This pure dephasing manifests itself as a quick non-exponential drop of the excitonic polarization, see e.g. the inset in Fig. 5.For small g, i.e. in the limit τ JC τ IB , this quick non-Markovian dynamics of the polaron formation is spread over the much slower JC evolution, so that on any time interval ∆t ∼ τ IB , only a very little change of lattice deformations around the QD occurs, and the whole process of the polaron formation or disappearance (due to the JC Rabi rotations) can be treat as adiabatically slow.In other words, the non-Markovian part of the overall dynamics of the whole system, which would manifest itself as a non-exponential behaviour, is negligible, apart from a very short period of time at the beginning during which the systems reaches a "steady state" in this adiabatic process.This can be clearly seen in the quality of the multi-exponential fit (green curves in Fig. 7).As g increases and the timescales τ JC and τ IB become comparable, this period of time to reach the steady state is getting longer, as it is clear from Figs. 9 and 11 (green curves), and the quasi-periodic process of the polaron formation and disappearance can no longer be treated as adiabatic.Finally, in the c-c FWM polarization for g = 0.8 meV, this steady state cannot be reached at all, see Fig. 11(b).In fact, even though the quality of the multi-exponential fit looks reasonable for sufficiently long times (t > 50 ps), the fit is not unique in this case and does not bring much value to the understanding of the spectrum, see the green line and the inset in Fig. 12(b).In this regime, the quantum dynamics observed in the FWM and even in the linear polarization becomes essentially non-Markovian, with a much stronger effect in the FWM, due to a larger number of quantum transitions involved and owing to the quantum nonlinearity of the JC system.

C. Temperature effects and spectral asymmetry
Here we concentrate on the temperature effects in the FWM polarization.Figures 13 and 14 show, respectively, the zero-delay FWM c-c polarization and its spectrum, for two different temperatures (T = 4 K and T = 50 K), and for three different coupling strengths (g = 0.05 meV, 0.3 meV, and 0.8 meV), in each case comparing LN results with the PA.For small g, the LN shows good agreement with the PA, as previously discussed (see Sec. III A), and the agreement is better at low temperatures.Note a significant change with temperature of the doublet splitting in the spectrum in Fig. 14(a), corresponding to the change in the beating frequency in Fig. 13(a), in agreement with the PA, see Eqs. (61) and (62).In fact, the coupling constant g is renormalized (reduced) stronger at larger temperatures, as dictated by the temperature dependence of the Huag-Rhys factor S.
For an intermediate exciton-cavity coupling strength of g = 0.3 meV [panels (b) in Figs. 13 and 14], the impact of phonons on the FWM is much stronger.The excitoncavity dynamics occurs faster (in accordance with a larger Rabi splitting), and phonons do not have enough time to adapt to the changes in the QD-cavity system.In this regime, the behavior is more complex and in general requires a non-Markovian treatment, which is provided by the LN solution.For both temperatures, the spectrum for g = 0.3 meV is asymmetric due to phonon-assisted transitions between the polariton states [27].Surprisingly, for T = 4K, the PA, despite its fundamental sim- plicity, still gives a reasonable agreement with the LN result but fails to predict the spectral asymmetry.At T = 50K the spectrum shows some drastic changes in shape as compared to the PA: The frequencies of all six transitions and their linewidths are significantly modified by the presence of the phonon environment.
At even larger exciton-cavity coupling strength of g = 0.8 meV [panels (c) in Figs. 13 and 14] the deviation between the PA and LN approaches becomes dramatic.The exciton-cavity dynamics now occurs on a faster timescale than the exciton-phonon dynamics.Note that comparable timescales τ JC ≈ τ IB are observed in the linear polarization at g ≈ 0.6 meV (see Fig. 3 in Ref. [27]).The spectral asymmetry becomes more pronounced for both temperatures.For T = 5 K, the PA still agrees with the LN results surprisingly well, apart from the asymmetry.This can be understood from the fact that at low temperatures, there is a very little polaron renormalization of the coupling g, hence the transition energies match quite well in both approaches.This is not true at higher temperatures for large g, see Fig. 3 in Ref. [27] demonstrating a strong deviation of the actual Rabi splitting in the first rung of the JC ladder from that predicted by the PA.This effect becomes even more pronounced in the second rung involved in the FWM dynamics.

IV. CONCLUSIONS
We have presented an asymptotically exact solution for the FWM response of a QD embedded in a photonic cavity, taking into account the LA phonon contribution in full.We have implemented a real-time path integralbased approach, combining Trotter's decomposition and linked-cluster expansion.Our work presents an important generalization of the method recently developed for calculating the linear polarization [27].This generalization, applied here to the FWM response, clearly demonstrates how the method can be used for other elements of the density matrix and applied in principle to any observable.This semi-analytic approach unites the behavior of two exactly solvable models into a common framework, where they can influence each other by means of a memory kernel responsible for temporal correlations within the system.Our solution has a microscopic origin and is capable of addressing intrinsically quantum dynamics beyond perturbative regimes.We were thus able to explore the effect of phonons on linear and nonlinear QD-cavity dynamics with comparable exciton-cavity and excitonphonon timescales, as well as arbitrary temperatures.
The exactly solvable JC and IB models constituting the system are characterised by the timescales τ JC and τ IB , respectively.In the limit of fast phonon environment, τ JC τ IB , our microscopic approach allows us to develop and evaluate the validity of two useful approximations: the nearest-neighbor approach and the fully analytic polaron approximation following from it.For τ JC τ IB , both approximation describe the dynamics of the system well but are computationally much simpler, with the nearest-neighbor approach being able to capture even the non-Markovian dynamics at the initial times (t τ IB ).In the polaron approximation, the effect of phonons is reduced to the polaron shift of the exciton frequency and a renormalization of the exciton-cavity coupling strength by the Huang-Rhys factor, so that the quantum dynamics is described as a superposition of exponentials.Inspired by the polaron approximation, we have applied a complex-valued multi-exponential fit to the exact solution, that allowed us to precisely quantify the individual quantum transitions contributing to the FWM signal.Moreover, the multi-exponential fit separates the full quantum dynamics into a short-time non-Markovian and a long-time multi-exponential behavior.The latter represents optical spectra as superpositions of complex Lorentzian lines, corresponding to quantum transitions between phonon-dressed states of the JC ladder, in agreement with the polaron approximation.We have demonstrated this by focusing on the excitation of the optical nonlinearity via the cavity mode, in which case the second-rung states of JC ladder are directly excited and a more complex dynamics (compared to the linear polarization) arises.We have also studied other excitation and measurement channels corresponding to the third-order optical nonlinearity.
We have studied the regime of τ JC τ IB for the exciton-cavity coupling strength of g = 0.05 meV which is close to the one measured in Ref. [13].For this coupling strength, the calculated FWM spectra are qualitatively similar to those simulated in Ref. [13] for g = 0.035 meV without phonons.In this regime, phonons around the QD can quickly adapt to any changes caused by the Rabi rotation.For larger values of the coupling strength, g = 0.3 meV (close to the one measured in Ref. [75]) and g = 0.8 meV, the similarity of the timescales τ JC and τ IB opens up a possibility of phonon-assisted transitions between different polariton states, previously studied in the linear polarization [27].These transitions cause spectral asymmetry which we have also demonstrated in FWM spectra.
In the regime of comparable system and environment timescales, τ JC ∼ τ IB , we have found an increased deviation from the multi-exponential behavior, clearly demonstrating a non-Markovian character of the exciton-cavity FWM dynamics, not seen in the linear polarization.The phonon cloud around the QD is unable to adiabatically adapt to a varying optical state that results in non-Markovian effects seen in the quantum dynamics.We have found a significant deviation from the polaron approximation and multi-exponential fit, which becomes even more pronounced at elevated temperatures.
In a model of increased complexity, where linear and nonlinear effects, as well as environmental interactions are combined, there is a richer variety of interesting physical behavior that can be seen.While the dynamics of the full system exhibits characteristics of its individual parts, when they are combined the actual behavior is different that results in emerging phenomena, such as those demonstrate above.These phenomena can only be captured in a fully microscopic treatment, which is presented in this work.decay of the phonon autocorrelation function Eq. (B3).Note that in our numerical calculations, we use slightly different values of τ IB which are determined by the condition D(τ IB ) − D(t → ∞) = D(τ IB ) + iΩ p t + S ≈ 10 −4 .As for the above given timescales τ JC1 and τ JC2 , they express the exact values of the periods of the Rabi rotations in the first and second rungs of JC ladder at zero detuning and in the absence of phonons, taking into account the √ n increase of the Rabi splitting with the rung number n.
Figure 16 shows all three timescales as functions of g.For the linear polarization, the maximum phononinduced dephasing is observed in the regime of comparable JC1 and IB timescales, which is at g ≈ 0.6 meV [27].We have also observed similar behavior of the first rung transitions in the FWM.In Section III we have shown the FWM results when both the excitation and the measurement are done in the same mode.Here we consider the other possibilities, when the system is excited via the cavity, while the response is measured via the exciton channel (c-x FWM, see Fig. 17) and vice versa (x-c FWM, see Fig. 18).The former is similar to the c-c polarization shown in Figs.7(b) and 8(b), although the results are quantitatively different.In particular, a faster rise time of the signal is seen in Fig. 17, since the cavity excitation needs to be converted to the excitonic polarization in order to be observed in this nonlinearity channel.This takes a shorter time as compared to a further conversion of the excitation back to the cavity, which is required in the c-c FWM, see a discussion at the end of Sec.III A. Another feature is a larger relative phase difference between the "inner transitions" [16], compare the insets in Figs.8(b) and 17(b).Also, an increased contribution of the secondrung transitions is seen in Fig. 17(b).All of these features are well reproduced in the PA.The x-c FWM signal in Fig. 18(a) demonstrates an even shorter rise time, and its deviation for the bi-exponential fit is similar to the x-x FWM polarization, shown in Fig. 7.The spectrum also consists of only two transitions, but this time with very similar phases.
So far, we have considered the situation, where both excitation pulses Q (I) and Q (II) correspond to the same mode.In general, this is not always the case, and we give here examples of the xc-c and xc-x FWM polarizations, for which the first (second) pulsed excitation Q   which to a more pronounced BB in the spectrum, see the green lines in Fig. 20.Because of the stronger BB and the presence of all six transitions at the same time, below we explore this FWM channel further for different values of g and T , see Figs.21 and 22.
For g = 0.3 meV and T = 4 K, the xc-x FWM spectrum is wider, and the individual transitions are already well resolved as opposed to the c-c FWM, compare Figs.21(b) and 14(b).At T = 50 K there is still a considerable overlap of the transitions in the spectrum and only a small spectral asymmetry.For g = 0.8 meV (Fig. 22) there is a deviation from a regular periodic behavior in the time domain.The deviation from a multi-exponential behavior appears to be also very significant, so these results cannot be fitted at all.The transitions are spectrally better resolved, having less overlap and rather similar weights, both in the LN and PA.The LN spectrum has an interesting profile, with more uneven distribution of weights caused by phonon assisted transitions.
Another possibilities of excitation and measurement channels are cx-x and cx-c, produced by Q x .However, no signal is produced in these channel for zero delay between the pulses, since Q (II) x Q (I) c = 0, as can be seen from Eqs. ( 16) and (21).Any other excitation and measurement conditions can be expressed as linear combinations of those discussed here and in the main text.
becomes visible only on a log scale, when the signal is small.
We have experienced difficulties in obtaining FWM results with good precision using TEMPO.The number of combinations of different convergence parameters we could sample was limited by a long computational time.Therefore we only show a visual comparison of the FWM results with TEMPO in Fig. 25.As stated in the main text, it is sufficient to consider six elements of the density matrix in order to calculate the FWM response.We note that in the TEMPO results shown here, we have used the full exciton-cavity density matrix for the first two rungs of the JC ladder (36 elements), which may have resulted in a poorer convergence.

FIG. 2 .
FIG.2.An example of the time grid for the FWM polarization with non-zero delay in the NN approach with NI = 4 and NII = 3, showing the delay time (τ ) region (I), the observation time (t) region (II) and the mixed region (I-II).Only the self interaction (main diagonal) and the NN interaction (nearest off-diagonal elements) are included.The evolution along this grid consists of L-shaped regions (these are bounded by a dashed line), which correspond to an element of a matrix G. Grid squares are colored to distinguish different cumulant elements [see Eq. (59)].
FIG. 5.x-x FWM polarization at τ = 0, g = 0.05 meV (g = 0.039 meV), and T = 50 K, calculated in the LN (L = 9, red solid), NN (red dotted), and PA approach (black dashed line).The LN result is fitted with two complex exponentials [see Eq. (69)], to separate the long-time from the shorttime behavior, here shown as a relative error of the fit (green curve).The inset shows the initial dynamics.
FIG.6.x-xFWM spectrum, i.e. the Fourier transform of the FWM polarization, calculated in the LN and PA approaches shown in Fig.5.The full numerical LN result (red), is separated into ZPLs (the inset) and a phonon broad band (green line), corresponding, respectively, to the long-time and short-time parts of the signal.The upper (lower) part of the spectrum has a logarithmic (linear) scale.The inset shows the complex frequencies of the bi-exponential fit Eq. (69) (blue crosses) and PA (black dots) along with the corresponding amplitudes |Ai| given by the circle area and their phases (all shifted by π) color coded.

FIG. 7 .
FIG. 7. As Fig. 5 but for (a) linear and (b) FWM c-c polarizations.For the FWM polarization the multi-exponential fit Eq. (69) consists of six exponentials, corresponding to the six transitions between the rungs of the JC ladder involved in the quantum dynamics.

FIG. 8 .
FIG. 8.As Fig. 6 but for (a) linear and (b) FWM c-c polarizations.The inset in (b) shows six exponentials of the fit (blue crosses and dotted circles) and PA (black dots and dashed circles).The amplitudes (circle areas) of the two transitions in the linear response in the inset in (a) are multiplied by a factor of 0.415 to match the amplitudes of the corresponding transitions in the FWM.

FIG. 10 .
FIG. 10.As Fig.8but for g = 0.3 meV.The amplitudes (circle areas) of the two transitions in the linear response in the inset of (a) are multiplied by a factor of 0.95 to match the amplitudes of the corresponding transitions in the FWM.

FIG. 12 .
FIG.12.As Fig.10but for g = 0.8 meV and the two amplitudes in the inset of (a) are scaled by a factor of 0.993.

FIG. 14 .
FIG. 14.Zero-delay c-c FWM spectrum for (a) g = 0.05 meV, (b) g = 0.3 meV, and (c) g = 0.8 meV, the Fourier transform of the time dependent polarization in Fig.13, calculated in the LN approach (red solid lines) and PA (black dashed lines), for T = 4 K (thin) and T = 50 K (think lines).

FIG. 17 .FIG. 18 .
FIG. 17.(a) c-x FWM polarization at τ = 0, g = 0.05 meV, and T = 50 K, calculated in the LN (L = 9, red solid), NN (red dotted), and PA (black dashed line).The LN result is fitted with six complex exponentials [see Eq.(69)], to separate the long-time from the short-time behavior, here shown as a relative error of the fit (green curve).(b) Fourier transform of c-x FWM polarization, calculated in the LN and PA approaches and shown in (a).The full numeric LN result (red), is separated into ZPLs (the inset) and a phonon broad band (green line), corresponding, respectively, to the longtime and short-time parts of the signal.The inset shows the complex frequencies of the multi-exponential fit (blue crosses) and PA (black dots) along with the corresponding amplitudes |Ai| given by the circle area and their phases color coded.
c ) occurs in the exciton (cavity) mode and the measurement is done, respectively, in the cavity and exciton mode.Figures19 and 20show, respectively, the

FIG. 19 .
FIG. 19.As Fig. 17but for the xc-c FWM channel, in which the system is excited by Q (I) x and Q (II) c .