Non-Markovian Dynamics of Open Quantum Systems via Auxiliary Particles with Exact Operator Constraint

We introduce an auxiliary-particle field theory to treat the non-Markovian dynamics of driven-dissipative quantum systems of the Jaynes-Cummings type. It assigns an individual quantum field to each reservoir state and provides an analytic, faithful representation of the coupled system-bath dynamics. We apply the method to a driven-dissipative photon Bose-Einstein condensate (BEC) coupled to a reservoir of dye molecules with electronic and vibronic excitations. The complete phase diagram of this system exhibits a hidden, non-Hermitian phase transition separating temporally oscillating from biexponentially decaying photon density correlations within the BEC. On one hand, this provides a qualitative distinction of the thermal photon BEC from a laser. On the other hand, it shows that one may continuously tune from the BEC to the lasing phase by circumventing a critical point. This auxiliary-particle method is generally applicable to the dynamics of open, non-Markovian quantum systems.


I. INTRODUCTION
Experimental platforms coupling a set of photonic cavity modes to an ensemble of dye molecules, comprised of two-level systems (TLS) of electronic excitations and local vibrational degrees of freedom, are relevant for applications ranging from Bose-Einstein condensates (BEC) of photons [1][2][3][4][5], exciton polaritons [6,7] and plasmonic lattices [8] to single-photon sources for quantum information [9].Despite general non-equilibrium field theory being available [10], the stationary states and dynamics of such open, driven-dissipative quantum gases are unexplored to a large extent.Recently, measurements of the secondorder coherence in photon BECs [11,12] have opened the door to a deeper understanding.In fact, the experiments revealed a hidden, complex structure in the temporal correlations which undergo a non-Hermitian phase transition as the system is driven away from equilibrium, while the spectral mode occupation still follows an equilibrium Bose-Einstein distribution.
As open quantum systems are driven far from equilibrium, reservoir frequency scales like the spectral substructure [13] and relaxation rates can no longer be considered fast compared to the photonic cavity loss and tunneling rates.In particular, in multiple coupled cavities [14], fast Josephson oscillations can approach the time scales of such reservoir processes.In these cases, quantum coherence and non-Markovian memory effects [9,15] are expected to become important.The theoretical treatment of this non-Markovian regime is complicated by the vast number of reservoir degrees of freedom and the necessity of describing their full quantum dynamics.Direct numerical time-evolution methods are limited to rather small systems, and the extension to larger systems requires truncation of the quantum entanglement between the subsystems [16].On the other hand, approximate techniques based on the Born-Markov approximation and rate equations [17][18][19][20][21], well established for Markovian dynamics even in large systems, inherently cannot capture non-Markovian memory or quantum coherence effects.Standard field-theoretic techniques suitable for a large number of degrees of freedom, like the quantum Langevin approach [22], are hampered by the non-canonical statistics of the pseudospin operators involved in the electronic TLS dynamics, which have, for this reason, been approximated by bosons [23].A recent mean-field model of organic polartions [24] rests on the assumption of symmetry breaking for the cavity mode and is thus limited to large photon numbers.Few-photon effects such as decondensation cannot be described by this approximation.Cumulant expansions such as employed in [25] have been shown recently to be difficult to control [26].
In this work, we develop a general field-theoretic technique based on an auxiliary-particle representation, which captures the non-canonical reservoir dynamics exactly, to describe a large class of coherent, light-matter coupled, open quantum systems such as the Holstein-Tavis-Cummings Hamiltonian or arbitrary TLS in a cavity, and generalize it to strong, time-dependent nonequilibrium.In this approach, an individual canonical quantum field is assigned to each matter state (electronic and vibronic molecule excitations), with the constraint that the system occupies exactly one of these states at any instance of time.Because of the canonical commutation relations, the auxiliary-field dynamics are amenable to standard field-theoretic techniques, and the exact confinement to the physical Hilbert space defined by the above-mentioned constraint is performed analytically as long as correlations between different molecules can be neglected.The auxiliary-particle method has been pioneered by Abrikosov [27] and Barnes [28,29] and is successfully used [30][31][32][33] for strongly correlated electron systems like the Anderson impurity model [34].The extension to stationary non-equilibrium has been developed in [35][36][37][38], and a more recent exposition for correlated electrons in time-dependent non-equilibrium can be found in [39].Generalizing this to quantum optics, we find an exceptional point in the dynamics of the second-order coherence for large losses, and uncover a non-equilibrium decondensation phase transition indicated by critical slowing down.We demonstrate how non-Markovian effects become relevant in the strongly driven-dissipative regime.

We consider a general class of Hamiltonians of the form
where describes a set of cavity modes with dispersion ω k alongside M molecules comprised of an electronic TLS with splitting ω D and a phonon mode Ω (the generalization to multiple phonon modes is straightforward).The Jaynes-Cummings (JC) coupling in the rotating-wave approximation is as usual given by In these expressions, a † k , b † m are the creation operators for a photon or a phonon quantum, respectively, and σ z m , σ ± m are Pauli operators describing the electron dynamics.Electron-phonon coupling results from the phononicoscillator displacement xm = b m + b † m depending on the electronic TLS state and reads

A. Auxiliary-Particle Representation
A single molecule has quantum states |σ, n⟩, where σ = g, e refers to the electronic ground and excited state, respectively, and n denotes the vibrational state.For each of these states, we introduce auxiliary or pseudoboson operators The molecular operators may then be expressed as where we have dropped the subscripts m.This representation is faithful within the Hilbert space spanned by the states |σ, n⟩ (no product states of these), i.e., under the operator constraint that the total auxiliary particle number obeys In the frame rotating with the electronic frequency ω D , the Hamiltonian for M = 1 thus reads where δ k = ω k − ω D is the resonator detuning.We note in passing that in pseudo-boson representation, the molecular part of this Hamiltonian is quadratic.That is, the polaron transformation diagonalizing this part is straightforward and commutes with the constraint Q = 1, but transforms the molecule-photon coupling into mn γ mn (a † k d † g, m d e, n + h.c.), where the coupling matrix γ mn may be calculated analytically and is proportional to the Franck-Condon integrals.In the following, we will use the non-diagonal representation Eq. ( 7) for quantitative evaluations.

B. Hilbert-Space Projection
Since the Hamiltonian (7) conserves the operator Q, any physical expectation value of operators acting on the molecular Hilbert space (e.g.combinations of b, b † or σ ± that annihilate the empty sector Q = 0) can be projected onto the sector Q = 1 by first inserting a factor of ζ Q at the beginning of the Schwinger-Keldysh contour [40], i.e. by replacing the initial density matrix ρ0 by ζ Q ρ0 , where the fugacity is ζ = exp (−µ), and the chemical potential µ can be gauged time-independent.To this end, we write the canonical partition function on a Hilbert space H Q with any fixed value of the operator constraint as Tr H Q ρ0 = Z C (Q).Then the grand-canonical partition function becomes where the decisive factor ζ Q weighs each canonical term according to the operator constraint.The expectation value associated with Z G is When X annihilates the empty sector Q = 0, its unphysical contribution is removed, and we find the correct average in the physical subspace H 1 as It is now instructive to look at two special operator averages at t 0 = 0, i.e.
Defining the greater and lesser Green functions for the molecular ground states, we conclude that these scale as at all times as long as the interactions are treated within a conserving approximation as done below.The same holds for the excited-state Green functions Note that, while the greater functions by themselves do not annihilate the ground state, they still occur in Wicktype expansions of physical expectation values such as ⟨σ + (t)σ − (t)⟩.For the special case of a TLS without vibrational states (i.e. the Jaynes-Cummings model discussed in Appendix E), we may calculate physical expectation values such as ⟨σ z (t)⟩ by applying the Hilbertspace projection according to where we have dropped the vibrational subscripts (m, n = 0) and explicitly pulled out the fugacity.Now, since any product of a greater and a lesser auxiliaryparticle Green function is of order ζ, Eq. ( 14) recovers the correct result.Note that arbitrary products of timeordered Pauli matrices can be treated in the same way.
For the equal-time average, we hence find the expected result which is plotted in Fig. 11 of Appendix E. Note also that we have used the fact that the equal-time greater Green functions follow Eq. ( 25).The projection of the molecule dynamics is thus imposed by the simple rule that, in any auxiliary-particle Feynman diagram, only terms of leading order in ζ are retained.The projected Dyson equations then reduce to where h G is the non-interacting part of the Hamiltonian (7) in the molecular ground-state sector, and products of bold-faced quantities are understood as matrix products.Analogous equations hold for the excited-state propagators E ≶ (t, t ′ ).Observe how in the first of Eqs. ( 16), terms containing two auxiliary-particle lesser functions have vanished, while in the second only terms without lesser functions remain.The photon Green functions, defined as do not participate in the ζ-scaling, since they do not operate in the auxiliary-particle space.Thus, their Dyson equations have the usual form, where we have defined We also introduce the occupation of the photon ground mode as with steady state N = N (t → ∞).

C. Conserving Approximation
As mentioned, a conserving approximation is necessary to implement the Q conservation.It can be generated from a Luttinger-Ward functional [33,41,42].As we show in Appendix E, other approximation schemes such as heuristic truncations of the cumulant hierarchy [20,25] can lead to unphysical effects excluded here by construction.To second order in the coupling g, one obtains the analog of the non-crossing approximation [33], where the self-energies become which is depicted diagrammatically in Fig. 1 (c).The factor ζ −1 in Σ ≷ D stems from the normalization of the physical average.In any self-energy insertion to D ≷ appearing in an auxiliary-particle diagram, such a normalization is not present.Consequently, such insertions from the same molecule vanish by the projection ζ → 0. For arbitrary numbers of molecules M , the projection onto the physical subspace Q = 1 is done as described above separately for each molecule m since Q is conserved locally.Coherence between different molecules may nevertheless be mediated by the photon modes, which couple to all the molecules.For multiple molecules, the right-hand side of Eq. ( 18) acquires an additional prefactor M , and the photon propagators within the self-energies of any particular molecule m become renormalized by self-energy contributions from all the other molecules m ′ ̸ = m.

D. Driven-Dissipative Photon-Molecule Dynamics
While the molecule-bath dynamics are treated fully coherently, we incorporate the external drive and loss by adding Lindblad terms to the master equation for the density matrix, where L[X]ρ = XρX † − {X † X, ρ}/2.The cavity loss is κ and L M = L ↑,↓ + L relax , where the external drive and radiationless decay of electron excitations are described by (c.f.Fig. 1 (b)) while the process of phonon relaxation due to molecular solvent collisions is captured via where λ is the relaxation rate of the molecular vibrations, and the solvent temperature enters implicitly through the average phonon occupation number n(Ω).To our knowledge, the exact projection method employed here has not, until now, been applied to open quantum systems.Therefore, we provide more details on how this generalization is to be performed.Note also that it was not a priori clear that such a generalization should be possible.
The Lindblad operator (23) describing the relaxation of the vibrational states can be understood by investigating the structurally equivalent operator where Γ is now an arbitrary constant.The projection is again performed by removing all terms with too many lesser functions, where one should keep in mind that for auxiliary particles there holds for all times T .We remark again that these properties are not an approximation, but exact identities under the projection method.We arrive at These equations evidently possess the correct fugacity scalings.For the equal-time equation of G < , once more we remove all terms with the wrong scaling.This yields where we have introduced an anti-commutator.Again, this equation has the correct fugacity scaling.As mentioned above, we do not need an equation for the equal-time evolution of G > .Within both the groundand excited-state manifolds, the Lindblad operators of Eq. ( 23) couple neighboring vibrational states according to a structure which for the self-energy matrices of Eqs.(26) and n + 1 states in total looks like We still have to consider the external pumping and electronic loss terms The relevant matrices for this read In terms of these, the contribution to the equations of motion can be written rather compactly as For the forward dynamics, we find i Ė< (T, 0) The equations for the electronic loss Γ ↓ follow analogously.The functional form of the total external pumping reads where t P varies with the attainment of the steady state, and A = 10 −2 , λσ P = 1/2.The initial conditions are FIG.2. Relaxation of the vibrational degrees of freedom in the electronic excited state.The dashed line shows − Im pmax n=0 nE < nn (t, t) → n.The parameters are g = 0, i.e. we have decoupled the excited states from the rest of the system, ωD/λ = 2.0, Ω/λ = 0.1, S = 0.5, n = 1.0,Γ ↑ = Γ ↓ = 0.The initial conditions is iE < 00 (0, 0) = 1 with all other lesser functions being zero.The number of vibrational states considered is 10, which is equivalent to a maximal number of phonons pmax = 9.The results were calculated with a step size of 10λ • 2 −9 for a number of 2 9 steps.N (t) = 0, with D > 00 (0, 0) following from the commutator, G < mn (0, 0) = −iδ 00 and E < mn (0, 0) = 0.The greater functions G > mn (0, 0) = E > mn (0, 0) = −iδ mn are fixed by the projection.
Before concluding this section, consider Fig. 2 which shows the relaxation of the vibrational degrees of freedom in the electronic excited state induced by the vibrational relaxation in Eq. ( 23).The dashed line indicates the statistical average − Im n nE < nn (t, t) which correctly approaches n.We therefore conclude that our master equation is capable of describing the phonon thermalization adequately.

III. RESULTS AND DISCUSSION
We solve the Kadanoff-Baym Eqs.(16,18) selfconsistently with Eq. ( 20) using a multi-step predictorcorrector method [43,44] developed earlier to facilitate adaptive two-time evolution [45].We truncate the time integrals at the memory time τ mem , which is set by the inverse, dissipative relaxation rates κ, λ (Appendix D).This allows for an efficient yet accurate simulation.For small system sizes, we find quantitative agreement with quasi-exact numerical time evolution [16].As a further consistency check, in the Markovian limit (phonon relaxation rate λ faster than any other scale in the system) our approach reduces to the previously studied semi-classical rate equations [17], with effective absorption and emission coefficients as derived in Appendix A. Γ(±δ k ) encodes the molecular absorption and emission spectra shown in Fig. 4 as a function of δ ≡ δ k .Around the zero-phonon line, δ = 0, these spectra satisfy the Kennard-Stepanov relation at the temperature fixed implicitly via n(Ω) (see parameter values in Fig. 3).For the numerical evaluations, we consider a low solvent temperature such that n(Ω) = 0.25 which justifies truncating the phonon occupation numbers at n ≤ n max = 4.
Previous experiments [11] revealed that drivendissipative photon condensates possess hidden, nontrivial dynamics in the second-order coherence despite the spectrally resolved photon number following an equilibrium Bose-Einstein distribution.The g (2) (t) oscillation frequency and decay rates as functions of Γ ↑ and κ, i.e. of the distance from true equilibrium, show a non-Hermitian phase transition marked by an exceptional point [11].Deviations from this picture, obtained on the basis of the rate equations [12,17], are to be expected when the system is strongly out of equilibrium.Our formalism is ideally suited for studying such non-Markovian effects.
We start the time evolution with an empty singlemode cavity filled with an ensemble of molecules in the ground state and drive the system into a steady state by a constant optical pumping Γ ↑ (inset of Fig. 3(a)).Once stationarity is reached, a short Gaussian pulse is added (Eq. ( 32)) to trigger a response of the photon field as the system is slightly displaced from the steady state.Via quantum regression [46], the ensuing relaxation is then equivalent to the spontaneous intensity fluctuations described by g (2) (t).As shown in Fig. 3(a), the photon relaxation changes qualitatively from oscillatory to bi-exponential as a function of κ (observe that κ/Γ ↑ is kept fixed) [11,12].This behavior is generic and independent of the specific parameters chosen.For each set of parameters, the transition between the two behaviors is characterized by an exceptional point where the eigenvalues of the linearized regression dynamics [12] coalesce and switch from a pair of complex frequencies ±iω − τ −1 to two real relaxation rates τ −1 , τ ′ −1 > τ −1 [47].We extract these parameters by fitting the sum of two complex exponentials to the numerical results Appendix C. We find that the phase boundary between both behaviors retraces qualitatively the one obtained from the rate equation model in Appendix A. Therefore, in Fig. 3 pare in Fig. 5 the full quantum dynamics, Eqs. ( 16)-( 20), with the rate-equation results along the red, dashed line shown in Fig. 3 (b).Evidently, the rate-equation model agrees with the full quantum dynamics only near equilibrium (κ/λ ≪ 1, Γ ↑ ≪ 1), while for strong drive and dissipation deviations occur, indicating a breakdown of the Markov assumption inherent to the rateequation approach: the system becomes non-Markovian in this stronger sense, as opposed to merely a frequencydependent photo-molecular coupling [13].As seen from Fig. 5 (a), the steady-state occupation N collapses beyond a critical loss rate κ in spite of constant κ/Γ ↑ [4].This seemingly counterintuitive effect is due to the loss of coherence with increasing κ and Γ ↑ .This decondensation is accompanied by critical slowing down, i.e., a diverging relaxation time, τ −1 → 0 (Fig. 5 (c)), a clear sign of a non-equilibrium phase transition , while the fast relaxation rate τ ′ remains finite.Since τ −1 → 0 is only possible when both eigenvalues are real (ω = 0) [12], the decondensation threshold must always occur in the biexponential phase, c.f. Fig. 5 (a), (b).
We furthermore notice that the strongly driven, nonequilibrium regime (upper right in Fig. 3 (b)) represents a laser.That is, the near-equilibrium photon BEC phase is qualitatively separated from the laser phase by a hidden g (2) (t) transition, which can be circumvented via the path in parameter space shown in the inset of Fig. 3 (b).
Finally, in Fig. 6, we study the non-equilibrium dynamics of the effective molecular emission and absorption spectra, as seen by the photons, for stronger coupling and at lower temperature.These results are obtained without further approximation such as memory truncation.Our method is thus capable of temporally resolving the emergence of characteristic spectral features (s. also Appendix F).  16)-( 20); red, dashed curves: rate-equation model (Appendix A). λΓ − 0 = 0.796g 2 , Γ + 0 /Γ − 0 = 0.741 (obtained from Fig. 4).The vertical line marks the exceptional point.

IV. CONCLUSION.
We have introduced a non-equilibrium auxiliary field theory which faithfully describes the time-dependent quantum dynamics of general coupled system-bath setups, here generalized to open, driven-dissipative quantum systems such as photon BECs coupled to dye-molecule reservoirs [1] or exciton-polariton systems [48].Our method may also be applied to the full quantum dynamics of multiple qubits coupled to sources of non-Markovian noise [49].
For the open photon-BEC system, we find significant non-Markovian memory effects in the strongly driven regime where system loss κ and reservoir relaxation λ become comparable.We uncovered the global shape of the phase diagram partially explored in [11].These calculations establish that the near-equilibrium photon BEC is separated from the lasing regime by a hidden phase transition of the photon density response.vibrational frequencies and relaxation.It has the property g nn ′ (0) = G < nn ′ (t, t)/G < (t, t), such that g nn (0) = p n .The analogous expression for the excited states is e nn ′ (0) = E < nn ′ (t, t)/E < (t, t).We furthermore introduce the total number of excited molecules as and similarly for the ground state where such that the constraint Q = 1 ensures total number conservation.Finally, the above assumptions ensure that the greater functions do not possess a forward-time dependence, i.e.
With these definitions, and using that under the approximation described above one can pull out the photon propagators from the time integrals and normalise by the respective occupation number of the ground and the excited states, we may transform Eq. (A4) into which is the established result [17,18], and we have identified Γ ± k = Γ(±δ k ) = 2 Re K(±δ k ).In our case, the emission and absorption coefficients are then given via where (A11) This is to be compared with the original expression in the quantum master equation derived via the usual Born-Markov approximation [17,18]: where f (t) is the polaron correlation function.An alternative and more practical way of defining Eq. (A10), which is also leading directly to Eq. (A9), follows by letting Defining K(δ) in this way occurs with the thought in mind that it is effectively the right-hand sides of Eq. (A13) acting instead of their previous definitions [17,18] when parameters are chosen such that the auxiliaryparticle dynamics effectively recovers the rate equations.Finally, the molecular rate equation follows completely analogously and reads Appendix B: Driven-Dissipative Processes in Schwinger-Keldysh Formalism To be self-contained, here we give a brief introduction to the Schwinger-Keldysh formalism as it applies to open systems.Subsequently, we show how to derive those parts of the equations of motions for the greater and lesser Green functions deriving from the introduction of driven-dissipative processes.We emphasize that we make no approximations here other than the formal diagrammatic expansion in the auxiliary particles.The projection as such is exact.Furthermore, the resulting equations of motion turn out to be effectively quadratic, which hints at the fact that their validity is rather robust also for strong drive and dissipation.
As the simplest example, consider a cavity mode of frequency ω 0 , coupled to an environment at a low temperature such that ℏω 0 ≫ β −1 , is described by the Lindblad master equation with the non-Hermitian Hamiltonian Ĥ = (ω 0 − iκ/2) a † a. (B2) In its most general formulation as given by Schwinger [50], the formalism easily captures this kind of dynamics.The Schwinger action that appears in the coherent-state expansion of the non-equilibrium partition function Z is  Note how this contains a contribution across the two branches of the contour.The respective equations of motion for the greater and lesser Green functions are then 0 where now the time-ordered Green function appears explicitly.The anti-time-ordered Green function would appear, for instance, when coupling to a bath at finite temperature.In the equal-time limit, these equations become which serves to illustrate how the commutator is preserved over time by the quantum jumps.
To have a consistent diagrammatic expansion for the four vertices g, Γ ↑, ↓ and λ of our theory, we expand the corresponding vertices to second order in ℏ.For the incoherent couplings, such a two-loop expansion amounts to working in Hartree-Fock approximation.The density responses in the main panel of Fig. 3 (a) show the evolution of iD < (t, t) following the Gaussian pulse (32), where the point t = 0 of the plot corresponds to the maximum of the curve after the Gaussian pulse has been injected.The frequency ω and decay rate τ in the oscillatory phase and the decay rates τ 1 , τ 2 are extracted by least-squares fitting the functions to the numerical density responses.The quality of the fit is then estimated by the standard error, and the oscillatory or bi-exponential fit is accepted, whichever has the smaller error.This determines whether the response is classified as oscillating or bi-exponentially relaxing.An illustration of the fitting procedure is given in Fig. 7 and Tabs.I, II.One can see that even for these two points directly next to the transition, the fit classification is always unique.

Appendix D: Memory Truncation
Solving the equations of motion for interacting systems numerically on the two-time grid becomes computationally expansive when the number of grid points needs to be large.This happens whenever the product of the fastest system frequency and the required final time is not small.Any integral over relatively fast decaying Green functions, however, will be computable at sufficient precision over a small support that does not grow as the two-time grid expands towards it final size.In the same spirit, computing points far off the "time diagonal" t = t ′ will be superfluous because they are negligible.In this way, by truncating the number of steps one moves away from the time diagonal, and by restricting the computation of the integrals to that same narrow "band", it is indeed possible to achieve a quasi-linear scaling in the number of grid points without losing accuracy.
The quantitative consequences of different truncation parameters n mem are studied in Fig. 8.Note that technically, n mem is not the number of points included away from the diagonal moving in the direction (t − t ′ ), but rather moving in the vertical and horizontal directions t and t ′ .A short-valued memory n mem = 24 results in a certain number of points inaccurately remaining zero.With a longer-valued memory n mem = 64, these points attain their proper values.However, looking to Fig. 8, we see that the influence of this is not considerable.
These are, however, not closed because of the occurrence of ⟨a † aσ z ⟩ in the last of Eqs.(E2).In the liter-ature, in generalizations of H JC to many spins, this term has been treated by means of the heuristic truncation ⟨a † aσ z ⟩ ≈ ⟨a † a⟩⟨σ z ⟩ [20], which is not a controlled approximation.As shown in Fig. 9, this can even result negative particle numbers (not a numerical artifact).While this problem is mild for the parameters chosen here, it becomes worse in other regimes.In particular, the factor of M coming from the sum over many molecules in the equation for ⟨a † a⟩ can lead to unphysical results.This is a manifestation of the difficulties involved in truncating the expectation-value hierarchy when phase coherence between photons and molecules is important.As explained above, the auxiliary-boson technique provides a resolution of this problem.Here, this can be appreciated from Fig. 10.The inset highlights that the unphysical negative particle numbers have been removed by the auxiliary-particle method and that it reproduces numerically exact results.
Appendix F: Time-Dependent Emission And Absorption Spectra The photon self-energies as defined in the main text, i.e.Σ ≷ D (t, t ′ ) = ig 2 ζ −1 Tr E ≷ (t, t ′ )G ≶ (t ′ , t), (F1) can be used to define the effective emission and absorption spectra via which the molecules act on the photons.This can also be understood by comparing the above expression for the self-energies with Eq. (A13), which we defined as (F2) Here we show more data complementing Fig. 4 in the main text.The evolution toward the steady state of the occupation numbers if shown in Fig. 12, while additional cross-section through the effective spectra presented in the main text's Fig. 4 are given in Fig. 13.To obtain these curves, the self-energies are first evolved with their full dependence on (t, t ′ ) and then transformed to socalled Wigner coordinates T = (t + t ′ )/2, τ = t − t ′ .The latter relative-time variable is then Fourier transformed to obtain the spectra shown in Fig. 13.
Appendix C: Determining the Oscillation Frequency and Decay Rates of the Second-Order Correlations

6 FIG. 7 .FIG. 8 .
FIG. 7.Comparison of the fits directly to the left and right of the transition at large κ, Upper panels: Data and best fits for several different ansatz functions.A uni-exponential decay is also fitted to underline that it is not a possible best fit, as can indeed be judged by eye.Lower panels: The absolute error (difference of data and fit) for the two viable options fosc and f bi-exp .

TABLE I .
Standard errors for the damped-oscillating ansatz.

TABLE II .
Standard errors for the bi-exponentially decaying ansatz.