1.

INTRODUCTION

Back when many-body physics and Green’s function methods were being established, it was noticed that the self-energy provides the lifetime for the equilibrium Green’s function.¹ In particular, it is the imaginary part of the self-energy evaluated at the pole of the Green’s function that lies closest to the real axis in the lower half plane that determines this relaxation rate. This result leads one to infer that for small deviations from equilibrium, the imaginary part of the self energy should continue to provide the relaxation rate for the electrons. But trying to make this work immediately leads to a complication. Namely, in equilibrium there is only one time, while in nonequilibrium there are two times associated with the Green’s function. The decay rate determined from the linear-response analysis governs decay in the relative time direction, but the change in the population of the electrons in a given momentum state is governed by decay in the average time direction. While those decay rates might be related to one another, it is by no means obvious that they must be the same, and indeed, we often find they are not.

Hence, there is a need to re-evaluate how relaxation occurs in nonequilibrium, since it is not governed by the same behavior that drives lifetimes for equilibrium systems.

One thing that is always relevant for relaxation is the constraints provided by the phase space for scattering and the Pauli exclusion principle. This leads to the so-called “phonon window effect,”^{2, 3} where scattering by an optical phonon is sharply reduced as one gets to energies that lie below the phonon frequency. In particular, the relaxation outside the phonon window ω > Ω is quite rapid, but within the window ω < Ω, it becomes slow. As the system is excited, the Pauli blocking is reduced within the window and enhanced outside the window, so the relaxation rates move closer to each other, but, speaking quantitatively, they rarely become too similar.

These effects have already been seen in experiment. A direct comparison of the relaxation time for a population and the lifetime, as measured in linear-response angle-resolved photoemission spectroscopy (ARPES), showed large differences, sometimes more than an order of magnitude difference in the relaxation rates.⁴ The phonon window effect and its change with pump fluence have also been observed.⁷ Theory has started to examine these effects too. In particular, an equation of motion technique was used to determine the initial contributions to the relaxation rate at long times, and comparison with the imaginary part of the self-energy also showed significant differences, but not yet as large as seen in experiment.^{5, 6}

In this work, we approach the problem from a different perspective. It is well known from linear-response theory, that the relaxation of the current comes from the transport relaxation time, which is related to the imaginary part of the self-energy, but is clearly different as well. Here, we compute the generalization of that relaxation time to see whether it gives a relaxation rate that is closer to the relaxation rate that can be extracted directly from the time dependence of the populations. While we find some improvement, it is not significant, and it requires us to adjust an overall scaling factor for the rate, so it is not a complete determination of the relaxation rate. These results are intimately related to a breakdown of Mathiessen's rule. That rule says that when we have multiple scattering mechanisms, the relaxation rates for each mechanism add together to create a net relaxation rate. But in nonequilibrium relaxation, we often see multiple rates arise in different time ranges due to different relaxation processes and bottlenecks to energy transfer; we do not go into full detail of that phenomenon here, as we instead focus on trying to identify the primary relaxation mechanism for the electron-phonon interaction.

We work with the Holstein model, which involves a single band of uncorrelated electrons that interact with an optical phonon (Einstein mode) via a density-coordinate coupling. The Hamiltonian is

where c_iσ destroys (creates) an electron with spin σ at lattice site i, b_i are the phonon lowering (raising) operators for the optical phonon at site i, –t_ij is the hopping integral that connects site i with site j, µ is the electron chemical potential, g is the electron-phonon coupling, and Ω is the phonon frequency. The dimensionless electron-phonon coupling λ is given by the slope of the real part of the self-energy in equilibrium; here we have g² = 0.02 eV, Ω = 0.1 eV and λ ≈ 0.34. We work on a square lattice with only a nearest neighbour hopping of 0.25 eV. The system is at half-filling, which results when we choose µ = 0 and ignore the Hartree term in the perturbation theory (which we do because the Hartree term only shifts the chemical potential). The initial temperature is chosen to be T = 0.025 eV (room temperature). The electric field is described in the Hamiltonian gauge, where E(t) = –dA(t)/dt, where we set ℏ = c = 1. Then, we use the Peierls' substitution in the bandstructure to incorporate the field into the Hamiltonian. If we write

with k the momentum and R_j the position vector for the lattice site j. Then the Peierls substitution is ϵ(k) →ϵ(k - A(t)). We start the system off in equilibrium at a temperature T, and then turn on an electric field, which we assume to be spatially uniform. We ignore all magnetic field effects. The vector potential is a sinusoidal oscillating wave with a Gaussian envelope. The amplitude of the vector potential is 0.5, the standard deviation of the Gaussian is 10 eV^–1 and the oscillation frequency of the sinusoidal wave is 0.5 eV. This wave has 4-5 visible periods and is visibly nonzero over a total range of about 50 eV⁻¹.

The contour-ordered Green's functions depend on two times, each lying on the Kadanoff-Baym-Keldysh contour shown in Fig. 1. The local Green's function is defined via

Figure 1.

Kadanoff-Baym-Keldysh time contour, which runs from a minimum time to a maximum time along the real time axis, then backwards to the minimum time, and then parallel to the imaginary axis for a length given by the inverse of the initial equilibrium temperature.

where the angle brackets denote an average with respect to the initial equilibrium distribution

β = 1/T is the inverse temperature of the initial equilibrium distribution, and is the partition function for the initial equilibrium state. The time-dependence of the operators is expressed in the Heisenberg picture. The symbol T_c is the time-ordering operator on the contour, which places later objects, according to where they sit on the contour, to the left.

We employ Migdal-Eliashberg theory as the impurity solver for the nonequilibrium version of dynamical mean-field theory (DMFT)^{8, 9} to solve this problem by employing the Kadanoff-Baym-Keldysh formalism.^{10, 11} The strategy is briefly summarized in the next section, where we discuss the perturbation theory and the techniques employed to solve Dyson's equation. We also show the explicit formulas used for the data analysis.

2.

FORMALISM

Migdal-Eliashberg theory is employed to determine the local self-energy because the phonon energy scale is much smaller than the electron energy scale, implying we can neglect vertex corrections. Instead of how Migdal-Eliashberg theory is employed for linear-response and equilibrium calculations, where the phonons are the dressed phonons, so they are not renormalized, we self-consistently dress the phonons for the Holstein model by solving for the electron and phonon self-energies from a common conserving approximation. We simultaneously perform a self-consistent perturbation theory for the electronic Green's function; this approach allows us to take into account the finite heat capacity of the phonons and have their properties transiently change as they absorb energy. Since Migdal-Eliashberg theory involves a local self-energy, it is a form of DMFT, and we employ the NEDMFT approach to solving the problem. Note that all objects are contour-ordered continuous matrix operators, which depend on two times, each lying on the Kadanoff-Baym-Keldysh contour.

We start with a guess for the self-energy and then determine the local Green's function by solving the Dyson equation via

where the Green's function and self-energy are continuous matrix operators in time, and we take the (t, t′) matrix elements after the inverse. This equation is solved by discretizing it on the contour and employing the method of solving the equation of motion differential equation via direct integration.¹² When properly formulated, the algorithm for doing this is highly efficient. Eq. (5) includes the noninteracting nonequilibrium Green's function on the lattice, which can be found analytically.¹³ Once we have the local Green's function, we are ready to solve the impurity problem. In this case, we do not need to determine the effective medium and solve the full impurity problem because the expression for the self-energy depends only on the local Green's function, hence we have

which has to be solved self-consistently because the local Green's function depends on Σ^c, and the dressed phonon propagator depends on the Green's function. The dressed phonon propagator is constructed from the initial equilibrium propagator, which given by

Here, n_B (Ω) = 1/[exp(βΩ) – 1] is the Bose distribution function and θ_c (t, t′) is equal to one if t is ahead of t′ on the contour and is zero otherwise. We employ a conserving approximation where the electron and phonon self-energies are derived from a single functional. This results, in addition to the Dyson equation for the electrons, in one for the phonon:

The phonon self-energy Π^c (t, t′) is obtained from the electron Green's functions through:

where the product is to be evaluated through the Langreth rules. More details can be found elsewhere.^{2, 3}

To summarize how the Green's functions are calculated, we do the following: (i) we decide what the initial temperature T is of the system before the field is turned on; (ii) we incorporate the field via a spatially uniform vector potential that is oriented along the diagonal direction; and (iii) we iterate the nonequilibrium Migdal-Eliashberg theory (or equivalently, the NEDMFT) until both self-energies converge. At this stage, we have both the contour-ordered self-energy and the contour-ordered Green's function for both the electrons and the phonons. Note that we work with the full local Green's function here, summing over all momenta in the Brillouin zone.

We extract the data in a similar fashion to what would be done in experiment. We first construct the time-resolved photoemission spectroscopy signal, neglecting matrix element effects,¹⁴ which is formed from the lesser Green's function and the envelope of the probe pulse, , centered at t₀, with a spread (standard deviation) given by σ_pr (we use σ_pr = 25 eV^-1 in this work). The formula is

There is a small subtlety in what we do to determine the PES signal, and we do this to match results closer to how some experiments have been performed.¹⁵ Namely, we first construct the gauge-invariant Green's function and then we perform a partial summation over momentum, by summing only over the diagonal direction for k. This is then used to approximate the TR-PES signal, which should be a good approximation because any anisotropy is washed out by the electron-phonon scattering, as in low-temperature superconductors. After computing the TR-PES signal, the frequency axis is divided into contiguous bins (with a width of 0.01) and we integrate the total signal within each bin and plot as a function of the probe time t₀. These populations then decay as a function of time, with a typical exponential decay. To find the exponent, we either fit the tail of the curve to an exponentially decaying curve, or we extract the time-constant for the decay directly by numerically calculating the derivative and dividing the derivative by the function to give the decay rate (under the assumption that the system is decaying exponentially). As we will see below, both methods give similar results for the decay rate at a given time.

We end this section with a discussion about linear response. Within DMFT, the optical conductivity has no vertex corrections for the linear-response regime.¹⁶ Hence, we can evaluate the dc conductivity in terms of a many-body transport relaxation time^{17, 18}

where f (ω) = 1/[1 + exp(βω)] is the Fermi-Dirac distribution function and in two-dimensions, we have

where the R superscript means retarded, and the transport Green's function is defined with an extra v² (square of the band velocity) in its definition, so that

Here v = ∇_k ϵ(k) is the band velocity.

Unfortunately, this formula does not easily generalize to the nonequilibrium limit. But in the spirit of this result, we examine a relaxation time similarly constructed from the retarded transport Green's function and the retarded self-energy at the same average time; we neglect the second and third terms in Eq. (12) since this is likely to only be semiquantitative result at best, and we assume the derivative of the Fermi-Dirac distribution is so sharp it can be approximated by a delta function for each population at a given energy above the chemical potential. Our ansatz is that we compare

against the calculated decay rates extracted from analyzing the photoemission spectra, with an overall normalization factor still to be determined. We call this transport-based relaxation rate η.

3.

DISCUSSION AND RESULTS

In Fig. 2, we plot the time-resolved photoemission spectra for different frequency bins as a function of time. The grayed region is the region where the pump is on. We take the data for different frequencies and times and extract an effective relaxation time in one of two different ways. The first way is to fit the data according to an exponential fit for the data in a time window about the given time. The second method is to extract an effective exponential relaxation time by computing the numerical derivative and dividing by the function at a given time; that is, by computing the logarithmic derivative. This satisfies

Figure 2.

(Color online.) Time-resolved photoemission spectra for different frequency bins as a function of time. One can see that the higher frequencies decay much faster than the lower frequencies due to the phonon window effect (here, Ω = 0.1 eV). The grayed region is where the field is on, which has a total width of about 50 eV⁻¹.

The results are plotted in Fig. 3. One can see that both methods for extracting the instantaneous decay rate agree to high accuracy, but the overall relaxation rate differs from the equilibrium self-energy result.

Figure 3.

(Color online.) Relaxation rate for different frequency bins at t_ave = 50 eV⁻¹. The imaginary part of the retarded equilibrium self-energy is shown for comparison. While there is semiquantitative agreement here, they clearly differ.

We next compare these relaxation rates, extracted with the two different methods, for different average times in Fig. 4. We can see a phonon window effect for short times, which disappears at longer times, and we can see a generic behavior that produces a fairly flat response for the relaxation rate over these frequency values.

Figure 4.

(Color online.) Left: relaxation rate as determined from a fit to an exponential decay for different average times. The phonon frequency is marked by a dashed line to help identify the phonon window effect. Right: a similar plot for the relaxation rates as extracted from a logarithmic derivative.

Finally, we form the result from the transport Green's function and compare with the relaxation rate data in Fig. 5. The overall factor is about 6. One can immediately conclude that this transport Green's function-based relaxation time is much flatter in frequency, in agreement with the data, but it does not show the strong average time-dependence within the phonon window, where the phonon window effect disappears for long average times.

Figure 5.

(Color online.) Transport-based relaxation rate (renormalized for best fit) compared to the relaxation rate extracted from the TR-PES data. One can see it is much flatter outside the phonon window, in good agreement with the data, but it does not show a disappearance of the phonon window effect at long times like the data do (compare the blue lines for the longest times).

4.

CONCLUSIONS

Nonequilibrium relaxation is complicated. Indeed, the simple notions for how to determine it from the self-energy are known to fail when the electron correlations become strong. Here, we have made an initial attempt to remedy this problem by considering the modification of the relaxation rate due to a generalization of the transport relaxation time to nonequilibrium. We find that while this approach does do better in modeling the weak frequency dependence of the data outside the phonon window, it does not properly show the evolution of this dependence inside the window, especially for long times. So, this might be a step in the right direction, but, it unfortunately has an adjustable parameter for the overall relaxation rate, which needs to be determined for this system, and which provides less predictive power than if we had a prediction on an absolute scale.

In the future, we hope to be able to find an even better ansatz for the nonequilibrium relaxation time in order to find the microscopic origin of this relaxation. The nonequilibrium behavior is far more complex than the linear-response regime and we need more work to both have good data to compare with and to determine the proper microscopic basis for the relaxation, including all bottleneck effects.