Extracting kinetic information from short-time trajectories: relaxation and disorder of lossy cavity polaritons

1 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; and (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.

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.

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 many-body 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 Section 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 Section 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 Section 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 Section 5, we compare the three standard cavity models: the Pauli–Fierz (PF) Hamiltonian, the Dicke Model, and the Tavis–Cummings (TC) model [28]. We show that, in the weak coupling limit, all models converge to similar results, whereas their behavior diverges in the strong coupling limit. We show that lifting the resonance between cavity and atoms can restore the similarity between all three models.

2 Cavity models and transfer tensor method

2.3 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 the density matrix ρ _f (t) of the full system including the cavity

(9) ρ ̇ f ( t ) = − i L f ρ f ( t ) .

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 P

P ρ f ( t ) = Tr C ρ f ( t ) ⊗ ρ C = ρ ( t ) ⊗ ρ C ,

where Tr _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. (9) with the projected identity superoperator of the full system I f

P I f = n 2 E 0 ⊗ ρ C ,

and propagating it to times t _k , P I f ( t k ) , the necessary dynamical maps are simply

E k = Tr C P I f ( t k ) / n 2 ,

from which the transfer tensors may be derived as usual (Eq. (5)).

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 Eq. (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.

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 H TC = h g σ + a + σ − a † . 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.

3 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. (6) and the unilateral z transform of the transfer tensors as defined by

   T z = ∑ k = 1 K z − k T k ,

   with K the index of the last transfer tensor considered,

2. or the Laplace transform of the generator of the Nakajima–Zwanzig Eq. (7)

   L ̃ s ≡ − i L 0 + K ̃ s ,

   where f ̃ ( s ) = ∫ 0 ∞ e − s t f ( t ) d t 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. (8). 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.

3.1 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.

3.1.1 Continuous Laplace transform

To begin, via the Laplace transformation of the quantum master equation, the density matrix is solved formally as:

s − L ̃ ( s ) ρ ̃ ( s ) = ρ ( 0 ) ,

Then, the steady state can be directly computed via the final value theorem

ρ ss = lim s → 0 s ρ ̃ ( s ) ,

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).

3.1.2 Discrete transfer tensors

An equivalent representation of the Laplace formalism is provided by the unilateral z transform of the TTM propagation in Eq. (6) as

ρ z = T z ρ z + ρ 0 ,

where we used T 0 = 0 . The final value theorem takes the form ρ ss = lim z → 1 z − 1 ρ z , so

ρ ss = lim z → 1 z − 1 1 − T z − 1 ρ 0 .

This result is consistent with the fact that the steady state density matrix must not be affected by propagation with transfer tensors by definition

ρ ss = ∑ k T k ρ ss ,

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.

3.1.3 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).

Figure 3:

Steady states inferred via the simple matrix calculation from Section 3.1. All diagrams are constructed for the TC Hamiltonian H TC = ℏ 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, and (d) triply excited N = 4. Under the partially excited states, since the TC Hamiltonian preserves symmetry, the system is prevented from fully decaying to the ground state, instead relaxing to a state that remains partially excited.

3.2 Relaxation rate and lifetime

Next, we extract the rate of convergence to the steady state from the transfer tensors. Formally, the relaxation lifetime τ of a given observable ⟨ o ̂ ⟩ is defined by the zeroth moment of the time-dependent observable relative to its steady state value, i.e.,:

(10) τ = ∫ 0 ∞ 〈 o ̂ 〉 ( t ) − 〈 o ̂ 〉 ( ∞ ) d t 〈 o ̂ 〉 ( 0 ) − 〈 o ̂ 〉 ( ∞ ) .

where ⟨ o ̂ ⟩ ( ∞ ) = ⟨ o ̂ ⟩ ss . Intuitively, this expression can be understood by assuming that ⟨ o ̂ ⟩ ( t ) obeys a simple exponential decay expression of the form ⟨ o ̂ ⟩ ( t ) = ⟨ o ̂ ⟩ ( ∞ ) + [ ⟨ o ̂ ⟩ ( 0 ) − ⟨ o ̂ ⟩ ( ∞ ) ] e − t / τ .

3.2.1 Relaxation matrix

In this section, we generalize the concept of relaxation rate to that of the relaxation matrix

τ ̂ = ∫ 0 ∞ Δ ρ t d t ,

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 o ̂ . Specifically, we define

τ ̂ o ̂ = Tr o ̂ τ ̂ = ∫ 0 ∞ o ̂ t − o ̂ ∞ d t ,

which is directly related to the relaxation timescale τ (Eq. (10)) via normalization, i.e.,

τ = τ ̂ o ̂ o ̂ 0 − o ̂ ∞ .

The relaxation matrix τ ̂ contains 2 ^N real eigenvalues, describing multiple relaxation timescales of the system. We may normalize by the initial condition

(11) τ ̂ ′ = Δ ρ ( 0 ) − τ ̂ ,

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 nonzero 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 nonzero 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.

3.2.2 Continuous laplace transform

In analogy to the spectral method presented in Refs. [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

τ ̂ = lim s → 0 s + L ̃ s − ρ 0 − ρ ∞ ,

where again the symbol − implies the pseudo inverse.

3.2.3 Discrete transfer tensors

Numerically, this calculation can be performed by means of the z-transform of the transfer tensors

τ ̂ ≃ lim z → 1 1 − T z − ρ 0 − ρ ∞ δ t ,

and the limit of z → 1 simply involves the sum of all transfer tensors ∑ _k T _k

τ ̂ ≃ 1 − ∑ k T k − ρ 0 − ρ ∞ δ t .

The transfer tensors are computed based on a time step δt, and the relaxation matrix τ ̂ converges as δt → 0.

3.2.4 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.

Figure 4:

Plot of the relaxation of an individual spin, ⟨σ _z ⟩, of the 2-TLS, 3-level cavity Tavis–Cummings model with H = ℏ g σ + a + σ − a † , g = 10 κ 0 / 2 . (a) For the underdamped case (κ = κ₀), the estimated relaxation timescale τ accurately fits the overall dynamics with the function e^t/τ − 1/2 (green). (b) We consider three increasing values of κ, corresponding to the underdamped (κ = κ₀), critically damped (κ = 15κ₀), and overdamped κ = 1 0 3 κ 0 regimes, respectively.

We now proceed to analyze the relaxation time-scale associated with the dynamics presented in Figure 3. The operator associated with the relaxation measurement o ̂ 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.

Figure 5:

Relaxation rate τ⁻¹ 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 κ.

3.3 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 analyzing their behavior 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 prior knowledge of oscillatory behavior in Section 3.3.1. We continue in Section 3.3.2 by 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.3, we show that an enhanced analysis of oscillatory behavior can be achieved by this method without prior knowledge of the dynamics.

3.3.1 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 o ̂ = σ 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).

Figure 6:

Dynamics of the initially excited TLS σ z t for N = 2 (blue). Time is expressed in units of κ⁻¹. (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 .

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.

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 .

3.3.2 Resonance analysis

The resonance effect as a function of δt is tightly connected to the Fourier transform of Δ o ̂ t = o ̂ t − o ̂ ∞ . Let us explore this connection with an example of an arbitrary observable ⟨ o ̂ ⟩ featuring an oscillatory decay described by

o ̂ t = cos ( ω t ) e − r t o ̂ 0 − o ̂ ∞ + o ̂ ∞ ,

Using transfer tensors with timestep δt to estimate τ, we obtain

τ ( δ t ) = ∑ k cos ( ω k δ t ) e − r k δ t δ t = R e δ t e δ t ( r + i ω ) − 1 = δ t 2 cos ( ω δ t ) − e − r δ t cosh ( r δ t ) − cos ( ω δ t ) ,

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.

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

τ ( δ t ) = ∫ Ξ d t Δ o ̂ t d t Δ o ̂ 0 .

By the properties of the Dirac comb under Fourier transformation, we have

∫ Ξ d t Δ o ̂ t d t = 1 δ t ∑ k F Δ o ̂ 2 π k δ t ,

where F f ω is the Fourier transform of f t . Therefore, the estimate τ(δt) provides a sampling of the Fourier transform of the deviation Δ o ̂ 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 δ o ̂ ω at the value ω = 2 g 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.

Figure 8:

Analysis of the dependency of the relaxation timescale as a function of the timestep. (a) Estimate τ for σ _z as a function of δt for N = 2 and different values of κ. As κ increases, the relaxation rate τ⁻¹ at the peak increases too. Parameters: g = 10 κ 0 / 2 . (b) Rate estimate τ⁻¹ as a function of κ for N = 2, showing a good agreement with the curve κ/4. Parameters: g = 10 κ 0 / 2 .

3.3.3 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 6 g . 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 = 6 g and ω₂ = 2ω₁. By inspecting the time dependence of the system, Figure 9a, it becomes clear that ω₁ matches exactly the oscillation frequency of the dynamics and hence maximizes the estimate τ. The double frequency ω₂ 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.

Figure 9:

Analysis of resonances of the relaxation timescale as a function of the timestep for the probability of the fully excited state. (a) Estimate τ of the probability of the fully excited state as a function of δt for N = 2 (blue) and a fully excited initial state. Two types of resonances can be recognized: one at odd multiples of δ t κ = π / 10 3 (vertical red dotted line) and another one at multiples of δ t κ = π / 5 3 (vertical black dash-dotted line). (b) Probability of the initial state as a function of time. The period of the oscillations matches exactly δ t κ = π / 5 3 , as shown by the black plus sign. Half the period is indicated by red crosses. Parameters: g = 10 κ / 2 .

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.

Figure 10:

Analysis of resonances of the relaxation timescale as a function of the timestep for the expectation value of Pauli z operator. (a) Estimate τ for operator σ _z as a function of δt for N = 2 (blue) and a fully excited initial state. Beyond the resonances already present in Figure 9a (vertical gray dotted lines), two additional resonances can be observed at δ t κ = 0.17 π / 3 (vertical black dash-dotted line) and at δ t κ = 0.18 π / 3 (vertical red dash-dotted line). (b) Expected value of operator σ _z as a function of time. The two extra resonances roughly correspond to the beating of the dynamics: the first, destructive resonance, corresponds to the half-period of the beating (black dots), whereas the second resonance corresponds to the period of the beating (red crosses). Parameters: g = 10 κ / 2 .

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.

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, T ̄ k is the kth transfer tensor used for propagation of the average dynamics, and T k i is the kth transfer tensor computed from the ith realization. (ii) We average the dynamical maps computed via each individual realization, i.e.,

where E ̄ k is the average kth dynamical map used for computing the transfer tensor, and E k i 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 recovers 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 E ̄ 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 T ̄ k .

Figure 11:

Comparison of simulated spins for various initial states by averaging over 50 realizations. We want to reproduce the behavior resulting from averaging the density matrices of each realization, and two attempts to capture that behavior are given by averaging the dynamical maps or the transfer tensors. All simulations have disorder given by choosing ω _j /κ ∈ U(40, 50), for the 2-TLS 3-level cavity TC model with H = ℏ g ( σ + a + σ − a † ) + ℏ ω a ∑ j = 1 N σ j x + ℏ ω c a † a and g = 10 κ / 2 , ω _c = 50κ. (a) Fully excited state. (b) Singly excited state. In this scenario, the true dynamics of the system are markedly different from those deduced via the disorder TTM, and this can be interpreted as the TTM being unable to recover any coherent effects from asymmetric initializations. (c) Coherent superposition initial state, i.e., all spins are in the coherent superposition ( | ↑ 〉 + | ↓ 〉 ) / 2 . Once again, averaging the density matrices and transfer tensors give good agreement with both each other and the symmetric Hamiltonian dynamics, just as in the fully excited initial state. Tensors are computed with learning time κt = 4 and with 10× lower resolution than the density matrices, leading to the figures not matching exactly in (c).

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 allow us to deduce the averaged behavior of the system using fewer realizations.

Finally, we compute the largest eigenvalue of the normalized relaxation matrix τ ̂ ′ Eq. (11) 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.

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 = h g ( σ + a + σ − a † ) + h ω a ∑ j = 1 N σ j x + h ω c a † a and disorder given by ω _a ∈ U(45 − δ, 45 + δ), ω _c = 50, and g = 10 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.

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.

5 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, 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 τ⁻¹ 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.

Figure 13:

Relaxation rate of σ _z as a function of cavity decay rate κ for three models: TC (Eq. (2)) in blue, the Dicke model (Eq. (1)) in black, and the PF model (Eq. (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) Nonresonant cavity regime with ω _a = 10κ₀ and ω _c = 15κ₀.

Let us first discuss the case of ω = 100κ₀. 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κ₀ 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 toward 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κ₀, but increase the cavity frequency to ω _c = 15κ₀. 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.

6 Conclusions

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 toward 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 nonmonotonic 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 partially 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, and nonadiabatic quantum dynamics [3], [8], [9], [10], [11], [12], [14]. Additionally, while we have shown the immediate application toward 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.