Extracting Kinetic Information from Short-Time Trajectories: Relaxation and Disorder of Lossy Cavity Polaritons

The emerging field of molecular cavity polaritons has stimulated a surge of experimental and theoretical activities and presents a unique opportunity to develop the many-body simulation methodology. This paper presents a numerical scheme for the extraction of key kinetic information of lossy cavity polaritons based on the transfer tensor method (TTM). Steady state, relaxation timescales and oscillatory phenomena can all be deduced directly from a set of transfer tensors without the need for long-time simulation. Moreover, we generalize TTM to disordered systems by sampling dynamical maps and achieve fast convergence to disordered-averaged dynamics using a small set of realizations. Together, these techniques provide a toolbox for characterizing the interplay of cavity loss, disorder, and cooperativity in polariton relaxation and allow us to predict unusual dependences on the initial excitation state, photon decay rate, strength of disorder, and the type of cavity models. Thus, we have demonstrated significant potential in the use of the TTM towards both the efficient computation of long-time polariton dynamics and the extraction of crucial kinetic information about polariton relaxation from a small set of short-time trajectories.


Introduction
Recent experiments of molecular systems in optical cavities have brought much excitement in chemical sciences, material engineering, and quantum physics, and have inspired many numerical simulations of the coherent light-matter interactions in cavity systems.[1,2] These studies often push the boundary of computational capacities because of the coherent nature of many-body interactions and the presence of various dissipative channels.[3][4][5][6][7][8][9][10][11][12][13][14][15][16] In this sense, cavity systems present a unique opportunity to develop and extend the many-body simulation methodology.Specifically, early studies have revealed the complex relaxation dynamics of cavity polaritons but also opened questions for further research: (i) the interplay of static disorder and various dissipation channels; (ii) characterization of the relaxation pattern and timescale; (iii) dependence on the initial excitation state; (iv) comparison of standard polariton models.In this paper, we aim to develop novel computational techniques based on the transfer tensor method (TTM) to address these challenging questions numerically.
First introduced in 2014, the TTM is a novel black-box method for extrapolating long-time dynamics from short-time trajectories for any dynamical system with a time-translationally-invariant correlation [17].As shown in Figure 1a, starting from dynamical maps computed via an input-output analysis, one can construct transfer tensors and obtain the complete information of a reduced quantum system, and can thus propagate its dynamics to arbitrarily long times at an error rate that is independent of the propagation time [18][19][20].In the field of open quantum dynamics, several methods have emerged recently to compute the longtime dynamics or to predict the steady-state and kinetic information [21][22][23][24][25][26][27].These predictions can be reliable if the calculations are converged, but the propagation to long times can be numerically expensive and can potentially accumulate errors.In comparison, the TTM takes advantage of short time trajectories that are computed with high accuracy and predicts not only the steady state but also relaxation rates.More importantly, the TTM is a black-box technique, which can be applied to any dissipative environment, not limited to Gaussian bosonic baths, and can be combined with any quantum or classical dynamics methods.In this sense, SMatPI [26,27] also adopts the concept of transfer tensors and can be regarded as the application of the TTM to QUAPI.
While the TTM is a powerful tool for reconstructing non-Markovian dynamics from short-time trajectories, novel technical developments are needed in its application to cavity polaritons and other manybody quantum systems.Specifically, we will explore how to extract kinetic information from TTM and how to generalize the TTM concept to disordered systems.The layout of the paper is described as follows: In Sec. 2, we begin with a brief introduction to three cavity Hamiltonians and review the basics of the TTM.Specifically, we will introduce the Pauli-Fierz (PF) Hamiltonian, the Dicke Model, i.e., the coupling of N Two-Level Systems (TLS) to a photonic cavity, and its rotating wave approximation, the Tavis-Cummings (TC) model [28].In these model Hamiltonians, the reduced density matrix has dimensions 2 N by 2 N , and each transfer tensor has dimensions 4 N by 4 N : the computation thus grows exponentially with the size of the system.To resolve this issue, we adopt the numerical technique to use the full identity matrix of the system as the initial condition to learn the short-time quantum dynamics in a single learning trajectory.
While the transfer tensors contain all the dynamic information, they have been used primarily as propagators.In Sec. 3, we present our approach to directly extract crucial kinetic information including the steady-state and relaxation rate from the transfer tensors without propagating the density matrix.Further analysis suggests a numerically robust approach to identify oscillatory modes and their decay rates.Beyond rates, we can also evaluate high-order moments using transfer tensors and this characterize the deviation from the single exponential decay.
Together, these techniques provide a tool box for characterizing the polariton relaxation profile and predict the dependences on the initial state, system size, and photon decay rate.
Though the TTM is conceptualized for a given Hamiltonian, its generalization to disordered systems is of broad interest but has not been explored.In Sec. 4, we explore a novel application of the TTM to disordered systems by incorporating the random distribution of system parameters into the dynamical maps.In this way, the disorder-averaged TTM (i.e., DA-TTM) is capable of predicting the averaged relaxation behavior of the disordered system without repeated density matrix propagation for an extended simulation time.In fact, the DA-TTM can even converge faster than directly averaging the density matrices over the realizations, since it can extrapolate averaged results from a relatively small sample size.We apply the approach to predict the average relaxation behavior of disordered cavity polaritons and analyze its dependence on the strength of disorder and on the symmetry of the initial state.
In Sec. 5, we compare the three standard cavity models: the Pauli-Fierz (PF) Hamiltonian, the Dicke Model, the Tavis-Cummings (TC) model [28].We show that, in the weak coupling limit, all models converge to similar results, whereas their behaviour diverge in the strong coupling limit.We show that lifting the resonance between cavity and atoms can restore the similarity between all three models.

Cavity Polariton Hamiltonians
The Dicke model is widely used in quantum optics for describing a variety of physical phenomena, such as superradiance [29].The standard Dicke model describes the dynamics of N TLSs coupled to a single-mode cavity [30] with Hamiltonian where ω c is the frequency of the cavity, ω a is the frequency of each of the identical TLS, and g is the coupling constant between the cavity and the TLS.The operator a is the annihilation operator of the photonic cavity and the σ j operator is the Pauli operator for the jth TLS.Note that while the number of cavity levels is physically infinite, to make calculations tractable, we truncate the cavity basis to n levels, so that the annihilation operator a is simply the rank n operator.
In the weak light-matter coupling regime, we often adopt the Tavis-Cummings (TC) model, which applies the rotating wave approximation (RWA) to the Dicke model and omits the counter-rotating terms in the light-matter interaction Hamiltonian [31].This approximation is valid in particular whenever ω c = ω a ≫ g.In this case, the Dicke Hamiltonian reduces to: In the strong light-matter coupling regime, we often use a generalization of the Dicke Model, known as the Pauli-Fierz (PF) Hamiltonian [32] where the term added to the Dicke Hamiltonian H D is a dipole self-energy (DSE) that accounts for cavity-mediated interactions between TLS of the form σ x j σ x k .It is worth noting that the transfer tensor method below is not affected by the choice of the cavity Hamiltonian, and can produce the dynamics of any model with identical computational effort.To simplify the analysis and validation of results below, we will present most of the numerical calculations using the TC model and will compare the three models in Sec. 5.
Beyond the coherent dynamics, we will consider the effect of a lossy cavity by means of a master equation where ρ f is the density matrix of the full system composed by the N TLS and the cavity and κ is the decay rate of the cavity.We often seek to estimate the relaxation rates of systems as a function of κ, fixing ω c = ω a = 0.As part of our simulations we also examine a variety of initial states: the singly and multi-excited manifolds of N TLS.For example, the singly excited initial state for 3 TLS may be represented as |↑↓↓⟩, and the fully excited as |↑↑↑⟩ or |N/2, N/2⟩.Concerning the initial state of the light, we restrict ourselves throughout to the vacuum state.

Transfer Tensor Method
The simulation of the full cavity photon and matter system ρ f features an unfavorable scaling that can be partly mitigated if the description is reduced to the matter subsystem alone.In general, the associated reduced density matrix ρ will feature non-Markovian and strong coupling effects that preclude the use of a simple time-local master equation.Further, different models of light-matter interactions define different non-Markovian features, which can characterize the light-matter entanglement.Here we propose the application of a general-purpose method to resolve this issue.First introduced in 2014, the transfer-tensor method (TTM) provides a computational process for extrapolating the complete dynamics of a quantum system from short time trajectories [17].As shown in Figure 1a, transfer tensors are extracted from the dynamical maps of the system E k associated with a uniform discretization of time t k = kδt, with δt the time step of the discretization.Dynamical maps are defined by their action on the initial density matrix of the system that is, the kth dynamical map propagates the initial density matrix to the kth density matrix according to the dynamics.By definition, E 0 is the identity superoperator acting on system operators.With the dynamical maps, one can then perform a linear transform to compute the transfer tensors via: with T 1 = E 1 .Finally, the transfer tensors are used to propagate the dynamics of the system, following the relation: As pointed out in [17], the success of this approach relies on the time-translational invariance of the underlying dynamics.This is guaranteed when (i) the total Hamiltonian is time independent, (ii) the initial total state is a product state, and (iii) the initial environment state is stationary.In this paper, one will find that (i)-(iii) all hold for our study of the cavity models.It has been shown that condition (ii) may be relaxed [33].Further, in the limit δt → 0, the TTM can be directly connected to the Nakajima-Zwanzig quantum master equation, which allows for generic formulation of open quantum processes [34,35]: In this sense, the transfer tensors are viewed as containing both the Liouvillian, L 0 , and the memory kernel, K(t).That is, the first transfer tensor T 1 encodes precisely the Liouvillian via the relation T 1 = (1 − iL 0 δt).Similarly, the other transfer tensors T i for i ≥ 2 encode the convolutional memory kernel via the simple relation T k = K(t k )δt 2 .In total, this can be summarized with the following relation: Although the above expression brings the TTM into the standard quantum master equation formalism of open quantum systems, the TTM is a self-contained conceptual framework and a general computational strategy.In fact, the 3 conditions underlying the TTM are sufficient to deduce the Nakajima-Zwanzig equation without using the projection operator formalism.
Theoretically, the TTM has many applications, namely as a dynamic propagator whose accuracy is not determined by propagation time, but rather only by its learning time [20,22,36].Analysis has shown TTM is especially promising in scenarios where the propagation time is much longer than the correlation time of the environment [37].

Demonstration of TTM
In this section, we demonstrate elementary computational results from applying the TTM directly to the TC model, i.e., the RWA of the Dicke Model.For this demonstration, as well as future computations, we take advantage of QuTiP, a python package for simulating quantum systems [38].
To apply the TTM, one must begin by computing the requisite dynamical maps from short-time dynamics for the N TLS.This may be done by considering the time-local Liouvillian eq. ( 4) acting on Figure 2: Simulation of the 4-TLS TC Model with a 5-level cavity.The system is parametrized g = 10κ/2, where κ is the decay rate of the cavity.Notice that the resonant TC model is simply Throughout the manuscript, κ or related quantities are taken as the inverse unit of time.The dynamics are computed for two different learning times, κt = 1.5 and κt = 3.5; the latter is found sufficient to accurately capture the reduced dynamics of the system until arbitrary time length.the density matrix ρ f (t) of the full system including the cavity In practice, we will need to restrict ourselves to a truncation of the cavity to the lowest n Fock levels.
For the TC model it is sufficient to consider n = N + 1.For other models, the particular value of n will depend on the strength of the cavity-matter coupling and in practice needs to be converged for each case.
The full density matrix ρ f may be reduced to the density matrix of the N TLS ρ by the action of an appropriate projection superoperator where T r C is the partial trace on the degrees of freedom of the cavity and ρ C is the initial state of the cavity.By replacing the initial condition for the solution of eq. ( 10) with the projected identity superoperator of the full system and propagating it to times t k , PI f (t k ), the necessary dynamical maps are simply from which the transfer tensors may be derived as usual (eq.6).
To demonstrate the efficacy of this method, we apply the TTM to a 4-TLS and 5-level cavity TC Model in Figure 2. Unless otherwise stated, for this and subsequent calculations, we set the cavity and atom frequencies, ω c and ω a to be on resonance, and set the coupling constant, g, to be 10κ/ √ N , where N is the number of TLS.The only decay channel allowed is via the cavity annihilation operator, with rate κ.Throughout the manuscript, we take κ as a natural unit to define the inverse time and gauge the light-matter coupling.This facilitates the identification of coupling strength regimes.Additionally, we benchmark the TTM results with the exact results as computed by direct propagation of master equation (4).As shown in the figure, the TTM is very effective at propagating the dynamics of the TC model to times significantly longer than the learning time, although the choice of learning time must be reasonable for the TTM to accurately capture the dynamics.In practice, the connection between the transfer tensors and the memory kernel indicates that the learning time must be of the order of the correlation time of the bath.

Information Extraction via Transfer Tensors
As illustrated in Figure 1b, our first main result is to demonstrate how to use the transfer tensors to directly predict key kinetic information without propagation.There exist two equivalent paths to demonstrate this: 1. direct exploration of the TTM propagation rule eq. ( 7) and the unilateral z transform of the transfer tensors as defined by with K the index of the last transfer tensor considered, 2. or the Laplace transform of the generator of the Nakajima-Zwanzig equation ( 8) where f (s) = ∞ 0 e −st f (t)dt denotes the Laplace transform of the function f .In the limit of an infinitesimal discretization and by substituting z → e sδt , the Laplace transform of the generator is recovered from the transfer tensors which is equivalent to Eq. ( 9).We show in this section that either perspective allows for the deduction of kinetic information without requiring any propagation of the system.Beyond this connection, we show in Section 3.3 that further information can be extracted from T [z] as a result of the finite time step δt.This information is not readily available from L (s) alone.

Steady State
First, we demonstrate the ability to compute the final, infinite-time state (i.e.steady state) of the nonequilibrium dynamics directly from the transfer tensors.We establish the approach first in the continuous Laplace space and then in the discretized transfer tensor form, and finally present an application to lossy cavities.

Continuous Laplace Transform
To begin, via the Laplace transformation of the quantum master equation, the density matrix is solved formally as: Then, the steady state can be directly computed via the final value theorem where ρ(∞) = ρ ss .This formalism can be understood as the extraction of the overlap between the null subspace of the generator L(0) and the initial state ρ(0).

Discrete Transfer Tensors
An equivalent representation of the Laplace formalism is provided by the unilateral z transform of the TTM propagation in eq. ( 7) as where we used T 0 = 0.The final value theorem takes the form ρ ss = lim z→1 (z − 1) ρ [z], so This result is consistent with the fact that the steady state density matrix must not be affected by propagation with transfer tensors by definition which implies (1 − T [1])ρ ss = 0. Therefore, ρ ss is an overlap between the initial state and the null space of 1 − k T k .Alternatively, in the language of complex systems, it can also be interpreted as a fixed point of the transformation represented by the sum of all transfer tensors.Under the partially excited states, since the TC Hamiltonian preserves symmetry, the system is prevented from fully decaying, instead relaxing to a state that remains partially excited.

Application
As shown in Figure 3, we demonstrate the ability of the method to predict the steady states of both the fully excited and partially excited initial states of the TC model.Since this method requires only a handful of matrix arithmetic steps as opposed to fully simulating the system, it can efficiently compute the long-time dynamics of the system with only the short-time trajectories, while remaining agnostic to the "true" dynamics of the system.The steady-state of a lossy cavity depends on the symmetry of the initial state: The fully-excited initial state relaxes to the ground state (i.e., ⟨σ z ⟩ = −0.5),whereas the singly excited initial state or partly excited initial state can remain partly trapped in an excited state (i.e., ⟨σ z ⟩ ̸ = −0.5).

Relaxation Matrix
In this section we generalize the concept of relaxation rate to that of the relaxation matrix where ∆ρ (t) = ρ (t) − ρ ss is the deviation of the density matrix from its steady state.For systems that reach the steady state on a finite time-scale, this matrix is bounded.By definition, the relaxation matrix is Hermitian and in general not positive, and depends on the initial condition of the system.
The information contained in the relaxation matrix τ becomes especially intuitive when evaluated for the operator of interest ô.Specifically, we define which is directly related to the relaxation timescale τ (eq.22) via normalization, i.e., The relaxation matrix τ contains 2 N real eigenvalues, describing multiple relaxation timescales of the system.We may normalize by the initial condition where the symbol − implies the pseudo inverse such that just the operator minus its null space is inverted.The normalized relaxation matrix τ ′ contains a number of non-zero eigenvalues which describe the effective lifetimes of the different decay modes in the system.In particular, one of the eigenvalues correspond to the depletion timescale from the initial state.There exist as many other non-zero eigenvalues as the dimensionality of the null-space that the steady state ρ ss has overlap with.These eigenvalues correspond to the timescales needed to reach the steady-state.The explicit calculation and spectral analysis of the relaxation matrix will be left for a future study.
Based on TTM, we can obtain τ from direct analysis of the transfer tensors, thus avoiding a lengthy and demanding numerical propagation of the density matrix of a system and the subsequent integration.Below, we will briefly describe the procedure in both the continuous and discrete formalisms.

Continuous Laplace Transform
In analogy to the spectral method presented in [39,40], the eigenvalues of L(s) as a function of s describe all possible decay (real part) and oscillatory (imaginary part) behaviors of the system.Further, they contain all the non-Markovian effects of the bath.Nevertheless, for a system of N TLS, the analysis of 2 2N of these functions of s becomes convoluted.
Alternatively, the relaxation matrix τ represents a much more compact object that is directly related to the Laplace transform of the generator by where again the symbol − implies the pseudo inverse.

Discrete Transfer Tensors
Numerically, this calculation can be performed by means of the z-transform of the transfer tensors and the limit of z → 1 simply involves the sum of all transfer tensors k T k The transfer tensors are computed based on a time step δt, and the relaxation matrix τ converges as δt → 0.

Application
To demonstrate the effectiveness of this method, we consider an initially fully excited N = 2 TC model with light-matter coupling of g = 10κ/ √ 2 and compute the relaxation timescale τ of the σ z observable of either TLS.Since each TLS has initial spin of 1/2 and steady state spin of −1/2, the exponential decay fit takes the form of e −t/τ − 1/2.In Figure 4a, we plot this exponential decay fit, and find that our method provides a good estimate of the decay rate of the system purely from the transfer tensors.We now proceed to analyze the relaxation timescale associated with the dynamics presented in Figure 3.The operator associated with the relaxation measurement ô is the projector into the fully excited state |N/2, N/2⟩ ⟨N/2, N/2|. Figure 5 plots the inverse of relaxation timescale (i.e. the relaxation rate) as a function of the cavity decay rate κ of the fully excited initial state |N/2, N/2⟩.Three different system sizes of N : 2, 3 and 4 are considered.For sufficiently small κ, the rate grows approximately linearly (see Figure 5a).The rate of growth increases with N and, for N = 2, it coincides with the prediction of perturbation theory of 2κ/3.For larger values of κ, the rate does not follow a trend of proportionality (see Figure 5b).On the contrary, large cavity decay rates suppress the transfer of excitations from the atoms to the cavity and reduces the overall effectiveness of dissipation, to the point where it becomes inversely proportional to κ.Thus, the relaxation rate exhibits a turnover as a function of the cavity decay rate.These correspond to three relaxation regimes illustrated in Figure 4b.

Oscillatory Relaxation
In this section we show that, beyond the steady state and relaxation timescales, oscillatory information may be extracted directly from the transfer tensors or the relaxation matrix τ .By analizing their behaviour as a function of δt it is possible to make a Fourier transform analysis of the dynamics without actual propagation.We begin by analyzing a case study with a priory knowledge of oscillatory behaviour in Section 3.3.1.We continue in Section 3.3.2by analytically demonstrating the relationship between the Fourier transform of the dynamics and the analysis of τ as a function of δt.Finally, in Section 3.3.3we show that an enhanced analysis of oscillatory behaviour can be achieved by this method without prior knowledge of the dynamics.

A Case Study
The estimation of relaxation timescale τ can further be used to detect oscillatory modes in the system dynamics.This becomes especially apparent in the relaxation dynamics of a single excitation initial state.In this case we monitor the operator ô = σ z , the z Pauli operator (i.e., population difference), of the initially-excited TLS.As shown in Figure 3b, its dynamics (i.e., population evolution) is oscillatory around the steady state value.Positive and negative parts of the dynamics cancel each other, so the overall value of the integral τ becomes much smaller than what would be obtained from the decaying envelope and thus does not constitute an appropriate estimate of the decay rate (see Figure 6a).This   can be fixed by adjusting the time-step δt of extraction of the transfer tensors.By tuning δt to the period of the oscillations, a resonance effect takes place by which the decaying envelope is reproduced and the true relaxation timescale is captured (see Figure 6b).In general, an estimate τ can be computed as a function of δt.As δt coincides with the period of oscillations of the system or its multiples, larger estimate τ (δt) is observed.This is shown in Figure 7, where a resonance is observed once δt coincides with a multiple of the period of the oscillations in Figure 6a.An increase of δt involves an error of discretization in the calculation of the time integral, which is proportional to δt/2 and can be easily subtracted.

Resonance Analysis
The resonance effect as a function of δt is tightly connected to the Fourier transform of ∆ ⟨ô⟩ (t) = ⟨ô⟩ (t) − ⟨ô⟩ (∞).Let us explore this connection with an example of an arbitrary observable ⟨ô⟩ featuring an oscillatory decay described by Using transfer tensors with timestep δt to estimate τ , we obtain which is a function with peaks at multiples of the value δt = 2π/ω.Thus, the calculation of τ (δt) offers insight into the Fourier transform of the dynamics.As δt increases, an error proportional to δt/2 that is associated with the discretization of the integral calculation becomes apparent.Parameters: The effect of the discretization introduced by the TTM can be formally elucidated by means of a Dirac comb Ξ δt = k δ (t − kδt), so that By the properties of the Dirac comb under Fourier transformation, we have where F [f ] (ω) is the Fourier transform of f (t).Therefore, the estimate τ (δt) provides a sampling of the Fourier transform of the deviation ∆ ⟨ô⟩ at the frequencies ω k = 2πk/δt.Every time one of these frequencies matches a peak of the Fourier transform, it appears as a peak in τ (δt).In Figure 7, a single peak in the Fourier transform F [δ ⟨ô⟩] (ω) at the value ω = √ 2g produces a repetition of the peak at values δt = √ 2πk/g.The first peak may be extracted for several values of κ as shown in Figure 8a.The value of τ at the peak (corrected by the error δt/2) allows us to evaluate the change of the relaxation rate as a function of κ, which for the case N = 2 can be analytically proven to follow the linear dependence of κ/4.This is reproduced in Figure 8b.

Multiple Resonance
The proposed approach is also useful even when the system is initialized in the fully excited state |N/2, N/2⟩, which does not oscillate around the steady state but may exhibit an oscillatory decay for small enough κ.Let us first evaluate the probability of remaining in the initial state for N = 2, whose dynamics may be analytically solved in the perturbative limit and shown to oscillate at frequency √ 6g.By extracting the estimate of τ as a function of δt, we show in Figure 9a that two distinct resonances are detected, corresponding precisely to ω 1 = √ 6g and ω 2 = 2ω 1 .By inspecting the time dependence of the system, Figure 9a, it becomes clear that ω 1 matches exactly the oscillation frequency of the dynamics and hence maximizes the estimate τ .The double frequency ω 2 matches half the period of oscillation, so both the maxima and minima of the oscillations are sampled.Although the corresponding estimate for τ (δt) cannot be the optimal one, it still provides a higher value than other choices of δt.
For the same initial state, the expected value of σ z features a more complicated behavior as shown in Figure 4a.In particular, a beating is apparent between two close frequencies.This pattern shows up in the form of additional resonances in the estimate of the relaxation timescale τ , Figure 10a.Beyond the resonances already identified in Figure 9a, two more resonances appear that correspond to the beating signal in the dynamics, Figure 10b.
In conclusion, the relaxation timescale τ as a function of the time step δt features as a versatile tool to extract relevant oscillatory behavior of the dynamics of open quantum systems.

Moments and Poisson Indicator
It is known that the Laplace transform of a distribution function is also the moment generating function of time; as such, we can directly calculate all the moments of relaxation via the same formalism.Namely, we have: Thus, we can compute the first moment via: By evaluating this expression in the limit where s approaches 0, one obtains an estimate of E[t].While mathematically accurate, this method is numerically unstable, and small variations in either learning time (i.e., the size of the memory kernel) or computation method can lead to varying results.However, in order to look at decay dynamics, one must compute the moments of the density matrix relative to the steady state density matrix, i.e., ∞ 0 [t n ρ(t) − ρ(∞)]dt.From the derivations in Section 3.2.2,we use the relaxation matrix τ so that the first moment can be computed via: where we use the pseudoinverse.Similarly, we can compute the second moment via: Given an operator ô of interest, we may project the moment by and With high order moments, we can then characterize the deviation from the exponential decay using the Poisson indicator.The detailed numerical calculation will be left for a future study.Since all terms are given entirely by the information in the transfer tensors, we can evaluate this expression for any given short-time simulation from which we can extract transfer tensor information.Thus, given sufficient learning time, this presents a method to tractably compute moments of the dynamics of the system without requiring a full simulation of the system, instead of requiring finite matrix multiplications.With high order moments, we can then characterize the deviation from the exponential decay using the Poisson indicator.The detailed numerical calculation will be left for a future study 4 Disorder-averaged (DA) TTM Disordered systems require extensive sampling of initial conditions in numerical simulations.For example, we model disordered cavity systems by performing simulations on random realizations of the parameters (e.g., κ, g, ω c , ω a ), drawn from a given probability distribution.
As illustrated in Figure 1c,, we generalize TTM to disordered systems to uncover the disorder-averaged dynamics e.g., the effective dynamics by averaging over the random distribution.
Specifically, we consider the TC model with Hamiltonian given by Note that here, as opposed to eq. ( 2), each TLS may have its individual frequency ω j ; therefore, it is not possible to remove ω j uniformly via the rotating wave approximation.To introduce disorder into the system, we fix g = 10κ/ √ N (here, N = 2) and ω c = 50κ, but draw ω j /κ ∈ U (40, 50).Below, we use the TTM to extract the disorder-averaged dynamics of the system.
To do so, we explore two techniques: (i) We average the transfer tensors computed for each realization, over all realizations, i.e., where M is the number of samples, Tk is the kth transfer tensor used for propagation of the average dynamics, and T i k is the kth transfer tensor computed from the ith realization.(ii) We average the dynamical maps computed via each individual realization, i.e.
where Ēk is the average kth dynamical map used for computing the transfer tensor, and E i k is the kth dynamical map computed from the ith realization.Below, we compare the result of these two applications of the TTM to the disordered averaged density matrices, which we take to represent the average behavior of the TLS system.
Figure 11 summarizes the results of this comparison.As shown in all three plots, averaging the dynamical maps gives remarkably good agreement with the averaged density matrices, whereas averaging the transfer tensors exhibits noticeable deviations.In addition, the steady state predicted from the transfer tensors computed via the averaged dynamical maps recover the true, fully-decayed steady state of the system, as expected.Thus, applying the TTM to disordered systems via averaging the dynamical maps computed from each realization of disorder allows for successful reproduction of disorder-averaged behavior for the system.This approach, named DA-TTM', is physically justified, since averaging the dynamical map Ēk is equivalent to averaging the density matrix of the system at time t k over disorder, ρ(t k ), when applied to initial state ρ(0) common to all realizations.This is not the case for the averaged transfer tensors Tk .
The three plots in Figure 11 are for different initial conditions: (a) the fully-excited state, (b) a singly excited state, and (c) the coherent superposition state.In case (a), the initial state and subsequent dynamics is highly symmetric, thus the agreement is nearly perfect.Notably, in the asymmetric initial condition of case (b), there are significant fluctuations within the average of the 50 density matrices, which can be interpreted as originating from the increased variance in asymmetric initialization.In comparison, the DA-TTM yields better convergence, and thus allows us to deduce the averaged behavior of the system using fewer realizations.
Finally, we compute the largest eigenvalue of the normalized relaxation matrix τ ′ eq. ( 26) to extract the decay rate of our systems as a function of the disorder.For the TC model in the fully-excited initial state with disorder given by ω j ∈ 45κ ± δ, where δ is the disorder, we vary δ and compute the decay rates from the transfer tensors extracted via averaging 50 realizations.The result of this numerical experiment is shown in Figure 12, where the decay rate is plotted as a function of the logarithm of the disorder.For small disorders, the decay rate of the system is slightly reduced, as one may expect-the introduction of small disorder into the cavity frequency is only sufficient to break symmetry.However, for higher levels of disorder, the rate of decay decreases rapidly.Compared with our recent analysis of disordered cavities in the thermodynamic limit of the single excitation manifold, [41,42] the difference may arise from the fully excited initial state adopted in this calculation.
In short, we have shown that the DA-TTM is able to successfully recover the average behavior of a disordered system through the averaging of dynamical maps, which allow for the computation of transfer tensors that represent the disorder-averaged dynamics.Additionally, applying our technique for extracting the decay rate from the transfer tensors shows that increasing levels of disorder suppress the relaxation of the cavity system.Thus, the DA-TTM performs remarkably well for predicting the average behavior of disordered systems, especially for the fast convergence to a single, smooth solution as well as for the deduction of kinetic information about the underlying dynamics.

Comparison of Cavity Models
We now apply the tools presented in Section 3 to a comparison of the Dicke model ( 1), the TC model (2), and the Pauli-Fierz (PF) model (3).In this section we first explore the resonant cavity regime,   2 .The decay rate is plotted against log disorder for δ ∈ [0.001, 25].As shown, for small magnitudes of disorder there is little impact on the decay rate, but for larger disorders the decay rate quickly decreases.i.e. ω a = ω c = ω.While we can take ω = 0 in the rotating frame of the TC model, this is not the case in the Dicke or PF models due to the counter-rotating terms and, in addition for the PF model, the self energy term.
In Figure 13a we show the relaxation rate τ −1 for operator σ z as a function of the cavity decay rate κ for an initially fully excited state.We consider N = 2 and two values of ω.The turnover presented in Section 3 is reproduced in all cases, so both the overdamped and the underdamped limits are explored here.
Let us first discuss the case of ω = 100κ 0 .This corresponds to the regime of a weak cavity-atom coupling ω > g, since g = 10κ 0 / √ 2. In this regime, the counter-rotating terms in both the Dicke and PF models become negligible by virtue of the RWA.Their residual effect, as exposed by the curve corresponding to the Dicke model, is to reduce the overall efficiency of the transfer of excitations from the TLSs and into the cavity.This reduction becomes more critical as the cavity dissipates faster, so the difference between TC and Dicke/PF models becomes apparent only for sufficiently large κ (i.e., the overdamped regime).The role of the self energy term in the PF model is also negligible, but it has the ability to partially restore the excitation transfer efficiency.As shown, the relaxation rate of the PF model is slightly larger than that of the Dicke model.
The case of ω = 10κ 0 corresponds to the strong coupling regime where ω ≃ g.In this limit, neither the counter-rotating terms nor the self-energy term are negligible and they strongly affect the relaxation dynamics in the overdamped limit.In the Dicke limit, a polariton forms between the TLS and the cavity, such that the cavity is displaced from its equilibrium position that is proportional to the spin state of the TLS.[43] This new state relaxes more slowly towards the ground state.In the PF case, the effect of the self energy term is that of shifting and mixing the total spin states of the TLS.Resonance between the cavity and the TLS system is lost, and relaxation rate becomes even further suppressed than in the Dicke model case.
Polariton resonance can be restored in the PF model by increasing the frequency of the cavity mode.This is shown in Figure 13 (b), where we keep ω a = 10κ 0 , but increase the cavity frequency to ω c = 15κ 0 .The PF model increases its relaxation efficiency to a level similar to that of the Dicke model, while the TC model suffers a suppression of its ability to dissipate energy due to the lack of resonance between atoms and cavity.Their qualitative behavior in this limit is very similar, which informs the interpretation that the strong-coupling Dicke model must share traits with an off-resonant TC model and a resonant PF polariton.Our observation is consistent with the previous analysis of the Pauli-Fierz model, which displays the change of refractive index in the presence of matter and polarization fluctuations [44][45][46].
Recent studies have explored the roles of the dipole self-energy and the rotating wave approximation in these model Hamiltonians.[4,47,48] Our calculation of the relaxation rate complements these studies, which are mostly based on eigen-solutions and perturbative analysis.

Conclusion
In summary, we have developed a novel approach to directly extract kinetic information from the transfer tensors without requiring long time propagation and apply the approach to analyze the relaxation process of disordered and lossy cavity polaritons.Technically, we have exploited several aspects of the TTM: 1.The full identity matrix is employed as the initial condition to learn the short-time dynamics in a single learning trajectory.
2. The null space of the sum of all transfer tensors minus the identity determines the steady state of a given propagation.
3. The transfer tensors can be transformed into a relaxation matrix, which contains information about decay rates and oscillatory dynamics.
4. The TTM is also viable for extracting the average behavior of disordered systems via sampling the dynamical maps over realizations.
In particular, we have demonstrated that the information contained in the transfer tensors combined with the initial state of the system is sufficient to compute its corresponding steady state and decay rates.We first applied this technique to a cavity model for a variety of system sizes and initial states, finding that the TTM can accurately predict the long-term equilibrium regardless of the system parameters.Then, we constructed the relaxation matrix from the transfer tensors, which contains the information about the system's relaxation towards its steady state.The projection of the relaxation matrix to a particular measurement defines the relaxation timescale (i.e., the decay rate) and its high order moments, and the tuning of the relaxation timescale at variable time steps characterizes the oscillatory behavior in the relaxation process.
Equally important is the successful generalization of the TTM to disordered systems.Specifically, the DA-TTM can accurately reproduce disorder-averaged phenomena via averaging the dynamical maps over realizations of static disorder, thus allowing us to examine the effects of disorder on relaxation kinetics.Further, the DA-TTM can achieve faster convergence than the direct average of the density matrices, since it can extrapolate averaged results from a relatively small sample size.
The application of these novel numerical techniques to polariton relaxation in lossy cavities reveals the rich interplay of disorder, dissipation, and cooperativity in light-matter interactions.Specifically, we have characterized the complex dependence of relaxation kinetics on the initial excitation state, system size, cavity decay rate, strength of disorder, and the type of cavity models.
1.The steady-state of cavity polaritons depends not only on the equation of motion but also on the symmetry of the initial excitation state: The symmetric fully-excited initial state relaxes to the ground state, whereas asymmetric partly-excited initial state will be trapped in an intermediate state without complete relaxation to the ground state.
2. For the Tavis-Cummings model, the relaxation rate is linearly correlated with the cavity decay rate of the system in the weak cavity loss limit, and the coefficient of the linear dependence depends on the number of TLSs and the initial excitation state.However, in the strong cavity loss limit, the relaxation rate is inversely proportional to the cavity decay rate and is independent of the system size.
3. The non-monotonic dependence on the photon decay rate defines a turnover which corresponds to the most efficient relaxation.The turnover also signals a transition in the relaxation profile from the underdamped oscillations to overdamped decay.
4. While most studies have been carried out for disordered polaritons in the single excitation manifold, our DA-TTM explores the relaxation of the highly-excited initial state in a disordered cavity and reveals a strong dependence on the initial state.In general, the relaxation rate slightly decreases in the weak disorder regime, as disorder in the cavity frequency is only sufficient to break symmetry, but drops rapidly in the strong disorder regime.
5. A comparison of standard polariton models, including the Pauli-Fierz (PF) Hamiltonian, the Dicke model, and the Tavis-Cummings (TC) model, reveals a universal turnover as a function of the photon decay rate.Though reasonable agreement is achieved in the perturbative and on-resonance regime, significant differences are observed in the strong light-matter coupling regime which are lifted in the off-resonance case.
There are, however, multiple routes for significant further progress.For one, the TTM has the potential to be combined with any numerical methods for simulating realistic molecular systems in cavities, such as ab initial modeling, mixed quantum-classical dynamics, non-adiabatic quantum dynamics.[3,[8][9][10][11][12]14] Additionally, while we have shown the immediate application towards polariton relaxation, our toolbox can be applied generally to dissipative quantum processes including quantum transport and quantum thermodynamics,with or without couplings to cavity photons.
In terms of methodology, the TTM can be further developed for more efficient analysis of cavity systems.First, while the decay lifetime provides an estimator of the relaxation process of the system, it remains to be seen if the time-dependent relaxation matrix can be fully exploited as it contains all the dynamic information.Secondly, by combining symmetry reduction and information extraction from the transfer tensors, we can produce significant computational speed-up.Thirdly, one may consider alternative dimensionality reductions beyond taking advantage of symmetry: e.g., the contraction of the density matrices via an operator, followed by TTM propagation, allowing for prediction of a desirable property with significantly less computation.

Figure 1 :
Figure 1: Flow charts of TTM and its generalizations: (a) the original TTM, where {ρ(t k )} represents a set of short-time density matrices in the learning period, {E(t k )} and {T (t k )} are the corresponding dynamical maps and transfer tensors, respectively, and {ρ(t long )} represents a set of long-time density matrices predicted by TTM; (b) TTM-based kinetic analysis, where the last step exploits the transfer tensor to extract kinetic information such as the steady-state ρ ss , decay lifetime τ , resonance frequency ω res and its decay time τ res ; (c) DA-TTM, where a disordered ensemble of short-time density matrices is used to generate the disorder-averaged dynamical maps and corresponding transfer tensors, which then directly yield the disorder-averaged long-time dynamics or kinetics.

Figure 3 :
Figure 3: Steady States inferred via the simple matrix calculation from Section 3.1.All diagrams are constructed for the TC Hamiltonian H T C = ℏg(σ + a + σ − a † ) with g = 10κ.Depicted are steady state (red) and initial dynamics (blue) for: a) Fully excited N=2 initial state, b) Singly excited N=2, c) Fully excited N=4, d) Triply excited N=4.Under the partially excited states, since the TC Hamiltonian preserves symmetry, the system is prevented from fully decaying, instead relaxing to a state that remains partially excited.

Figure 5 :
Figure 5: Relaxation rate τ −1 of the probability of the fully excited state as a function of κ for N = 2, 3 and 4 for the fully excited initial state and g = 10κ 0 / √ N .(a) For smaller values of κ, the growth is linear.The linear growth for N = 2 follows the perturbative scaling 2κ/3.(b) For larger values of κ, the relaxation rate is upper bounded and eventually becomes inversely proportional to κ.

Figure 6 :
Figure 6: Dynamics of the initially excited TLS ⟨σ z ⟩ (t) for N = 2 (blue).Time is expressed in units of κ −1 .(a) The decay rate estimate τ in the function [3 exp(−t/τ ) + 1]/4 (orange) misses the decay of the envelope.(b) If the timestep δt is adjusted to match the oscillations observed in the left (δt = √ 2π/g), just the decaying envelope is observed by an oversampling effect.The decay rate estimate τ now approximately matches the decay of the envelope.Parameters g = 10κ/ √ 2.

Figure 7 :
Figure 7: Estimate τ for σ z as a function of δt for N = 2 (blue).Clear resonances occur when the timestep δt hits the period of the oscillations observed in Figure 6a, i.e. δtκ = π/5, and its multiples.As δt increases, an error proportional to δt/2 that is associated with the discretization of the integral calculation becomes apparent.Parameters: g = 10κ/ √ 2.

Figure 12 :
Figure 12: Results from simulating the decay rate as a function of disorder for the 2-TLS Tavis-Cummings Model in the fully excited initial state with H = hg(σ + a + σ − a † ) + hω a N j=1 σ x j + hω c a † a and disorder given by ω a ∈ U (45 − δ, 45 + δ), ω c = 50 and g = 10

Figure 13 :
Figure 13: Relaxation rate of σ z as a function of cavity decay rate κ for three models: TC (equation 2) in blue, the Dicke model (equation 1) in black and the PF model (equation 3) in red.We consider an initially fully excited state and g = 10κ 0 / √ 2. (a) Resonant cavity regime ω a = ω c .For the Dicke and PF model, we consider two values of ω a as indicated in the legend.(b) Non-resonant cavity regime with ω a = 10κ 0 and ω c = 15κ 0 .