An introduction to many-body Green's functions in and out of equilibrium
J. K. Freericks
cond-mat.str-el





§ INTRODUCTION


Green's function
This chapter is all about Green's functions. Many-body Green's functions are a marvelous tool to employ in quantum-mechanical calculations. They require a rather large set of mathematical machinery to develop their theory and a sophisticated computational methodology to determine them for anything but the simplest systems. But what do we really use them for?  Essentially, there are two main uses: one is to compute the thermodynamic expectation value of c^† c, which allows us to determine all interesting single-particle expectation values, like the kinetic energy, the momentum distribution, etc. and two is to determine the many-body density of states, which tells us how the quantum states are distributed in momentum and energy. In large dimensions, they can also be used to determine some two-particle correlation functions like the optical conductivity, because the vertex corrections vanish. It seems like this is a lot of work to end up with only a few concrete results. But that is the way it is. No one has figured out any way to do this more efficiently. If you determine how to, fame and fortune are likely to follow!
Green's function

Green's function!history
Before jumping into the full glory, we will give just a brief history. This is one that focuses on my opinion of where critical ideas originated. It is not intended to be exhaustive or complete, and of course I may be wrong about where the different ideas came from. The equilibrium theory for Green's functions was developed primarily in the 1950s. The Lehmann representation <cit.> was discovered in 1954. The Russian school developed much of the perturbative approach to Green's functions, which is summarized in the monograph of Abrikosov, Gorkov, and Dzyaloshinski <cit.>. Joaquin Luttinger and John Ward developed a functional approach <cit.> and Matsubara determined how to solve for thermal Green's functions <cit.>. But the reference I like the most for clarifying many-body Green's functions both in and out of equilibrium is the monograph by Baym and Kadanoff <cit.>. Their approach is one we will follow, at least in spirit, in this chapter. Of course, Keldysh's perpsective <cit.> was also important.  Serious numerical calculations of Green's functions (in real time and frequency) only began with the development of dynamical mean-field theory in the 1990s <cit.>. The nonequilibrium generalization began only in the mid 2000's <cit.>, but it was heavily influenced by earlier work from Uwe Brandt <cit.>.
Green's function!history

Now that we have set the stage for you as to where we are going, get ready for the ride. The theory is  beautiful, logical, abstract, and complex. Mastering it is a key in becoming a many-body theorist. We begin in section 2 with a discussion of equilibrium Green's functions focused on the Lehmann representation. Section 3, generalizes and unifies the formalism to the contour-ordered Green's function. Section 4 introduces the self-energy and the equation of motion. Section 5 illustrates how to include the external electric field. Section 6 introduces how one solves these problems within the nonequilibrium dynamical mean-field theory approach. Numerics are discussed in section 7, followed by examples in section 8. We conclude in section 9.




§ GREEN'S FUNCTIONS IN EQUILIBRIUM AND THE LEHMANN REPRESENTATION


Green's function!equilibrium
Green's function!retarded
Green's function!advanced
Green's function!lesser
Green's function!greater
We begin with the definition of four Green's functions in equilibrium: the retarded G_ijσ^R(t), advanced G_ijσ^A(t), lesser G_ijσ^<(t) and greater G_ijσ^>(t). These Green's functions are time-dependent complex thermal expectation values weighted by the thermal density matrix exp(-βℋ)/𝒵, with 𝒵= Trexp(-βℋ) the partition function. They are defined in terms of the (anticommuting)  fermionic creation (annihilation) operators for an electron at site i with spin σ: c_iσ^† (c_iσ^†). Time dependence is handled in the Heisenberg picture, where the operators depend on time and the states do not, given by 𝒪(t)=U^†(t)𝒪U(t). For a time-independent Hamiltonian, the time-evolution operator U(t) is simple and expressed as U(t)=exp(-iℋt), where we set ħ=1. The four Green's functions then become

    G_ijσ^R(t)   =   -iθ(t)1/𝒵 Tre^-βℋ{c_iσ^†(t),c_jσ^†(0)}_+,
    
    G_ijσ^A(t)   =   iθ(-t)1/𝒵 Tre^-βℋ{c_iσ^†(t),c_jσ^†(0)}_+,
    
    G_ijσ^<(t)   =   i1/𝒵 Tre^-βℋc_jσ^†(0)c_iσ^†(t),
    
    G_ijσ^>(t)   =   -i1/𝒵 Tre^-βℋc_iσ^†(t)c_jσ^†(0).

Here, θ(t) is the Heaviside unit step function, which vanishes for t<0, is equal to 1 for t>0 and we choose it to be 1/2 for t=0, and the symbol {·,·}_+ denotes the anticommutator. Note that the fermionic operators satisfy the canonical anticommutation relations, namely {c_iσ^†,c_jσ^'^†}_+=δ_ijδ_σσ^'.
These Green's functions can also be defined in momentum space, with the obvious changes of the subscripts from real space to momentum space—the only significant modification is that the Green's functions will be diagonal in momentum for translationally invariant systems—which is all we will study in this chapter.
Green's function!equilibrium
Green's function!retarded
Green's function!advanced
Green's function!lesser
Green's function!greater

Green's function!momentum-dependent
Green's function!retarded
Lehmann representation
density of states
Heisenberg representation
These definitions may look like they come out of the blue, but the key point to realize is that when t=0, the lesser Green's function provides precisely the expectation value we need to be able to calculate any single-particle expectation value. We will see in a moment that the retarded Green's function also leads to the momentum-dependent spectral function and the density of states. The best way to understand and make sense of these Green's functions is with the Lehmann representation <cit.>. This allows us to explicitly determine the Green's functions as functions of frequency via a Fourier transformation:

    G_ijσ(ω)=∫_-∞^∞ dt e^iω tG_ijσ(t).

The trick is to insert the resolution of unity (𝕀=∑_n|n⟩⟨ n| where |n⟩ is an eigenstate of ℋ such that ℋ|n⟩=E_n|n⟩) in between the creation and annihilation operators in the definitions of the Green's functions. We illustrate how this works explicitly for the retarded Green's function, where we have

    G_ijσ^R(t)=-iθ(t)1/𝒵∑_mn [ ⟨ m|e^-βℋc_iσ^†(t)|n⟩⟨ n|c_jσ^† |m⟩+⟨ m| e^-βℋc_jσ|n⟩⟨ n|c_iσ^†(t)|n⟩ ].

We now use the fact that these states are eigenstates of ℋ and the Heisenberg representation for the time dependent operators to find

    G_ijσ^R(t)=-iθ(t)∑_mne^-β E_m/𝒵 [ e^i(E_m-E_n)t⟨ m| c_iσ^†|n⟩⟨ n|c_jσ^†|m⟩+
    e^-i(E_m-E_n)t⟨ m|c_jσ^†|n⟩⟨ n|c_iσ^†|m⟩ ].

Next, we interchange m↔ n in the second term to give

    G_ijσ^R(t)=-iθ(t)1/𝒵∑_mn ( e^-β E_m+e^-β E_n )
    e^i(E_m-E_n)t⟨ m |c_iσ^†|n⟩⟨ n|c_jσ^†|m⟩.

It is now time to do the Fourier transform. We achieve this by shifting ω→ω+i0^+ to make the integrand vanish at the large time limit, which finally results in

    G_ijσ^R(ω)=1/𝒵∑_mne^-β E_m+e^-β E_n/ω+E_m-E_n+i0^+⟨ m |c_iσ^†|n⟩⟨ n|c_jσ^†|m⟩.

This is already of the form that we can discover some interesting things. If we recall the Dirac relation 1/ω+i0^+=𝒫/ω-iπδ(ω) (with 𝒫 denoting the principal value), then we find the local density of states via

    A_iσ(ω)=-1/π ImG_iiσ(ω)=∑_mnδ(ω+E_m-E_n) e^-β E_m+e^-β E_n/𝒵|⟨ m|c_iσ|n⟩|^2.

This density of states satisfies a sum rule. In particular,

    ∫_-∞^∞ dω A_iσ(ω)=1/𝒵∑_n e^-β E_n ( c_iσ^†c_iσ^†+c_iσ^† c_iσ^† )=1,

which follows from the anticommutation relation of the fermionic operators; note that you can also see this directly by taking the limit t→ 0^+ in Eq. (<ref>). This result holds in momentum space too:
the retarded Green's function becomes

    G_kσ^R(t)   =   -iθ(t)1/𝒵 Tre^-βℋ{ c_kσ^†(t),c_kσ^†(t)}_+
    
    G_kσ^R(ω)   =   1/𝒵∑_mne^-β E_m+e^-β E_n/ω+E_m-E_n+i0^+⟨ m|c_kσ^†|n⟩⟨ n|c_kσ^†|m⟩,

with the operators being the momentum space operators, determined by a Fourier transform from real space:

    c_kσ^†=1/|Λ|∑_j∈Λc_jσ^†e^-i k· R_j     and    
    c_kσ^†=1/|Λ|∑_j∈Λc_jσ^†e^i k· R_j.

Here, the lattice with periodic boundary conditions is denoted by Λ={ R_i} and |Λ| is the number of lattice sites in the system. We usually do not use a bold font for the momentum label in the subscript of the creation or annihilation operator; one should be able to figure out from the context whether we are working in real space or momentum space.
The momentum-space Green's functions have a similarly defined spectral function A_kσ(k,ω)=- ImG^R_kσ(ω)/π, which also satisfies a sum rule given by the integral over all frequency being equal to one. One can easily show that the advanced Green's function in frequency space is equal to the Hermitian conjugate of the retarded Green's function, because we need to shift ω→ω-i0^+ to control the integrand at the lower limit. Hence, we have G^A_ijσ(ω)=G^R*_jiσ(ω).
Green's function!momentum-dependent
Green's function!retarded
Lehmann representation
density of states
Heisenberg representation

Green's function!momentum-dependent
Green's function!lesser
Green's function!greater
Lehmann representation
spectral function
fluctuation-dissipation theorem
We also sketch how the calculation works for the lesser and greater Green's functions. Here, we have to control the integrand at both endpoints of the integral. To do so, we split it at 0  and introduce the appropriate infinitesimal convergence factor in each piece. The rest of the calculation proceeds as above and one finds that the two pieces have their imaginary parts add, so we finally obtain:

    G_ijσ^<(ω)=2iπ1/𝒵∑_mnδ(ω+E_m-E_n)e^-β E_n⟨ m|c_iσ^†|n⟩⟨ n|c^†_jσ|m⟩.

The greater Green's function similarly becomes

    G_ijσ^>(ω)=-2iπ1/𝒵∑_mnδ(ω+E_m-E_n)e^-β E_m⟨ m|c_iσ^†|n⟩⟨ n|c^†_jσ|m⟩.

Note that we can use the delta function to immediately infer an important identity, namely that 
G^<_ijσ(ω)=-exp(-βω)G^>_ijσ(ω). Combining this result with the fact that ImG_ijσ^R(ω)= Im[G_ijσ^<(ω)-G_ijσ^>(ω)]/2, we finally learn that

    G_ijσ^<(ω)=-2 if(ω) ImG_ijσ^R(ω)    and   
    G_ijσ^>(ω)=2 i[1-f(ω)] ImG_ijσ^R(ω).

This equilibrium relationship is sometimes called the fluctuation-dissipation theorem.
Green's function!momentum-dependent
Green's function!lesser
Green's function!greater
Lehmann representation
spectral function
fluctuation-dissipation theorem



§ GREEN'S FUNCTIONS OUT OF EQUILIBRIUM AND THE “CONTOUR”


-1 in


Kadanoff-Baym-Keldysh contour
Green's function!time representation
Now we move onto nonequilibrium where the Hamiltonian is time-dependent ℋ(t). The evolution operator depends on two times and becomes U(t,t^')=𝒯_t exp [-i∫_t^'^t dt̅ℋ(t̅) ]. It satisfies the semigroup property U(t,t^'')U(t^'',t^')=U(t,t^'), with U(t,t)=𝕀 and i∂_t U(t,t^')=ℋ(t)U(t,t^'). The Green's functions in nonequilibrium are defined by the same equations as before, namely Eqs. (<ref>–<ref>), except now they depend on two times, with t^' being the argument of the c^† operator. All of these Green's functions require us to evaluate one of two matrix elements. We examine one of them here and assume t>t^' for concreteness:

    1/𝒵 Tre^-βℋc_iσ^†(t) c_jσ^†(t^')    -0.1in =    -0.1in
    1/𝒵 Tre^-βℋU^†(t,-∞)c_iσ^†U(t,-∞)U^†(t^',-∞)c_jσ^† U(t^',-∞)
        -0.1in =    -0.1in 1/𝒵 Tre^-βℋU(-∞,t)c_iσ^†U(t,t^')c_jσ^† U(t^',-∞)

where we used the facts that U^†(t_1,t_2)=U(t_2,t_1) and the semigroup property to show that U(t,-∞)U^†(t^',-∞)=U(t,t^'). The time-evolution operators, including time evolution in the imaginary axis direction for the exp(-βℋ) term, can be seen to all live on the so-called Kadanoff-Baym-Keldysh contour, which is shown in Fig. <ref>. Starting with time at -∞, we evolve forward to t^', then apply the c^†, evolve forward to time t, apply the c operator, and then evolve backwards in time to -∞. The lesser and greater Green's functions are easily determined in this fashion. The retarded or advanced Green's functions can then be found by taking the appropriate differences of greater and lesser functions        with the convention that t>t^' for the retarded case and t<t^' for the advanced case. Now that we have these Green's functions, we can generalize the definition of the Green's function to allow both times to lie anywhere on the contour.  We take the time-ordered Green's function, with the time-ordering taking place along the contour, and we call it the contour-ordered Green's function:

    G^c_ijσ(t,t^')={
    -i1/𝒵 Tr e^-βℋ c_iσ^†(t)c_jσ^†(t^')   for t>_ct^'
     
    
    i1/𝒵 Tr e^-βℋ c_jσ^†(t^')c_iσ^†(t)      for t<_ct^'.
     .,

where the c subscript on the less than or greater than symbol is to denote whether one time is ahead of or behind the other on the contour, regardless of the numerical values of t and t^'.
This Green's function is the workhorse of nonequilibrium many-body physics—we use it to calculate anything that can be calculated with the single-particle Green's functions.
Kadanoff-Baym-Keldysh contour
Green's function!time representation

Note that in these lectures, we will work with the contour-ordered Green's function itself, which is a continuous matrix operator in time. In many other works, the Green's function is further decomposed into components where the times lie on the same or different branches of the real contour (2× 2 matrix) <cit.> or of the two real and one imaginary branches of the contour (3× 3 matrix) <cit.>. In the 3× 3 case, one has real, imaginary and mixed Green's functions, the last ones determining the so-called initial correlations. In the 2× 2 case, it is common to transform the matrix to make it upper triangular, expressing the retarded, advanced and so-called Keldysh Green's functions. In general, I find that these decompositions and transformations make the material more difficult than it  actually is and it is best to defer using them until one has gained some mastery of the basics. Here, we will work solely with the contour-ordered functions.



§ THE SELF-ENERGY AND THE EQUATION OF MOTION


Green's function!contour-ordered
equation of motion
self-energy
The next step in our analysis is to derive the equation of motion. We start from the definition of the contour-ordered Green's function in Eq. (<ref>) and we employ a unit step function on the contour θ_c(t,t^') which is one for t ahead of t^' on the contour and zero for t behind t^'. This theta function can be used to impose the two conditions. Its derivative with respect to time becomes a delta function on the contour ∂_tθ_c(t,t^')=δ_c(t,t^'). The delta function on the contour satisfies 

    ∫_c dt δ(t,t^')f(t^')=f(t),

for any time t on the contour (including those on the imaginary spur). Now a simple derivative of the Green's function with respect to time yields the equation of motion (we find it convenient to illustrate this in momentum space)

    i∂/∂ t G_kσ^c(t,t^')=1/𝒵 Tre^-βℋ{c_kσ^†(t),c_kσ^†(t)}δ_c(t,t^')
       +   iθ_c(t,t^')1/𝒵 Tre^-βℋ [ℋ(t),c_kσ^†(t) ]c_kσ^†(t^')-iθ_c(t^',t)1/𝒵 Tre^-βℋc_kσ^†(t^') [ℋ(t),c_kσ^†(t) ]
       =   δ_c(t,t^')+[ϵ_k(t)-μ]G_kσ^c(t,t^')
    +∫_cdt^''Σ_c(t,t^'')G_kσ^c(t^'',t^').

Here, we assumed that the Hamiltonian (in the Schrödinger picture) takes the form

    ℋ_s(t)=∑_kσ[ϵ_k(t)-μ]c_kσ^† c_kσ^†+V(t)

where ϵ_k(t) is a time-dependent bandstructure and μ is the chemical potential. The symbol V(t) is the potential energy operator of the Hamiltonian. The commutator of the potential with the fermionic destruction operator is complicated—we introduce the contour-ordered self-energy Σ_c to describe that cumbersome term—which can be viewed as the definition of the self-energy. The self-energy appears as a convolution with the Green's function (over the entire contour) in Eq. (<ref>). The equation of motion can be rearranged to look like

    ∫_c dt^''{ [ i∂/∂ t-ϵ_k(t)+μ ]δ_c(t,t^'')+Σ_c(t,t^'')}G_kσ^c(t^'',t^')=δ_c(t,t^').

As we mentioned before, the Green's functions are continuous matrix operators. Seen in that light, it should be clear that the object inside the curly brackets is the matrix inverse of the Green's function (since matrix multiplication, implied by the integration over the contour, has the matrix defined by the curly brackets multiplied by G^c and equaling the identity matrix).
Green's function!contour-ordered
equation of motion
self-energy

At this stage it may look like we are pulling the wool over your eyes. In some sense we are. We have hidden the difficult question of how to evaluate the commutator of the potential with the destruction operator, by simply introducing a self-energy which makes the complications go away. This really would lead us nowhere if we did not have another means by which we can determine the self-energy. Fortunately we do, and it will be introduced below. You will have to wait to see how that fits in. All of the challenge in many-body physics boils down to this question of “how can we calculate the self-energy?”

Now is also a good time to discuss the issue of how do we work with these complicated objects? If you look carefully, you will see we are working with matrices. Admittedly they are continuous matrix operator equations, but like any continuous object, we discretize them when we want to work with them on a computer. And indeed, that is what we will do here when we describe the DMFT algorithm below. We always extrapolate the discretization to zero to recover the continuous matrix operator results. It is important to recognize these facts now, as it will make your comfort level with the approach much higher as we develop the subject further.



§ INCORPORATING AN EXTERNAL FIELD AND THE ISSUE OF GAUGE INVARIANCE


I always found it odd in dealing with fields in quantum mechanics to learn that we do not always respect Maxwell's equations. This should not come as a surprise to you, because Maxwell's equations are relativistically invariant, and quantum mechanics only becomes relativistically invariant with a properly constructed relativistic quantum field theory. Nevertheless, we choose to neglect even more aspects of Maxwell's equations than just relativistic invariance. This becomes necessary because the quantum problem cannot be solved unless we do this (and, fortunately, the neglected terms are small).

The main issue we face is that we want the field to be spatially uniform but time varying. Maxwell says such fields do not exist, but we press on anyway. You see, a spatially uniform field means that I still have translation invariance and that is needed to make the analysis simple enough that it can be properly carried out. We justify this in part because optical wavelengths are much larger than lattice spacings and because magnetic field effects are small due the extra 1/c factors multiplying them.
Even with this simplification, properly understanding how to introduce a large amplitude electric field is no  easy task. 

Peierls' substitution
The field will be described via the Peierls' substitution <cit.>. This is accomplished by shifting k→ k-eA(t), which becomes k+eEt, for a constant (dc) field. This seems like the right thing to do, because it parallels the minimal substitution performed in single-particle quantum mechanics when incorporating an electric field, but we have to remember that we are now working with a system projected to the lowest lying s-orbital band, so it may not be obvious that this remains the correct thing to do. Unfortunately, there is no rigorous proof that this result is correct, but one can show that if we work in a single band (neglecting interband transitions) and have inversion symmetry in the lattice, then the field generated by a scalar potential, and the field generated by the equivalent vector potential will yield identical noninteracting Green's functions <cit.>. We then assume that adding interactions to the system does not change this correspondence, so we have an approach that is exact for a single-band model. One should bear in mind that there may be subtle issues with multiple bands and the question of how to precisely work with gauge-invariant objects in mutliband models has not been completely resolved yet. Fortunately, for model system calculations, working with a single band is often sufficient to illustrate the nonquilibrium phenomena. This is what we do here.
Peierls' substitution

Falicov-Kimball model
Peierls' substitution
It is now time to get concrete in describing the Hamiltonian we will be working with, which is the spinless Falicov-Kimball model <cit.>. This model is nice to work with because it does possess strongly correlated physics (Mott transition and order-disorder transitions), but it can be solved exactly within DMFT in both equilibrium and nonequilibrium. Hence, it provides a nice test case to illustrate the approach we take here. The Hamiltonian, in an electric field is

    ℋ_S(t)=∑_k [ ϵ(k-eA(t))-μ ]c_k^† c_k^†
    +E_f∑_i f_i^† f_i^†+U∑_i c_i^† c_i^† f_i^† f_i^†,

where we employed the Peierls' substitution in the bandstructure. Here, the conduction electrons (c) are spinless and do not interact with themselves. They interact instead with the localized electrons (f), when both electrons are situated on the same lattice site—E_f is the local site energy of the f-electrons. Note that we use k for momentum and i for position in ℋ. All of the time dependence is in the kinetic-energy term, because the field does not couple to the localized electrons (due to the fact that they do not hop). Hence, we have an effective single-band model and the discussion above about the applicability of the Peierls' substitution applies.
Falicov-Kimball model
Peierls' substitution

Green's function!nonequilibrium!noninteracting
It is instructive at this stage to determine the noninteracting Green's function exactly. We do this by determining the Heisenberg creation and annihilation operators from which we form the Green's function. In the case where U=0, one can immediately find that that the destruction operator (in the Heisenberg picture) becomes

    c_kσ^†(t)=e^i∫_-∞^t dt̅ [ϵ(k-eA(t̅)-μ]c_kσ^† c_kσ^†
    c_kσe^-i∫_-∞^t dt̅ [ϵ(k-eA(t̅)-μ]c_kσ^† c_kσ^†
    =e^-i∫_-∞^t dt̅ [ϵ(k-eA(t̅))-μ]c_kσ^†.

There is no need for time ordering here, because the Hamiltonian commutes with itself at different times. This simplification is what allows one to solve the problem easily.  By taking the Hermitian conjugate of Eq. (<ref>), one can substitute in the definitions of the Green's functions and immediately find that

    G_kσ^<(t,t^')=i1/𝒵 Tre^-βℋe^-i∫_t^'^t dt̅ [ ϵ(k-eA(t̅))-μ ]c_kσ^† c_kσ^†=ie^-i∫_t^'^t dt̅ [ ϵ(k-eA(t̅))-μ ]⟨ n_kσ⟩,

where ⟨ n_kσ⟩= Trexp[-βℋ]c_kσ^† c_kσ^† is the initial equilibrium momentum distribution. We similarly find that the noninteracting retarded Green's function in nonequilibrium is

    G_kσ^R(t.t^')=-iθ(t-t^')e^-i∫_t^'^t dt̅ [ ϵ(k-eA(t̅))-μ ];

there is no equilibrium expectation value here because the anticommutator of the equal-time fermionic operators is 1. To find the corresponding local Green's functions, we simply sum these expressions over all momenta.
Green's function!nonequilibrium!noninteracting

gauge invariance
Green's function!gauge-invariant
We next discuss the issue of gauge invariance. This issue is important, because we are required to perform calculations in a specific gauge but any physically measurable quantity must be gauge-invariant. So how does this come about? To begin, one should notice that the Peierls' substitution guarantees that any quantity that is local, like the current, is gauge invariant, because the sum over all momenta includes every term, even if the momenta are shifted due to the Peierls substitution. Momentum-dependent quantities, however, appear to depend directly on the vector potential. To properly formulate a gauge-invariant theory requires us to go back and formulate the problem (using a complete basis set, not a single-band downfolding) properly to illustrate the gauge-invariant nature of the observables. No one has yet been able to do this for a system with a nontrivial bandstructure. For a single-band model, we instead directly enforce gauge invariance using the so-called gauge-invariant Green's functions, introduced by Bertoncini and Jauho <cit.>. These Green's functions are constructed to be manifestly gauge invariant—the change in the phase of the creation and destruction operators is chosen to precisely cancel the change in the phase from the gauge transformation. We determine them in momentum space by making a shift

    k→ k(t,t^')=k+∫_-1/2^1/2dλ  A ( t+t^'/2+λ(t-t^') ).

Note that one needs to be careful in evaluating this integral in the limit as t→ t^'. In this limit, we find k→ k+A(t).
The explicit formula for the gauge-invariant Green's function, which uses a tilde to denote that it is gauge-invariant,  is

    G̃_kσ(t,t^')=G_k(t,t^')σ(t,t^').

The demonstration of the gauge-invariance is given in Ref. <cit.>. Note that the local Greens functions are always manifestly gauge-invariant because the map k→ k(t,t^') is a “one-to-one and onto” map for any fixed set of times.
gauge invariance
Green's function!gauge-invariant

current operator
We end this section with a discussion of the current operator and how it is related to gauge invariance.
The current operator is calculated in a straightforward way via the commutator of the Hamiltonian with the charge polarization operator. We do not provide the details here, but simply note that they can be found in most many-body physics textbooks. The end result is that

    j=∑_σ∑_k v_kc_kσ^† c_kσ^†.

Here, we have the band velocity is given by v_k=∇ϵ_k. The gauge-invariant form of the expectation value for the current is then

    ⟨ j⟩ =-i∑_σ∑_k v_k G̃_kσ^<(t,t).

We can re-express in terms of the original Green's functions, by shifting k→ k-A(t). This yields the equivalent formula

    ⟨ j⟩ = -i∑_σ∑_k v_k-A(t)G_kσ^<(t,t).

As you can easily see, these two results are identical by simply performing a uniform shift of the momentum in the summation. Generically, this is how different expectation values behave—in the gauge-invariant form, one evaluates them the same as in linear response, but uses the gauge-invariant Green's functions, while in the non-gauge-invariant form, one evaluates in a form similar to linear response, but appropriately shifts the momentum in the terms that multiply the Green's functions (which is a beyond-linear-response “correction”).
current operator



§ NONEQUILIBRIUM DYNAMICAL MEAN-FIELD THEORY


dynamical mean-field theory!nonequilibrium
The premise of dynamical mean-field theory (DMFT) is that the self-energy of an impurity in an appropriately chosen time-dependent field is the same as the self-energy on the lattice. In general, we know this cannot be true, because the lattice self-energy depends on momentum, and hence is not local. But, if we choose the hopping matrix elements to scale like t=t^*/2√(d) for a d-dimensional lattice, then one can rigorously show that in the limit d→∞, the self-energy does indeed become local <cit.>. So, in this limit, the DMFT assumption holds and it provides another limiting case where one can exactly solve the many-body problem (the other limit being that of d=1).
dynamical mean-field theory!nonequilibrium

dynamical mean-field theory!nonequilibrium
dynamical mean-field theory!iterative algorithm
The physical picture to keep in mind, for DMFT, is to think of looking at what happens just on a single site  of the lattice. As a function of time, we will see electrons hop onto and off of the site. If we can adjust the time-dependent field for the impurity in such a way that it produces the same behavior in time on the impurity, then the impurity site will look like the lattice. This motivates the following iterative algorithm <cit.> to solve these models: (1) Begin with a guess for the self-energy (Σ=0 often works); (2) sum the momentum-dependent Green's function over all momentum to determine the local Green's function

    G^c_iiσ(t,t')=∑_k { [ G_kσ^c(U=0) ]^-1-Σ}^-1_t,t^',

with G_kσ^c(U=0) the noninteracting momentum-dependent contour-ordered Green's function on the lattice and the subscript denoting the t,t^' element of the matrix inverse) ; (3) Use Dyson's equation to determine the effective medium G_0 via

    [G_0σ^c ]^-1_t,t^'= [G_iiσ^c ]^-1_t,t^'+Σ^c(t,t^');

(4) solve the impurity problem in the extracted effective medium (which determines the time evolution on the impurity); (5) extract the self-energy for the impurity using an appropriately modified Eq. (<ref>); and (8) repeat steps (2–7) until the self-energy has converged to the desired level of accuracy. The difference from the conventional iterative DMFT algorithm is that the objects worked with are now continuous matrix operators in time rather than simple functions of frequency. Sometimes this creates a roadblock for students to follow how the algorithm works, but there is no reason why it should. Finally, we note that because each inverse of a Green's function here implies a matrix inversion, we will want to organize the calculation in such a way as to minimize the inversions. This will be explained in more detail when we discuss numerical implementations below.
dynamical mean-field theory!nonequilibrium
dynamical mean-field theory!iterative algorithm

bandstructure!hypercubic lattice
One of the simplifications in equilibrium  DMFT is that the fact that the self-energy has no momentum dependence allows us to perform the summation over momentum via a single integral over the noninteracting density of states. Things are not quite as simple for the nonequilibrium case, as we now discuss. The bandstructure for a hypercubic lattice in d→∞ is given by

    ϵ_k=-lim_d→∞t^*/√(d)∑_i=1^dcos(k_i).

The central limit theorem tells us that ϵ_k is distributed according to a Gaussian distribution via

    ρ(ϵ)= 1/√(π)t^*e^-ϵ^2.

In nonequilibrium, we have a second band energy to work with 

    ϵ̅_k=-lim_d→∞t^*/√(d)∑_i=1^dsin(k_i),

because the Peierls substitution shifts cos(k_i)→cos(k_i)cos(A_i(t))+sin(k_i)sin(A_i(t)). The joint density of states becomes the product of two Gaussians, given by

    ρ(ϵ,ϵ̅)=1/π t^*2e^-ϵ^2/t^*2-ϵ̅^2/t^*2.

This second band energy can be thought of as the projection of the band velocity onto the direction of the electric field. Hence, the computation of the local Green's function in step (2) is complicated by requiring a two-dimensional integration rather than a one-dimensional integration. This does make the numerics require significant additional resources, especially because the integrands are matrix-valued. 
bandstructure!hypercubic lattice

dynamical mean-field theory!impurity solver
The most challenging part of the algorithm is step (4)—solving the impurity problem in nonequilibrium. Unfortunately, there are few techniques available to do this. Monte Carlo methods suffer from the so-called phase problem, which might be able to be tamed using expansions about a perturbative solution and restricting updates to be nearly local in time via the inchworm algorithm <cit.>. The other choice is perturbation theory either in the interaction (which often does not work so well, except for electron-phonon coupling) or the hybridization of the impurity model (which works well at half filling). Here, we make a different choice, and choose a simplified model that can be solved exactly, the so-called spinless Falicov-Kimball model <cit.>:

    ℋ=∑_ijt_ijc_i^† c_j^†-μ∑_i c_i^† c_i^†+E_f∑_i f_i^† f_i^†+U∑_i c_i^† c_i^† f_i^† f_i^†.

This model involves the interaction of conduction electrons (c) and localized electrons (f). The conduction electrons can hop on the lattice (we usually take the hopping only between nearest neighbors on a hypercubic lattice in the d→∞ limit) and have a chemical potential μ. The localized electrons have a site energy E_f. Both electrons interact with an on-site interaction U.
dynamical mean-field theory!impurity solver

Falicov-Kimball model
The Falicov-Kimball model describes a rich set of physics, but it should not be viewed as a paradigm for all strongly correlated electrons. This is because it has some aspects that are not seen in more common models like the Hubbard model. For example, it is never a Fermi liquid when U 0, the conduction electron density of states is independent of temperature in the normal state, and it is never a superconductor. But, it does display a lot of rich physics including a Mott metal-insulator transition, ordered charge-density-wave phases at low temperature and even phase separation when the densities of the two species are far enough apart. In addition, the f-electron density of states does exhibit strong temperature dependence, as expected in strongly correlated systems. It also displays orthogonality catastrophe behavior. The main interest in the nonequilibrium context has focused on the fact that it has a Mott transition. The most important aspect of the model is that it can be solved exactly within DMFT—both in equilibrium and in nonequilibrium.
Falicov-Kimball model

dynamical mean-field theory!impurity solver
We will not go through the details of the equation of motion for the Green's functions in the Falicov-Kimball model to show how one can solve it. The procedure is most efficiently performed using a path-integral formulation because the time-dependent field on the impurity cannot be easily expressed in terms of a Hamiltonian (unless one introduces a bath for the impurity which does this). Instead, one can employ functional methods to exactly compute the functional derivative of the partition function with respect to the dynamical mean-field, which then is employed to extract the Green's function. Then, because the f-electron number is conserved in the model, we obtain the Green's function of the impurity by summing over appropriately weighted combinations of the solution with no f-electrons and with one f-electron. The end result is that the impurity Green's function for the Falicov-Kimball model in an effective medium is the following:

    G_0^c(t,t^')=(1-w_1)G_0^c(t,t^')+w_1 [  ( G_0^c )^-1-U𝕀 ]^-1_t,t^'.

The symbol w_1=⟨ f^† f⟩ is the filling of the localized electrons (per lattice site). Since the f-electron number operator commutes with ℋ, the localized electron number is a constant of the motion and does not change with time, even in the presence of large electric fields. Note that we will work at half-filling where the electron densities for the conduction electrons and the localized electrons are each equal to 1/2. This point corresponds to μ=U/2 and E_f=-U/2. Details not presented here can be found in the original literature <cit.>.
dynamical mean-field theory!impurity solver



§ NUMERICAL STRATEGIES TO CALCULATE THE GREEN'S FUNCTION


Now we need to sort out just how we put this algorithm on a computer and implement it. This discussion closely follows Ref. <cit.>. As mentioned above, the challenge is that we are working with continuous matrix operators, with integration over the contour corresponding to matrix multiplication of these operators. Such objects cannot be directly put onto a computer. Instead, we have to first discretize the problem, replacing the continuous operators by 
discrete matrices. There are a few subtle points that are involved corresponding to the integration measure over the contour when we do this, but we  will precisely describe how we achieve this below. In order to recover the results for the continuous matrix operators, we need to extrapolate the discretized results to the limit where the discretization size goes to zero. This is done in the most mundane way possible—simply repeat the calculation for a number of different Δ t values and use Lagrange interpolation as an extrapolation formula to the Δ t→ 0 limit. We finally check sum rules of the data to verify that the extrapolation procedure worked properly. Details follow.

Kadanoff-Baym-Keldysh contour
We first need to be concrete about how the discretization is performed. This involves N_t points on the upper real branch (ranging from t_min to t_max-Δ t), N_t points on the lower real branch (ranging from t_max to t_min+Δ t), and 100 points along the imaginary axis (ranging from t_min to t_min-iβ+0.1i, with β=10); hence Δ t=(t_max-t_min)/N_t. We often find that fixing the number of points on the imaginary time branch rather than scaling them in the same fashion as on the real axis does not cause any significant errors, but it allows the calculations to be performed with fewer total number of points.
The discrete time values on the contour become

    t_j   =   -t_min+(j-1)Δ t,       1≤ j≤ N_t,
       =   t_max-(j-N_t-1)Δ t,    N_t+1≤ j≤ 2N_t,
       =   t_min-0.1i(j-2N_t-1),  2N_t+1≤ j≤ 2N_t+100,

where we used the fact that the discretization along the imaginary axis is always fixed at Δτ=0.1 in our calculations (and we pick the initial temperature to be T=0.1t^* or β=10).  We use a leftpoint rectangular integration rule for discretizing integrals over the contour, which is implemented as follows:

    ∫_c dt f(t)=∑_i=1^2N_t+100W_if(t_i),

where the weights satisfy

    W_j   =   Δ t,       1≤ j≤ N_t,
       =   -Δ t,      N_t+1≤ j≤ 2N_t,
       =   -0.1i,     2N_t+1≤ j≤ 2N_t+100.

The leftpoint integration rule evaluates the function at the “earliest” point (in the sense of time ordering along the contour, see Fig. <ref>) in the time interval that has been discretized for the quadrature rule (which is the left hand side of the interval when we are on the upper real branch and right hand side when we are on the lower real time branch).
Kadanoff-Baym-Keldysh contour

Green's function!boundary condition
One of the important aspects of many-body Green's functions is that they satisfy a boundary condition. This is what determines whether the Green's function is bosonic or fermionic. For example, you should already be familiar with the fact that the thermal Green's functions are antiperiodic on the imaginary time axis. The contour-ordered Green's function also satisfies a boundary condition where we identify the points t_min with t_min-iβ. One can show from the definition of the contour-ordered Green's function and the invariance of the trace with respect to the ordering of operators that G^c_iiσ(t_min,t^')=-G^c_iiσ(t_min-iβ,t^')
and G^c_iiσ(t,t_min)=-G^c_iiσ(t,t_min-iβ).  The proof is identical to how it is done for the thermal Green's functions. It involves cyclically moving a creation operator from the left to the right and commuting the e^-βℋ term through it which employs the fact that a cyclic permutation of elements in the product of a trace does not change the value of the trace.
Green's function!boundary condition

delta function!contour
The delta function changes sign along the negative real time branch, and
is imaginary along the last branch of the contour in order to satisfy the property that ∫_c dt^'δ_c(t,t^')f(t^')=f(t).  In addition, we find that the numerics work better if the definition of the delta function is done via “point splitting”  (when we calculate the inverse of a Green's function) so that the delta function does not lie on the diagonal, but rather on the first subdiagonal matrix (in the limit as Δ t→ 0 it becomes a diagonal operator). Because we identify the times t_min and t_min-iβ, the point splitting approach to the definition of the delta function allows us to incorporate the correct boundary condition into the definition of the discretized delta function. Hence, we define the discretized delta function in terms of the quadrature weights, in the following way

    δ_c(t_i,t_j)   =   1/W_iδ_ij+1,    for integration over j,
       =   1/W_i-1δ_ij+1,  for integration over i,

where t_i and t_j are two points on the discretized contour as described in Eq. (<ref>), and W_i are the quadrature weights described in Eq. (<ref>). We have a different formula for integration over the first variable versus integration over the second variable because we are using the leftpoint quadrature rule. Note that the formulas in Eqs. (<ref>) and (<ref>) hold only when i 1. When i=1, the only nonzero matrix element for the discretized delta function is the  1,j=2N_t+100 matrix element, and it has a sign change due to the boundary condition that the Green's function satisfies.
The discretization of the derivative of the delta function on the contour is more complex. It is needed to determine the inverse of the G_0^c for the impurity.  The derivative is calculated by a two-point discretization that involves the diagonal and the first subdiagonal.  Since all we need is the discrete representation of the operator [i∂^c_t+μ]δ_c(t,t^'), we summarize the discretization of that operator as
follows

    [i∂_t+μ]δ_c(t_j,t_k)=i1/W_jM_jk1/W_k,

with the matrix M_jk satisfying

    M_jk= (
    [         1         0         0                                               ...                                          1+iΔ t μ; -1-iΔ t μ         1         0                                               ...                                                 0;         0 -1-iΔ t μ         1         0;                                       ⋱;                                       0  -1+iΔ tμ         1         0;                                                 0  -1+iΔ tμ         1;                                                                               ⋱;                                                                                    -1-Δτμ         1;                                                                                                             ⋱;                                                                                                                  -1-Δτμ         1 ] );

here Δτ=0.1. The top third of the matrix corresponds to the upper real branch, the middle third to the lower real branch and the bottom third to the imaginary branch.
Note that the operator [i∂_t^c+μ]δ_c is the inverse operator of the Green's function of a spinless electron with a chemical potential μ.  Hence the determinant of this operator must equal the partition function of a spinless electron in a chemical potential μ, namely 1+exp[βμ]. This holds because detG^-1_non=𝒵_non for noninteracting electrons. Taking the determinant of the matrix M_jk (by evaluating the minors along the top row) gives

    M   =   1 +(-1)^2N_t+N_τ(1+iΔ tμ)(-1-iΔ tμ)^N_t-1
       ×   (-1+iΔ tμ)^N_t(-1-Δτμ)^N_τ,
       ≈   1+(1+Δτμ)^N_τ+O(Δ t^2),

which becomes 1+exp[βμ] in the limit where Δ t,Δτ→ 0 (N_τ is the number of discretization points on the imaginary axis).  This shows the importance of the upper right hand matrix element of the operator, which is required to produce the correct boundary condition for the Green's function. It also provides a check on our algebra. In fact, we chose to point-split the delta function when we defined it's discretized matrix operator precisely to ensure that this identity holds.
delta function!contour

We also have to show how to discretize the continuous matrix operator multiplication and how to find the discretized approximation to the continuous matrix operator inverse.  As we mentioned above, the continuous matrix operator is described by a number for each entry i and j in the discretized contour. We have to recall that matrix multiplication corresponds to an integral over the contour, so this operation is discretized with the integration weights W_k as follows

    ∫_c dt̅ A(t,t̅)B(t̅,t^')=∑_k A(t_i,t_k)W_kB(t_k,t_j).

Thus we must multiply the columns (or the rows) of the discrete matrix by the corresponding quadrature weight factors when we define the discretized matrix.  This can be done either to the matrix on the left (columns) or to the matrix on the right (rows). To calculate the inverse, we recall the definition of the inverse for the continuous matrix operator

    ∫_c dt̅ A(t,t̅)A^-1(t̅, t^')=δ_c(t,t^'),

which becomes the following

    ∑_k A(t_i,t_k)W_kA^-1(t_k,t_j)=1/W_iδ_ij,

in its discretized form. Note that
we do not need to point-split the delta function here.  Hence, the inverse of the matrix is found by inverting the matrix W_iA(t_i,t_j)W_j, or, in other words, we must multiply the rows and the columns by the quadrature weights before using conventional linear algebra inversion routines to find the discretized version of the continuous matrix operator inverse. This concludes the technical details for how to discretize and work with continuous matrix operators.

Green's function!local
The next technical aspect we discuss is how to handle the calculation of the local Green's function. Since the local Green's function is found from the following equation:

    G_ii^c(t,t^')=1/π t^*2∫ dϵ∫ dϵ̅ e^-ϵ^2/t^*2-ϵ̅^2/t^*2 [ (𝕀-G_ϵ,ϵ̅^c,nonΣ^c )^-1G_ϵ,ϵ̅^c,non ]_t,t^'

with the noninteracting Green's function discussed earlier and given by

    G_ϵ,ϵ̅^c,non(t,t^')=i[ f(ϵ-μ)-θ_c(t,t^') ]
    e^-i∫_t^'^t dt̅ [ cos ( A(t̅) )ϵ-sin (A(t̅) )ϵ̅ ].

Note that we must have A(t)=0 before the field is turned on. Of course an optical pump pulse also requires ∫_-∞^∞ A(t)dt=0, because a traveling electromagnetic wave has no dc field component. We did not compute the matrix inverse of the noninteracting Green's function because re-expressing the formula in the above fashion makes the computation more efficient (because matrix multiplication requires less operations than a matrix inverse does).
Green's function!local

This step is the most computationally demanding step because we are evaluating the double integral of a matrix-valued function and we need to compute one matrix inverse for each integrand. Fortunately, the computation for each (ϵ,ϵ̅) pair is independent of any other pair, so we can do this in parallel with a master-slave algorithm. In the master-slave algorithm, one CPU controls the computation, sending a different (ϵ,ϵ̅) pair to each slave to compute its contribution to the integral and then accumulating all results. The sending of the pairs to each slave is simple. One has to carefully manage the sending of the results back to the master, because the data is a large general complex matrix and they are all being received by one node. This is precisely the situation where a communication bottleneck can occur.

Gaussian integration
In order to use as few integration points as possible, we employ Gaussian integration, because the bare density of states for both the ϵ and the ϵ̅ integrals is a Gaussian. We found it more efficient to average the results with two consecutive numbers of points in the Gaussian integration (like 100 and 101) instead of just using 200 points (which would entail twice as much computation). 
Gaussian integration

The rest of the DMFT loop requires serial operations and is performed on the master to reduce communications. One does have to pay attention to convergence. In metallic phases, and with strong fields, the code will converge quickly But when the field is small or the interactions large, one might not be able to achieve complete convergence. Often the data is still good, nevertheless. One also should not make the convergence criterion too stringent. Usually 4 digits of accuracy is more than enough for these problems. In most cases, one will need to iterate as many as 50-100 times for hard cases. But many results can be obtained  with ten iterations or less. 

Finally, one has to repeat the calculations for different discretizations and extrapolate to the Δ t→ 0 limit. In order to use as much data as possible, it is best to first use a shape-preserving Akima spline  to put all data on the same time grid and then use Lagrange interpolation as an extrapolation method to the Δ t→ 0 limit. It is critical that one employs a shape-preserving spline, otherwise it will ruin your results. In addition, we  find quadratic extrapolation usually works best, but sometimes had to look at higher-order extrapolations. It is best to extrapolate the data after the desired quantity has been determined as a function of t for all times on the grid. For example, one would first determine the current by using the lesser Green's function and summing over all momentum for each time and then extrapolate the final current values as a function of time instead of extrapolating the Green's functions, since doing the latter requires a two-dimensional Lagrange extrapolation formula and uses large amounts of memory.



sum rules
Finally we mention that the accuracy can be tested by calculating spectral moment sum rules. These are known exactly for the retarded Green's function (through third-order).  At half-filling they simplify tremendously. We will explore how this works when we discuss some examples in the next section. It is important to perform this test on all of your data. It is easy to do and ensures the data is high quality.
sum rules

We end this section with a discussion of other methods. Tran showed how one can reduce the size of the matrices in half by working with the retarded and lesser Green's functions only <cit.>, but such an algorithm is harder to code and you need to solve for twice as many Green's functions, so it reduces the computational time by only a factor of two (it may allow for slightly larger matrices to be handled as well). There is an alternative way to discretize and solve this problem that integrates the Dyson equation in time in a causal fashion, producing the Green's function along the boundary of the matrix as time increases <cit.>. It is very efficient, but requires you to be able to determine the self-energy in a causal fashion from your data at all earlier times. This is always possible when perturbation theory is used to determine the self-energy. But, while there should be a way to do this for the Falicov-Kimball model, the explicit method has not yet been discovered.  For the Hubbard model, we have no exact methods available.



§ EXAMPLES


We now illustrate how all of this machinery works. The idea is not to provide an exhaustive review of the research work that has been completed, but rather to illustrate how these ideas are concretely put into action. As such, we will pick and choose the examples more for their didactic features than for their importance, or perhaps even their physical relevance.

self-energy

self-energy

self-energy
To start, we focus on equilibrium. The DMFT for the Falicov-Kimball model was solved exactly in equilibrium by Brandt and Mielsch <cit.> and is summarized in a review <cit.>. We take the exact result as a function of frequency and Fourier transform to time. Then we compare that result with the result that comes out of the nonequilibrium formalism, to understand its accuracy. In all of this work, we solve the Falicov-Kimball model on a hypercubic lattice at half filling. In Fig. <ref>, we plot the (local) lesser self-energy Σ^<(t). The calculation used a fixed Δ t=0.05 and has U=0.5. The different curves are for various different integration schemes over the noninteracting density of states. In one case, we average the N=54 and N=55 Gaussian formulas. In another, we do the same, but with N=100 and N=101. We also compare to a much denser trapezoidal formula with 1000 points. Additionally, we plot the exact result. One can see the different integration schemes are quite close to each other for times near zero, but they begin to deviate at large times. The inset focuses on large negative time and one can quickly conclude that there is a maximum absolute value of time for which the results are accurate for any integration scheme. Beyond that, the system starts to generate increasing amplitude oscillations. The disagreement at short times from the exact result stems from the fact that these calculations are at a fixed Δ t—no scaling to the Δ t→ 0 limit were taken. This gives a sense of the accuracies we can hope to attain with this approach.
self-energy

Green's function

Green's function

self-energy
In Fig. <ref>, we plot the Fourier transform of the time data as a function of frequency. The number of Gaussian points used is N=54, 55. The time steps are varies and one can see that they are approaching the exact result. If we use Lagrange extrapolation here, we would get quite close to the exact result, but we do not include that plot here because it would be too many lines close to each other. The oscillations in the tail of the data in Fig. <ref> is responsible for the oscillations at negative frequency seen in the data. Those oscillations will remain even after scaling to the Δ→ 0 limit. One can clearly see how the extrapolation method works for this case.
self-energy

Green's function
In Fig. <ref>, we do a similar plot in the time domain for the lesser Green's function. Note that the real part of the lesser Green's function (inset) is an odd function when we are at half-filling due to particle-hole symmetry. Similarly, the imaginary part in the main panel is even for the exact result. One can clearly see how the extrapolation will work to approach the exact result if we did the full extrapolation of this data. 
Green's function

density of states!equilibrium

density of states!equilibrium

Mott transition
Now that we have gotten our feet wet with the numerics, we are ready to discuss some physics. At half-filling, there are enough conduction electrons to fill half of the sites of the lattice and similarly for the localized electrons. In this case, if the repulsion between the two species is large enough, they will avoid each other and the net effect is that the system becomes an insulator because charge motion is frozen due to the high energy cost for double occupancy. We can see this transition occur in Fig. <ref>, which is the local density of states at half-filling for different U values. As U is increased we evolve from the initial result, which is Gaussian for no interactions to results where a “hole” is dug into the density of states until if opens a gap at a critical value of U called the Mott transition. Beyond this point, the system has a gap to charge excitations.
Mott transition

density of states!nonequilibrium

density of states!nonequilibrium

density of states!nonequilibrium
Next, we illustrate how the density of states evolves in a transient fashion in time after a dc field is turned on. In Fig. <ref>, we show a series of frames at different time snapshots that plot the transient density of states (imaginary part of the local retarded Green's function) as a function of frequency. One can see that the system starts off in a near Gaussian density of states and then develops features that are quite complex. The infinite-time limit is solved using Floquet nonequilibrium DMFT theory <cit.>; unfortunately, we are not able to describe the details for how that problem is solved. But there are a few important features to notice. Both before the field is applied and long after it has been applied, the system has a positive semidefinite density of states. This is the same as in equilibrium and it allows the density of states to be interpreted probabilistically. But in the transient regime, it often becomes negative. Furthermore, because the density of states is determined via a Fourier transformation, one cannot just say it corresponds to a spectra at a specific average time. Instead it senses times near the average time governed by how rapidly the Green's function decays in the time domain.
density of states!nonequilibrium

density of states!nonequilibrium
Wannier-Stark ladder
There is a lot of physics in these figures. If there was no interactions, the system would undergo Bloch oscillations, because the electrons are not scattering. The density of states then is the Fourier transform of a periodic function, which leads to delta function peaks forming the so-called Wannier-Stark ladder, with the separations given by the dc field E that is applied to the system. When interactions are added in, the Wannier-Stark ladder is broadened, but also split by the interaction U. This occurs because a delta function has a zero bandwidth and hence is highly susceptible to the Mott transition, even for relatively small interactions.
density of states!nonequilibrium
Wannier-Stark ladder

current
Bloch oscillations

current
Bloch oscillations

current
Bloch oscillations
Next, we move on to examining the transient current. We use the same  case we have been examining throughout this brief summary—U=0.5 and β=0.1. Here the dc field is E=1. One can see the current starts of as a weakly damped Bloch oscillations (underdamped). It dies off and remains quiescent for some time and then starts to recur at the longest times. The characteristic Bloch oscillations occurs for both metals and insulators, but it is damped much more rapidly in insulators because they interact more strongly. This is one of the common observables measured in a nonequilibrium experiment. But it is not measured directly, because oscilloscopes are not fast enough to see them.
current
Bloch oscillations

sum rules

sum rules

sum rules
Finally, we show how the sum rules hold and illustrate why it is important to scale results to the Δ t→ 0 limit. We pick a case which is challenging to calculate E=1 and U=2. This is in the Mott insulating phase, where the numerics are much more difficult to manage. We primarily use the sum rules to indicate whether the calculations are accurate enough that they can be trusted. As shown in Fig. <ref>, we can see that the raw data can be quite bad, but the final scaled result ends up being accurate to 5% or less! Note how the results are worst on the equilibrium side (to the left) than on the nonequilibrium side (to the right) and they approach the exact results on the nonequilibrium side. Hence, these calculations are most accurate in moderate to large fields.
sum rules

This ends our short representative tour of some numerical results that can be calculated with this nonequilibrium DMFT approach.



§ CONCLUSIONS


These lecture notes have been somewhat brief due to the page constraints of this contribution (and the time constraints of its author). Nevertheless, I have tried to present enough information here that you can follow the logic, reasoning, and even develop the formalism for yourself if you want to engage in these types of calculations. The field of nonequilibrium many-body physics is wide open. There are new and interesting experiments performed every day and the theory to describe them still needs further development. We do have a number of successes. The approach can be used for optical conductivity <cit.>, time-resolved angle-resolved photoemission <cit.>, Raman scattering <cit.>, x-ray photoemission spectroscopy and x-ray absorption spectroscopy <cit.>, and resonant inelastic x-ray scattering <cit.>. The challenge is always about how far out in time can a simulation go. 
If you have  good idea for a nonequilibrium solver for a Hubbard impurity, I encourage you to give it a try. We really need it. I also hope you will enjoy working with many-body Green's functions in the future. They are truly wonderful!



§ ACKNOWLEDGMENTS

This work was supported by the Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering under contract number DE-FG02-08ER46542. It was also supported by the McDevitt bequest at Georgetown University. I would also like to thank the organizers for the invitation to contribute to this school.


99
lehmann 
H. Lehmann, Nuovo Cim. 11, 342 (1954).
agd
A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Prentice Hall, New York, 1963).
luttinger_ward
J. M. Luttinger and J. C.  Ward, Phys. Rev. 118, 1417 (1960).
matsubara
T. Matsubara, Prog. Theor. Phys. 14, 351 (1955). 
baym_kadanoff
L. P. Kadanoff and G. Baym, Quantum statistical mechanics (Benjamin, New York, 1962).
keldysh
L. .V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964) in Russian; [Sov. Phys. JETP 20, 1018 (1965)].
dmft_review
A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
freericks1
 J. K. Freericks, V. M. Turkowski, and V. Zlatić, Phys. Rev. Lett. 97, 266408 (2006).
freericks2
J. K. Freericks,  Phys. Rev. B 77, 075109 (2008).
noneq_rmp_rev
H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner,
Rev. Mod. Phys. 86, 779 (2014).
brandt_mielsch1
U. Brandt and C. Mielsch, Z. Phys. B: Condens. Matter 75, 365 (1989).
brandt_mielsch2
U. Brandt and C. Mielsch, Z. Phys. B: Condens. Matter 79, 295 (1990).
brandt_neq
U. Brandt and M. P. Urbanek, Z. Phys. B: Condens. Matter 89, 297 (1992).
wagner
 M. Wagner, Phys. Rev. B 44, 6104 (1991).
peierls
 R. E. Peierls, Z. Phys. 80, 763 (1933).
joura
A. Joura, Static and dynamic properties of strongly correlated lattice models under electric fields (Dynamical mean field theory approach), Ph. D. thesis, Georgetown University (2014).
falicov_kimball
L. M. Falicov and J. C. Kimball, Phys. Rev. Lett. 22, 997 (1969).
bertoncini_jauho
R. Bertoncini and A.-P. Jauho, Phys. Rev. B 44, 3655 (1991).
metzner
W. Metzner, Phys. Rev. B 43, 8549 (1991).
jarrell
M. Jarrell, Phys. Rev. Lett. 69, 168 (1992).
inchworm
G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis,
Phys. Rev. Lett. 115, 266802 (2015).
freericks_ugc
 J. K. Freericks, V. M. Turkowski, and V. Zlatić, “Real-time formalism for studying the nonlinear response of `smart' materials to an electric field,” in Proceedings of the HPCMP Users Group Conference 2005, Nashville, TN, June 28–30, 2005 edited by D. E. Post (IEEE Computer Society, Los Alamitos, CA, 2005), pp. 25–34.
tran
M.-T. Tran, Phys. Rev. B 78, 125103 (2008).
contour_alt
G. Stefanucci and R. van Leeuwen, Nonequilibrium many-body physics of quantum systems: A modern introduction (Cambridge University Press, Cambridge, 2013).
freericks_review
J. K. Freericks and V. Zlatić, Rev. Mod. Phys. 75, 1333 (2003).
freericks_turk_rev
V. Turkowski and J. K. Freericks, “Nonequilibrium dynamical mean-field theory of strongly correlated electrons,” in Strongly Correlated Systems: Coherence and Entanglement, edited by J. M. P. Carmelo, J. M. B. Lopes dos Santos, V. Rocha Vieira, and P. D. Sacramento (World Scientific, Singapore, 2007), pp. 187–210.
freericks_neq_hubbard
A. V. Joura, J. K. Freericks, and T. Pruschke, Phys. Rev. Lett. 101, 196401 (2008). 
freericks_joura
J. K. Freericks and A. V. Joura, “Nonequilibrium density of states and distribution functions for strongly correlated materials across the Mott transition,” in Electron transport in nanosystems, edited by J. Bonca and S. Kruchinin (Springer, Berlin, 2008) pp. 219–236.
aoki_floquet
N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B 78, 235124 (2008).
ugc2
J. K. Freericks, V. M. Turkowski, and V. Zlatić, “Nonlinear response of strongly correlated materials to large electric fields,” in Proceedings of the HPCMP Users Group Conference 2006, Denver, CO, June 26–29, 2006, edited by D. E. Post (IEEE Computer Society, Los Alamitos, CA, 2006), pp. 218–226.
lex_oc
 A. Kumar and A. F. Kemper, preprint arXiv:1902.09549.
freericks_trarpes
J. K. Freericks, H. R. Krishnamurthy and T. Pruschke, Phys. Rev. Lett. 102, 136401 (2009).
freericks_raman
O  P. Matveev, A. M. Shvaika, T. P. Devereaux, and J  K. Freericks,
Phys. Rev. Lett. 122, 247402 (2019).
freericks_xas
A. M. Shvaika and J. K. Freericks, Cond. Mat. Phys. 15, 43701 (2012); unpublished.
freericks_rixs
Y. Chen, Y. Wang, C. Jia, B. Moritz, A. M. Shvaika, J. K. Freericks, and T. P. Devereaux,
Phys. Rev. B 99, 104306 (2019).





























