Here we consider a system on a hypercubic lattice in the limit of infinite spatial dimensions (d → ∞) with a constant electric field E applied at a specific time and oriented along a diagonal of the lattice. The dynamical mean-field theory employs a local self-energy such that

In the case of a system initially in equilibrium at an initial temperature T = 1/β, the symbol 〈A(t)B(t′)〉 in Equations (2–7) is synonymous with Tre-βHeqA(t)B(t′)/Zeq, where Heq=H(tmin) is the initial Hamiltonian and Zeq=Tre-βHeq is the corresponding partition function.

The contour-ordered Green's function obeys the Dyson equation given by

where Gkσ0(t,t′) is the non-interacting (U = 0) Green's function [also defined by Equations (2–7), but with a Hamiltonian that has U = 0]. The integrals each range over the entire Kadanoff–Baym–Keldysh contour. Prior to the electric field being turned on, the non-interacting Hamiltonian is

Here, ϵ_k is the band structure for the lattice electrons. When the electrons move on a hypercubic lattice in infinite dimensions, the band energy is given by:

where a is the lattice spacing (and is set equal to 1); the nearest-neighbor hopping is rescaled via J=t*/2d and t^* is the unit of energy. In the paramagnetic phase, the spin index may be dropped.

The system is driven out of equilibrium by an applied electric field E(r, t) which can be expressed in terms of a scalar potential ϕ(r, t) and a vector potential A(r, t) via

We choose the Hamiltonian gauge (Φ = 0), so that the electric field is completely determined by the vector potential. The electric field then enters into the Hamiltonian via the Peierls' substitution which modifies the hopping amplitude with a time-dependent phase [33]:

We will be working with a spatially uniform, but time-varying electric field, which is given by a spatially uniform vector potential A(t). In this case, the problem remains translationally invariant, but note that a spatially uniform but time-varying electric field does mildly violate Maxwell's equations. We ignore those magnetic-field effects since they are small.

The non-interacting Green's function obeys the equation:

where δc(t,t′) is a generalization of the Dirac delta function onto the contour.

When the field is constant and is turned on at time t = 0, the vector potential becomes A(t) = −cEtθ(t) and the Peierls' substitution results in the transformation k→k-eA(t)ℏc=k+eEtθ(t)ℏ in the Hamiltonian. The non-interacting Green's function then depends on the band energy ϵ_k and on the band velocity ϵk¯=-limd→∞t*d∑i=1dsin(kia) in the direction of the electric field. Setting ℏ = 1 and c = 1, this leads to the following expression [34]:

which becomes

The dressed lattice Green's function is constructed from the non-interacting Green's function and the self-energy following Equation (10). The local Green's function Gσ(t,t′) is found by summing the momentum-dependent dressed Green's function Gkσ(t,t′) over all momentum, or, equivalently integrating over the joint density of states, ρ2(ϵ,ϵ¯). In the non-equilibrium DMFT formalism, the hybridization of the impurity to the dynamical mean field is given by

where we suppressed the spin indices.

To complete the self-consistency loop, we calculate the local Green's function from the dynamical mean field and use Dyson's equation (on the contour) to extract the self-energy. This is usually the bottleneck in the DMFT formalism. But in the case of the Falicov–Kimball model, this step can be determined analytically. We start from the spinless Falicov–Kimball model in the absence of a field

where w_i = 0, 1 is the occupation number operator of a localized fermions at site i. The corresponding impurity problem is solved by

where 〈w〉=∑i〈wi〉/N is the initial (equilibrium) filling of the localized fermions. Note that the Green's function is represented by sum of the matrix inverses of two continuous matrix operators.

The main shortcoming of the above non-equilibrium DMFT solution is that even in cases such as that of the Falicov-Model, where we can write an exact expression for the impurity Green's function, the non-equilibrium DMFT implementation remains computationally expensive both in terms of memory and time (the compute-intensive step is determining the local Green's function from the self-energy). It is for this reason that most studies are limited to examining the short-time transient only. One may alternatively be interested in studying the system after it has undergone its relaxation from an initial state and settled into its steady state. In this case, a steady-state non-equilibrium DMFT formalism can be formulated directly to study the steady-state regime. In this situation, we can use a non-interacting initial state with its non-interacting contour-ordered Green's function given by:

Switching to Wigner coordinates, as illustrated in Figure 3, with the relative time trel=t-t′ and average time tave=(t+t′)/2, we obtain:

in the long-time limit (t ≫ t_on and t′≫ton, with t_on the time the field is initially turned on). We denote the Green's function in Wigner coordinates by gk0(trel,tave). From Equation (22), it is easy to show that

or, equivalently

when the system is driven by a DC field.

Using the steady state assumption (in the long-time limit) that Σ(t, t′) = Σ(t − t′), one obtains that the dressed Green's functions also satisfy the expressions in Equations (24) and (25). For the choice of τ = (t + t′)/2, we obtain Gk-Eτ(t,t′)=Gk((t-t′)/2,(t′-t)/2). The local Green's function is G(t,t′)=∑kGk((t-t′)/2,(t′-t)/2), which, in terms of the relative time, becomes

Note how g_k depends on k only through ϵ_k and ϵ¯k leading to

with the joint density of states [13] given by

This means that in frequency space, the local Green's function is obtained from

One can similarly show that in the long-time limit, we have Gk(t+2πE,t′+2πE)=Gk(t,t′) or equivalently, gk(trel,tave-2πE)=gk(trel,tave) i.e., g_k is periodic with respect to t_ave with period 2π/E; this is the Bloch-Zener periodicity and ν_n = nE (with n an integer) are the Bloch frequencies. There is a subtle issue going on here, because we are working in this long-time limit. The formulas stated here hold when both times are much larger than t_on. Nevertheless, the Green's function has a periodic average time dependence in this long-time limit, and this is what we are describing here. As a result, Fourier transforming g_k(t_rel, t_ave) in this long-time limit, gives

where ω is a continuous variable while ν_l is quantized and the integral over t_ave can be restricted to be taken over just one Bloch period. In frequency space, the Dyson equation can then be written in the form

The periodicity in this long-time limit, maps the system to behave exactly as a many-body Floquet system acts. We now describe how to solve the steady-state problem within such a Floquet formalism, using discrete matrices within the Dyson equation for the Green's function defined within a frequency range determined by the Floquet period. We start by choosing a finite frequency range over which the Green's function has non-zero spectral weight. Next, we choose an integer L ≫ 1 such that Σ(ω) ≠ 0 only for ω ∈ [−ν_L, ν_L] and let ϕ be a real variable allowed to vary continuously between −E/2 and E/2. Equation (31) for (ω, ν_l) = (ϕ + ν_{−L + p}, ν_p) can be rewritten as follows:

and one can immediately see that the frequency ϕ is coupled only to frequencies shifted by integer multiples of the Bloch frequency. Here we have defined

with i, j, p being integer indices. Equation (32) provides a convenient way to solve the Dyson equation. For this purpose, it is necessary to obtain the Fourier transform of the non-interacting Green's function in Equation (22) to frequency space. The result is

with

and

where J_a(x) is the Bessel function of the first kind and P denotes the principal value. The different variables are defined by u=ω+μE, r=ϵ2+ϵ¯2, r′ = r/E, cosα=ϵϵ2+ϵ¯2 and sinα=ϵ¯ϵ2+ϵ¯2.

Note that the integer indices i, j, and p in Equation (32), are chosen such that i, j, p = 0, 1, 2, ⋯L, for the quantized frequencies ν_i, ν_j, and ν_p, which combine with the real variable ϕ ∈ [−E/2, E/2] to produce a dependence on a single continuous frequency that varies within the range of non-zero spectral weight [−ν_L, ν_L]. This single variable is relabeled as ω.

To complete the self-consistency loop for this non-equilibrium steady-state DMFT algorithm, the only missing ingredient is the impurity solver. In the case of the Falicov–Kimball model [which is initially described by the Hamiltonian in Equation (19)], the dressed Green's function is related to the non-interacting Green's function via

Combining this with the Dyson equation yields

The self-consistency loop is now complete and proceeds as follows:

This formalism is equivalent to using Floquet theory to solve the steady-state problem for a field driven system [20, 21, 35, 36]. It allows the characterization of the steady state for a given electric field and interaction strengths. In particular, one can obtain the local density of states (DOS) from ρ(ω) = −Im[g(ω)]/π.

In order to be able to describe the effects of the electric field, it is instructive to start by looking at the equilibrium local density of states (DOS) in Figure 4 obtained by the DMFT algorithm in equilibrium [16]. As the interaction strength is increased, we see a transition from a broadened Gaussian-like DOS in the weak-interaction regime to a gapped spectrum with two peaks in the strong-interaction regime. The spectral gap at ω = 0 is characteristic of the Mott-insulating regime, with the critical U occurring at U=2. Note that the DOS of the Falicov–Kimball model, in the normal state, does not depend on temperature.

The non-equilibrium steady state DOS shows much richer behavior than in equilibrium. Figures 5, 6 show the steady-state DOS for electric fields E = 0.5 and E = 1.0, respectively, for the same interaction strengths as the equilibrium DOS in Figure 4. In the weak-interaction regime, the DOS displays the Wannier–Stark ladder [37–39] with peaks centered around ω = nE with n = 0, ±1, ±2, ... and the amplitude is sharply suppressed away from the central peak. The peaks are broadened by a width governed by the interaction strength and the central peak has a dip separated by two side peaks at ±U/2. As the interaction is further increased, the initial peaks eventually merge, while preserving the central dip that arises from strong correlations. With stronger interactions this dip does not however become the full width of the gap due to the presence of “in gap” states created by the electric field.

In both the equilibrium and the steady-state regime, the system obeys the fluctuation-dissipation theorem and the lesser Green's function satisfy the relation

Where f_T(ω) is the Fermi distribution function at temperature T. It is interesting to now ask how do the dynamics transiently take us from equilibrium to the steady state?

In earlier work, we examined the transient dynamics of correlated quantum systems when they are driven away from equilibrium by a DC electric field [9]. Key findings included the fact that the DOS is only constrained by causality and takes its steady-state value as soon as the field is turned on (seeing this in the frequency domain takes a little longer, because one needs the retarded Green's function behavior to hold for a long-enough relative time span that a Fourier transformation can be performed). Since the system is isolated and the electric field remains constant after it is switched on, Joule heating should lead to an infinite-temperature state, if the system thermalizes. Without assuming the system to be thermal, an effective temperature at a time t can be obtained by comparing the energy of the system to that of the corresponding equilibrium system (in other words, using a thermometer based on the energy as a scale, without determining if the system has actually thermalized). The evolution of the effective temperature as a function of time is presented in Figure 7. We found that for some parameter regimes, this infinite-heating scenario did take place. Sometimes the infinite-temperature limit was approached monotonically, sometimes it was approached in an oscillatory fashion, with the effective temperature alternating between positive and negative temperatures en route to T = ∞.

We focus on the case of monotonic relaxation toward the infinite-temperature thermal state as it exhibits interesting manifestations of the fluctuation-dissipation theorem. Besides the effective temperature from an energy thermometer, several other quantities are used in characterizing the relaxation. These include the real part of the frequency-dependent lesser Green's function (Figure 8, left panel), and the current, the kinetic, potential and total energy (Figure 8, right panel). If the system is in a thermal state, then the lesser Green's function and the retarded Green's function satisfy the fluctuation-dissipation theorem. Namely, in frequency space, they obey Equation (42).

In this case, the lesser Green's function has a vanishing real part both in time and in frequency space. Figure 8, left panel, shows that the real part of the lesser Green's function relaxes monotonically toward zero (its infinite-temperature value). Additionally, it is easy to evaluate the total energy at infinite temperature. This is given by the dashed magenta line in Figure 8, right panel. To evaluate the transient total energy of the system, one starts from the known equilibrium energy at the initial temperature T = 0.1. Next, we integrate the inner product 〈j(t)〉·E up to time t to obtain the energy added to the system by Joule heating [here j (t) is the current at time t]. Adding this Joule heating energy to the initial energy produces the total energy of the system at time t. In Figure 8, one can see that in this case the total energy in the transient calculation approaches this limit monotonically, as the current also decays monotonically toward zero.

This decay of the current toward zero is to be expected at infinite temperature where particles are equally likely to move in opposite directions. The potential and kinetic energy also monotonically approach their infinite-temperature values. In this monotonic thermalization scenario, we find that even before reaching the infinite-temperature thermal state, the system appears to follow an evolution through successive quasi-thermal states approximately satisfying the fluctuation-dissipation theorem. This is illustrated in Figure 9 for the lesser Green's function as a function of frequency and average time. After the field is turned on, g^< switches from its equilibrium value (green line) to different transient values. Their evolution is pictured in grayscale with lighter shades indicating later average times. At infinite average time, g^< is expected to adopt the infinite-temperature fluctuation-dissipation result (blue dashed line). It is seen in the graph that at the last available simulation times, the transient already follows the fluctuation-dissipation theorem: the red line is the fluctuation-dissipation theorem at the effective temperature (using energy as the thermometer scale) corresponding to this average time and shows good agreement with the transient.

Another window into the relaxation of the field driven system is through the distribution function [40]. For the non-interacting system in equilibrium, this is simply the Fermi-Dirac distribution function. This equilibrium distribution function is modified by many-body effects for finite U values although the lineshape is essentially preserved as shown in Figure 10. We envision this situation as corresponding to fermionic atoms of two different masses trapped in a mixture on an optical lattice [24].

For the interacting system away from equilibrium, the distribution function is given by the gauge-invariant Wigner distribution nk(t)=-iGk+A(t)<(t,t) [41]. Here, n_k(t) represents the occupation of states in momentum space. The dependence on k is recast into a dependence on the band energy and the band velocity. Prior to the electric field being switched on, the system is in equilibrium at a temperature T = 0.1 and the Wigner distribution in momentum space is shown in Figure 11 for U = 0.25 and is similar for other interaction strengths.

Once the field is turned on, it is expected that the onset of the electric current will lead to Joule heating that will increase the temperature until the system has reached an infinite-temperature state where all points in momentum space, or equivalently in [band energy:E(k) ≡ ϵ_k; band velocity:V(k)≡ϵk¯]-space, are equally likely to be occupied. In this state, the current vanishes and the Wigner distribution is now homogeneous for all k. This is translated into the false color plot being entirely white for the plotted region. It was previously reported that the evolution between the initial configuration of Figure 11 to this uniform configuration goes through the development of different patterns that depend on the field and interaction strength. For instance, for a strong electric field and weak interaction (E = 2.0, U = 0.25), we observe in the false-color plots of the Wigner distribution as a function of time, the formation of ring-shaped disturbances on a timescale of 2π/U (Figures 12a,c,d). These rings spiral around the center of the graph on a timescale of 2π/E (Figure 12b) and we see the formation of additional rings after every 2π/U time step (Figures 12c,d). 2π/E is the period of Bloch oscillations while 2π/U is the timescale for the collapse and revival of Bloch oscillations [1, 42–45]. As the interaction is increased (E = 2.0, U = 1.0), these two timescales merge and we no longer observe well-separated individual rings but rather the growth of a single spiral that gradually scrambles the occupation of the states. We also observe at larger interaction strengths (E = 2.0, U = 3.0) that long-lived features with both a stable region of high occupation and a region of low occupation rotate around the origin at the Bloch period.