Relationship between Population Dynamics and the Self-Energy in Driven Non-Equilibrium Systems

We compare the decay rates of excited populations directly calculated within a Keldysh formalism to the equation of motion of the population itself for a Hubbard-Holstein model in two dimensions. While it is true that these two approaches must give the same answer, it is common to make a number of simplifying assumptions within the differential equation for the populations that allows one to interpret the decay in terms of hot electrons interacting with a phonon bath. Here we show how care must be taken to ensure an accurate treatment of the equation of motion for the populations due to the fact that there are identities that require cancellations of terms that naively look like they contribute to the decay rates. In particular, the average time dependence of the Green's functions and self-energies plays a pivotal role in determining these decay rates.


Introduction
Non-equilibrium many-body physics is a complex problem because it requires the determination of two-time Green's functions within a system that does not have time-translation invariance. One pathway to simplify this approach has been to investigate the populations of electronic states (as a function of momentum and time) since they depend only on one time variable. Evaluating the differential equation of motion, which determines how the populations evolve with time, reveals that they are complicated by "memory effects", given by the state of the system in the recent past, and hence they involve integrations over past times. This non-locality in time arises entirely from the fact that we need to represent the two-particle averages (that typically give the potential energy of the system) via convolutions of the self-energy with the Green's function. One route to searching for a description that is local in time is to work on directly determining the two-particle Green's functions without resorting to this convolution. We do not pursue that approach here. Instead, we focus on examining the structure of the population dynamics and strive to understand as much as we can about the exact nature of these equations before we delve into approximate treatments. We also compare the approximate results of the equations of motion for the populations with direct numerical results for the Green's functions and the self-energies that are found by self-consistently solving the problem for the Green's functions within the Keldysh approach. This shows how the different scattering integrals behave as functions of time and how they ultimately determine the relaxation of the populations.
In recent years, non-equilibrium dynamical mean-field theory and closely related approaches [1,2], have solved numerous many-body physics problems that are driven by external fields to model ultrafast pump-probe experiments [1,[3][4][5][6][7][8][9][10][11]. In those solutions, it is often seen that the populations relax exponentially to their equilibrium values, and in some cases, one can identify the relaxation rate as given by the imaginary part of the equilibrium self-energy. Furthermore, a simple linear-response analysis of the relaxation of a Green's function shows that the relaxation rate is dominated by the imaginary part of the analytic continuation of the self-energy into the lower complex plane at the location of the pole in the Green's function that lies closest to the real axis; if the system is described by a Fermi liquid, then the relaxation rate for electrons lying near the Fermi energy is given essentially by the imaginary part of the self-energy at the excitation energy above the Fermi energy [12].
The fundamental question we address here is to determine the relationship between the relaxation rate of the populations and the self-energy for more general non-equilibrium cases. One point, which makes the analysis complex, is that it is the average time dependence that determines the relaxation of the system (because systems with no average time dependence have time translation invariance, and do not evolve with time). But it is the relative time dependence which is most closely related to the frequency dependence of the self-energy because the two are related via Fourier transformation. Reconciling these two issues is the primary goal of this work. While we shed light onto this problem, we do not completely determine an analytic form for the population decay for the general nonequilibrium case.  The analysis of the two-time Green's function is done on the Keldysh contour, illustrated in Fig. 1. The imaginary spur allows for the calculation of an equilibrium thermal state at t min , which is subsequently propagated forward in time along the contour using the equations of motion discussed below. For convenience, we will assume t min = 0. The system is then subjected to an external perturbation (a field) at some time t > 0. The Green's function for a state with quasi-momentum k is defined on the contour as

Equations of motion
where T C is the contour time ordering operator, andĉ † k /ĉ k are the usual creation and annihilation operators. The superscript C denotes a contour-ordered quantity. As discussed in Ref. [13], the Green's function satisfies the equation of motion on the Keldysh contour: as well as a similar (adjoint) equation for t . Here, the integration is carried out over the entire contour. Σ encodes the interactions in a diagrammatic sense, and is a functional of the Green's function. For simplicity, we have chosen a momentum-independent self-energy.
Based on the location of the two times, we can break up the contour-ordered Green's function into different temporal components. The real-time Green's functions are the lesser, greater, retarded, and advanced components(<, >, R, and A). They are related through: The mixed components, when one of the arguments (τ) is on the imaginary spur, and the other (t) on the real time branch, are denoted as G (t, τ) and G (τ, t). Finally, the component which has both arguments on the imaginary spur is the usual Matsubara Green's function G M (τ, τ ). These subdivisions are applied for all contour-ordered quantities. We can carry this through for the equation of motion, which splits into 3 distinct pieces, and similarly for the adjoint equation. For completeness, we list them here, following the notation of Ref. [13]. Since time-translation invariance applies on the imaginary axis, the equation of motion is significantly simplified: The lack of dependence on the real time integrations is a reflection of causality -the equilibrium state on the imaginary spur should not depend on any quantities that occur later in time. The equations involving real times are: where the derivative applies following the direction of the arrow. The terms on the right hand side are the scattering integrals, which we will focus on: The scattering integrals have two distinct sets of terms: those involving values on the imaginary axis, and those that do not. In the presence of any type of interactions, the Green's functions decay in the t − t direction. This leads to a decay in the contribution of the mixed pieces to the integrals as time gets further from t = 0. We will assume that t and t have advanced sufficiently far from t = t = 0 that we can neglect those mixed pieces in our analysis. The density for momentum k is given by the time-diagonal piece of the lesser Green's function, To propagate along the t = t (t ave ) direction, we need to combine Eqs. (8) and (9). This yields which can be simplified by noting that For single-band models on lattices without a basis the Green's function can be described by a set of scalars and the energy k (t) commutes with the Green's function. We can further assume that we have gone far enough along in t and t that we can extend the lower bound of the integral to −∞ due to the Green's function decay, and are left with where in the last line we have introduced the · (as a non-Abelian operator) to represent the integral for notational convenience.

General remarks on the scattering integrals
It is instructive to consider the contour-ordered quantities as functions of the average and relative time, rather than the explicit t and t : t ave = 1 2 (t + t ) , t rel = t − t . The population dynamics occurs along the t ave direction, whereas quasiparticle lifetimes arise from the t rel direction. The scattering integrals control the temporal dynamics of the Green's function along the t ave direction. Thus, an examination of the integrals can lead to some general insights regarding the dynamics. First, it can be shown that an average time dependence is necessary for any change to occur -that is, without any dependence on t ave , the scattering integrals (and thus the time derivative of the density) are 0.
It is illustrative to consider the equations in a form where the Fourier transform over t rel has been performed. For this transform to be well-defined, we must assume that any average time dependence is slow enough within the window of Fourier transformation set by the decay of the Green's function along t rel that we can replace the average time dependence with a single time only, namely t ave . Any applied field is assumed to have occurred sufficiently far in the past that it falls outside this window. This yields When there is no dependence on t ave , we can use Dyson's equation and substitute for G < k (ω) and ImG R k (ω). Long after t = t = 0 where the mixed components have decayed, we can re-write Dyson's equation [Eqs. (8) and (12)] in the following way: Making the same transformation to t ave , t rel and Fourier transforming over t rel , we can substitute this result into Eq. (19): The point here is that the lack of time dependence in the density comes from the cancellation of two terms, rather than each term being 0 individually. This will come up again when we consider specific cases.
The central result is even stronger than this: once the lesser Green's function can be represented by any distribution function of average time and frequency that multiplies the imaginary part of the retarded Green's function at the same average time and frequency, that is G < and similar for the self-energy), then the population no longer changes with time. The factor of two arises from the connection to equilibrium; this assumption is true if the system has equilibriated at an elevated temperature T el , where the fluctuation-dissipation theorem holds and f (ω, T el ) is the Fermi function. We can readily see this from Eq. (19) by substituting this relation: Note that this cancellation occurs for all types of scattering and for all strengths of the scattering. It does not require any so-called "bottlenecks" but simply follows from the dynamics of the populations and the exact form of the scattering integral. A further implication of this cancellation is that a simple hot electron model cannot fully describe the dynamics because approximating G < k (ω, t ave ) and Σ < (ω, t ave ) by the retarded components multiplied by the same average-time dependent distribution function will not work, because such systems do not evolve in time any further, even though the distributions are different from the equilibrium distributions.

Self-consistency
A point should be raised here regarding self-consistency -that is, using the density after the pump to determine the self-energy, rather than using the equilibrium density. This is equivalent to using the unrenormalized Green's function to evaluate the self-energy, which is often done in equilibrium approaches to evaluate first-order effects. As we will show, this is not the same in the time domain. Without self-consistency, the above approach cannot be applied in the same fashion because the relation between Σ and G k is not as clear. Instead, as above, let us assume that some time after the pump, the system can be described as a population occupying an equilibrium spectral function -that is, G < k (t ave , ω) = −2i f dist (ω, t ave )ImG R k,eq (ω). Then, if the self-energy is given by the original equilibrium one, we find where f eq (ω) is the equilibrium distribution used to evaluate self-energy. If the self-energy is thus determined based on the equilibrium Green's functions, there cannot be a balance of terms and the density will decay until it reaches the initial equilibrium state, independent of the types of interactions present. More generally, if the distribution is the same for G < k (t ave , ω) as for Σ < (ω), then there is no further relaxation.
One way this can be further understood is by breaking the self-energy into "dark" and "induced" pieces, following the work bySpicka et al. [14][15][16]. Starting from Eq. (18) we find At this point, we can identify the first two terms on the RHS as those captured by approaches lacking self-consistency -they capture the changes in the Green's function, but not the self-energy. If we limit the RHS to these terms and apply the approach discussed above, it becomes clear how the population dynamics have a simple connection to the equilibrium (or "dark") retarded self-energy: As long as the induced density is different from 0, we can observe a decay with an exponential decay constant τ(ω) −1 = −2ImΣ R D (ω). This is the relation between the population decay and self-energy reported previously [5,17]. It is the same relation as that between the decay of the Green's function in relative time t rel , but does not have the same origin. In that case, the self-energy gives rise to a line width in the spectrum of a single excitation, as measured e.g. by equilibrium photoemission spectroscopy. Here, the relation arrives only through copious approximations applied to population dynamics [14][15][16], which occur along an orthogonal direction to t rel . We require the population dynamics to be very slow for this to hold, or for the pump-induced population changes to be extremely small such that only these terms play a role. Once the approximations are violated, more complex dynamics arise, as discussed by several authors. [6,13,18,19] This can already be seen from Eq. (25), which bears a remarkable similarity to Eq. (27) except that the full pump-modified self-energy has to be accounted for. This leads to modification of the observed interactions by the pump, as discussed in. [6,20,21] Evidently, the lack of self-consistency is an uncontrolled approximation in the time domain, which can fail to capture important aspects of the population dynamics.

Specific examples
As a concrete illustration of the above, let us consider a simple form of impurity scattering which can account for the changes in the density due to the pump: the self-consistent Born approximation, where the self-energy is proportional to the local Green's function: , and similarly for the other Keldysh components. Inserting this into Eq. (18), and applying the equations above, we find This is a particularly striking result in light of the usual equilibrium connection between the quasiparticle lifetime and the self-energy. Here, the self-energy is clearly finite, whether calculated based on the pumped states or the equilibrium states, yet the population decay rate is 0 once the system has reached a state where the distribution function is independent of momentum. As a final note, we may apply the above methodology to a system of purely interacting electrons. The lowest-order term that leads to scattering of more than 1 electron is at 2 nd order, with Although the algebra is more complex, it can nevertheless be shown that once the electrons have been approximated as a thermal distribution occupying the spectral function, the RHS of the equations of motion are identically 0. This is as expected, since a thermalized state is the stable solution for a system of interacting electrons. Indeed, the specific form of the scattering mechanism is not important, the dependence of the lesser quantities on the distribution function is the important element.

Numerical evaluation of the scattering integrals
To evaluate the approximations and confirm the arguments made above, we perform numerical evaluation of the equations of motion on the Keldysh contour. The method has been described in some detail in Refs. [6,13]. We will follow the same separation used and motivated in the previous section, and consider separately the terms These are different from I 1 k and I 2 k , but rather are the full scattering integral separated into the terms that lead to two terms as they appear in Eq. (19). As noted above, these two terms should balance when there is no further time dependence to the density. The full scattering integral is then a difference between the two [see Eq. (18)].
As an illustrative model, we choose the Holstein-Hubbard Hamiltonian with local scattering: whereĉ † /ĉ andb † /b are the creation and annihilation operators for electrons and phonons, respectively, andn i ≡ĉ † iĉ i is the local density operator. The electron-electron interactions are treated at second order, with where the star ( ) denotes a product on the Keldysh contour. When present, U 2 = 0.09 eV 2 . The electron-phonon scattering is treated at first order within the Migdal-Eliashberg formalism, with where D C 0 (t, t ) is the propagator for a bare Einstein phonon with frequency Ω. When present, g 2 = 0.01 eV 2 and Ω = 0.2 eV. We drive with two-dimensional tight-binding model with nearest neighbor hopping V nn = 0.25 eV and a chemical potential µ = −0.2 eV. The field is applied using minimal coupling in the Hamiltonian gauge. We use a simple oscillatory field in the zone diagonal direction with A(t) = A max sin(ωt) exp(− 1 2 (t − t 0 ) 2 /σ 2 ), where ω = 0.5 eV, t 0 = 60 eV −1 , σ = 10 eV −1 and A max = 0.8.  Fig. 2 shows the scattering integrals for a case with both types of scattering (e-e and e-p) at t = 108 after the pump. The separated terms I Σ k and I G k are both large and appear featureless on these scales, while their sum has some definite structure. To elucidate this further, we consider the difference from equilibrium, ∆I Σ k and ∆I G k , shown in the bottom row of Fig. 2. The pump has clearly induced a difference from the equilibrium scattering integrals, giving an average time dependence to the self-energy and Green's functions and thus resulting in a finite summed scattering rate. There is some momentum dependence to the terms, although the details of the momentum dependence can vary based on the types of scattering present and their amplitudes. As time goes on, the system will reach some steady state where it has either returned to its equilibrium state (where ∆I Σ k and ∆I G k also go to 0) or where a new state is achieved, with a finite ∆I Σ k and ∆I G k . The latter can be achieved in the case of e-e scattering without any further mechanism for removing the energy from the system. These cases can be illustrated by considering the time dependence of ∆I Σ,G ≡ ∑ k ∆I Σ,G k . Fig. 3 shows the time traces of momentum-summed ∆I Σ and ∆I G . When only e-p scattering is present, the system behaves essentially as expected and reported previously [6] with an exponential, slow return to equilibrium together with a slow decrease of both scattering integrals. The sum also goes to zero on a similar time scale. It should be noted that the magnitude of ∆I is much larger than the two contributions ∆I Σ,G , confirming again that the balance between terms controls the dynamics. As e-e scattering is introduced, more energy is absorbed during the pump and the scale of the changes increases, but the salient features remain the same. Finally, when e-p scattering isn't present, a markedly different scenario occurs. The scattering integrals remain different from their equilibrium values, but as the system relaxes, they lose their average time dependence and balance out to give a net sum of 0 at long times. This state reflects a "thermalized" electron system, where the scattering has taken place until a balance (reflected in a Fermi distribution) is achieved. As noted above, the balance can also be achieved with a non-thermal electron distribution, so long as the same distribution is used in determining the self-energy.
To further illustrate the importance of cancellation between ∆I Σ and ∆I G , we can estimate the decay constant by considering the ratio In the weak pumping or non-selfconsistent limits, this is approximately (or exactly) equal to −2ImΣ(ω = k ). For e-p scattering only, ∆I k /∆n k shows a reasonable agreement with the equilibrium self-energy, although some deviation is present and in fact expected because we know that there is a contribution from the ∆I G term in the scattering integrals. Nevertheless, the step in the scattering rate near the phonon frequency is recognizable. Since the scattering integral ∆I does not show any signs of a balance before the full decay, we observe that ∆I k /∆n k remains relatively constant at early times. When e-e scattering is included, a competition between the scattering mechanisms arises that leads to a disagreement between 1/τ and −2ImΣ eq (ω = k )) [22]. This is reflected in the scattering integrals as a deviation from simple exponential decay, and a time dependence in ∆I k /∆n k . Although the phonon window can still be observed as the energy dissipation of the system is governed by the e-p process, the connection to the equilibrium self-energy is rapidly lost as the e-e scattering strives to achieve a thermal balance.

Discussion
We have considered a few aspects of the equations of motion for the population of an excited system. The major result is that a nontrivial average time dependence of the Green's function and self-energy is needed to have any population dynamics at all, and that a steady state is reached when there is a balance of scattering integrals after excitation rather than a full return to equilibrium of the same. Approximating the system with a hot electron model cannot capture the dynamics, as it was shown to lead to an exact cancellation of the scattering integrals unless the electrons are coupled to an external bath (e.g. of phonons). To be able to achieve this balance, a full self-consistent treatment is necessary. This naïvely makes sense because the redistributed population should be accounted for in the interactions, but stands in contrast to equilibrium physics where self-consistency often only gives a correction on top of the major feature obtained in 1 st order. For example, the kink and step in the spectral weight due to a strongly coupled boson in Migdal-Eliashberg theory arise readily ∆I k /∆n k , e-p only ∆t =120 eV −1 ∆t =180 eV −1 ∆t =240 eV −1 e-p −2ImΣ eq (ω) ∆I k /∆n k , e-e and e-p ∆t =120 eV −1 ∆t =140 eV −1 ∆t =160 eV −1 e-p −2ImΣ eq (ω) Figure 4. ∆I k /∆n k for times ∆t after t = 0 (the pump pulse occurs at t = 60 eV −1 ). The equilibrium electron-phonon self-energy is shown for reference.
whether the bare or renormalized electron propagator is used in the calculation of the self-energy. In the time domain, ignoring self-consistency does capture the main feature when the system is solely coupled to an external bath [5,6], but beyond this situation this approach is insufficient. This leads to the conclusion that the relation between the population scattering rate and the equilibrium self-energy τ(k, ω) −1 = −2ImΣ eq (k, ω) only holds exactly in very limited cases for non-equilibrium experiments, although it appears to be approximately true over a much wider range of cases.

Materials and Methods
The equations of motion for the Green's function are solved using the time-stepping algorithm detailed in Ref. [13]. The thermal problem is solved on the Matsubara axis, which is followed by a forward time-stepping procedure. At each time step, the self-energy is evaluated and the forward time step is made, followed by a correction step. This is repeated until the Green's function at the new time is converged, and the algorithm moves on to the next time step.