2.1. Equations of Motion
The analysis of the two-time Green’s function is done on the Keldysh contour, illustrated in
Figure 1. The imaginary spur allows for the calculation of an equilibrium thermal state at
, which is subsequently propagated forward in time along the contour using the equations of motion discussed below. For convenience, we will assume
. The system is then subjected to an external perturbation (a field) at some time
.
The Green’s function for a state with quasi-momentum
is defined on the contour as
where
is the contour time ordering operator, and
are the usual creation and annihilation operators. The superscript
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
. 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
and
. Finally, the component which has both arguments on the imaginary spur is the usual Matsubara Green’s function
. 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
direction. This leads to a decay in the contribution of the mixed pieces to the integrals as time gets further from
. We will assume that
t and
have advanced sufficiently far from
that we can neglect those mixed pieces in our analysis.
The density for momentum
is given by the time-diagonal piece of the lesser Green’s function,
. To propagate along the
(
) direction, we need to combine Equations (
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
commutes with the Green’s function. We can further assume that we have gone far enough along in
t and
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.
2.2. 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 : . The population dynamics occurs along the direction, whereas quasiparticle lifetimes arise from the direction. The scattering integrals control the temporal dynamics of the Green’s function along the 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 , 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
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
that we can replace the average time dependence with a single time only, namely
. 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
, we can use Dyson’s equation and substitute for
and
. Long after
where the mixed components have decayed, we can re-write Dyson’s equation (Equations (
8) and (
12)) in the following way:
Making the same transformation to
and Fourier transforming over
, we can substitute this result into Equation (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
(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
, where the fluctuation-dissipation theorem holds and
is the Fermi function. We can readily see this from Equation (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
and
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.
2.2.1. 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
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,
. Then, if the self-energy is given by the original
equilibrium one, we find
where
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
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 by S̆pic̆ka
et al. [
14,
15,
16]. Starting from Equation (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
. 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 , 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
. 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 Equation (22), which bears a remarkable similarity to Equation (24) 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 such as predicting the wrong decay rates and a return to equilibrium when the scattering process cannot dissipate energy.
2.2.2. 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 Equation (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 2nd 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.
The conventional “hot-electron” viewpoint, where the electrons rapidly relax to a Fermi-Dirac distribution at high temperature which then cools down to a low temperature Fermi-Dirac distribution, cannot take place in this fashion. Once the system is described by a quasi-equilibrium distribution function, it will no longer relax.
2.3. 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
and
. These are different from
and
, but rather are the full scattering integral separated into the terms that lead to two terms as they appear in Equation (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 Equation (18)).
As an illustrative model, we choose the Holstein-Hubbard Hamiltonian with local scattering:
where
and
are the creation and annihilation operators for electrons and phonons, respectively, and
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,
eV
. The electron-phonon scattering is treated at first order within the Migdal-Eliashberg formalism, with
where
is the propagator for a bare Einstein phonon with frequency Ω. When present,
eV
and
eV. We employ a two-dimensional tight-binding model with nearest neighbor hopping
eV and a chemical potential
eV. The field is applied using minimal coupling in the Hamiltonian gauge. We use a simple oscillatory field in the zone diagonal (
) direction with
, where
eV,
eV
,
eV
and
.
Figure 2 shows the scattering integrals for a case with both types of scattering (e-e and e-p) at
after the pump. The separated terms
and
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,
and
, shown in the bottom row of
Figure 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
and
also go to 0) or where a new state is achieved, with a finite
and
which balance. In the first case, the system is once again thermal (since it started from a thermal state). In the latter case, the system does not have to achieve a thermal state—it is sufficient that
and
are described by the same population distribution. A thermal state also satisfies this, but nonthermal states can as well. In the case we are studying with phonons, the system will approach a thermal state at the temperature of the phonon bath, but other situations may result in nonthermal states since both are naturally possible within the general theory.
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
.
Figure 3 shows the time traces of momentum-summed
and
. 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
is much smaller than the two contributions
, 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 is not 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. In this particular case, the state reflects a “thermalized” electron system, where the scattering has taken place until a balance (reflected in a Fermi distribution) is achieved. However, as noted above, the balance is not limited to thermalized states. It 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
and
, we can estimate the decay constant by considering the ratio
In the weak pumping limit this is approximately equal to
. The approximate relation becomes exact if self-consistency is neglected. For e-p scattering only,
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
term in the scattering integrals. Nevertheless, the step in the scattering rate near the phonon frequency is recognizable. Since the scattering integral
does not show any signs of a balance before the full decay, we observe that
remains relatively constant at early times. When e-e scattering is included, a competition between the scattering mechanisms arises that leads to a mismatch between
and
[
22]. This is reflected in the scattering integrals as a deviation from simple exponential decay, and a time dependence in
. Although the phonon window can still be observed as the energy dissipation of the system is governed by the e-p process, the exact connection to the equilibrium self-energy is rapidly lost as the e-e scattering modifies the relaxation of the system so that the electrons only thermalize when they reach the phonon temperature.