Characterizing the Non-equilibrium Dynamics of Field-Driven Correlated Quantum Systems

Recent experimental advances in ultrafast phenomena have triggered renewed interest in the dynamics of correlated quantum systems away from equilibrium. We review non-equilibrium dynamical mean-field theory studies of both the transient and the steady states of a DC field-driven correlated quantum system. In particular, we focus on the non-equilibrium behavior and how it relates to the fluctuation-dissipation theorem. The fluctuation-dissipation theorem emerges as an indicator for how the system thermalizes and for how it reaches a steady state. When the system thermalizes in an infinite temperature steady state it can pass through a succession of quasi-thermal states that approximately obey the fluctuation-dissipation theorem. We also discuss the Wigner distribution and what its evolution tells us about the non-equilibrium many-body problem.


I. INTRODUCTION
Strongly correlated systems include some of the most technologically promising materials of our time. The same quantum-mechanical complexity behind their most intriguing properties also renders them challenging to study. Much of the interest is motivated by recent experimental successes. For example, trapping and manipulating ultracold atomic gases in optical lattices provides a new platform for the controlled examination of strongly correlated systems in and out equilibrium 1,2 . In electronics, device miniaturization leads to the creation of nanoscale devices in which electrons experience strong electric fields and thus cannot be well approximated by linear-response theories 3,4 . Additionally, pump-probe spectroscopy offers a new avenue for the exploration of available electronic states in correlated materials 5 . These advances have revived interest in the fundamental behavior of quantum systems away from equilibrium.
There are many unanswered questions. How does an equilibrium quantum system that is suddenly driven out of equilibrium, subsequently relax to a thermal state [6][7][8][9] ? What governs the nonequilibrium driving of a system into a metastable state that is not found among known equilibrium phases 10,11 ? Answering these questions requires an accurate theoretical approach. Nonequilibrium dynamical mean-field theory 12 is a powerful tool to address these pressing questions.
In equilibrium, dynamical mean-field theory (DMFT) [13][14][15][16] treats spatial correlations in a mean-field fashion, while treating temporal correlations exactly. It is one of the most commonly used and successful methods for studying strongly correlated systems. It has also been extended to nonequilibrium 12,17,18 .
Field-driven correlated systems display nontrivial dynamics. When a DC electric field is applied, the system responds by creating an electric current that flows. That current also leads to Joule heating, which, if left unchecked, will heat the system to infinite temperature (where the current vanishes and the heating ultimately stops). The heating can also stop prematurely, with the system reaching a metastable (or possibly even thermal) state. Indeed, several scenarios for this relaxation process can occur depending on the specific details of the system 9 . In many of these situations, the fluctuation-dissipation theorem is no longer rigorously valid. However, it still remains an important concept in nonequilibrium, because it can be used to tell us when a system settles into a thermal steady state [19][20][21][22] . We review nonequilibrium DMFT manifestations of these properties in a mixture of heavy and light fermionic particles described by the Falicov-Kimball model 23,24 that starts in equilibrium and then has a DC electric field suddenly turned on.
The review is organized as follows. In section II, we present a brief description of the nonequilibrium Green's functions and introduce the Falicov-Kimball model as it relates to the Hubbard model. In section III, we describe the nonequilibrium DMFT solution first for the transient and then for the steady state. In section IV, we present a set of results that include manifestations of the fluctuation-dissipation theorem and provide insights into the rich dynamics of the field-driven strongly correlated quantum system.

II. NONEQUILIBRIUM GREEN'S FUNCTIONS
The many-body formalism for the nonequilibrium problem is similar to that of the equilibrium problem. All essential requirements for perturbation theory still hold away from equilibrium, except for the identification of the state at long times being identical to the state for the earliest times. Instead, it becomes necessary to evolve the system according to the Heisenberg representation for operators: start from the distant past, evolve forward to the time of physical interest and then evolve backward to the distant past (because of the Hermitian conjugate of the evolution operator in the Heisenberg representation) 25,26 . This gives rise to the Kadanoff-Baym-Keldysh contour illustrated in figure 1; the contour-ordered Green's functions has time ordering performed with respect to the advance along the contour. It is defined via where θ c (t, t ′ ) is the contour-ordered Heaviside function, which orders time with respect to the contour. It is equal to 1 if t is ahead of t ′ on the contour and is equal to 0 if t is behind t ′ on the contour. Hence, the nonequilibrium formalism automatically includes several different Green's functions depending on the location of each time argument on the contour. The lesser (t on lower, t ′ on upper), greater (t on upper, t ′ on lower), time-ordered (t and t ′ on upper), anti-time-ordered (t and t ′ on lower) Green's functions are respectively defined by the following when both time arguments of the contour-ordered Green's function are real: From these Green's functions, we can construct the so-called retarded and advanced Green's functions via Here c † kσ (t) and c kσ (t) are respectively the Heisenberg representation of the creation and the destruction operators for an electron of momentum k and spin σ at time t; θ is the ordinary Heaviside function, i. e., θ(t − t ′ ) = 0 if t < t ′ and θ(t, t ′ ) = 1 otherwise; {A, B} is the anticommutator of operators A and B. The symbol O is the expectation value taken of the operator O with respect to the initial thermal state: where H(t min ) is the initial Hamiltonian before any field is turned on. Several additional Green's functions are included in the formalism when at least one time argument is on the imaginary spur of the contour. These include the imaginary-time (or Matsubara) Green's functions with both times on the imaginary time vertical branch and the mixed-time Green's functions where one of the times is on the horizontal branch while the other is on the vertical branch (formulas not explicitly shown here). The nonequilibrium many-body formalism works directly with the contour-ordered Green's functions 12 . All other Green's functions can then be extracted from it. Note that these Green's functions are not all independent and are related by various identities and symmetries. One choice for two independent ones is the lesser and the retarded Green's functions, G < and G R . These determine most physical quantities. In particular, the retarded Green's function is related to the quantum states while the lesser Green's function is directly related to how those states are occupied by fermions.
We are interested in the dynamics of a strongly correlated system that starts initially in thermal equilibrium and then is driven out of equilibrium by turning on an external DC electric field. We describe how the fluctuationdissipation theorem manifests itself and how it can be utilized in characterizing relaxation through transient states to the (long-time) steady-state.
As an example, we consider a generalized Hubbard model for strongly correlated systems with spin-dependent hopping: FIG. 1. Kadanoff-Baym-Keldysh contour. The arrows indicate the direction of time ordering. In this example, t occurs before t ′ on the contour, even though t > t ′ , when thought of as real numbers.
The dynamical mean-field theory maps the lattice problem (here, a square lattice with hoping integral J between nearest-neighboring sites and Coulomb interaction U for doubly occupied sites) onto that of an impurity embedded into a selfconsistently determined bath with a hybridization function Λ(t, t ′ ). In equilibrium, the hybridization function only depends on the relative time t − t ′ .
Here ij represents unique nearest-neighbor pairs for sites i and j; J ijσ is the nearest-neighbor hopping integral; µ is the chemical potential; n iσ = c † iσ c iσ is the number operator for electrons of spin σ at site i; and U is the on-site repulsion for doubly occupied sites.
When J ijσ = J ijσ = J, we obtain the conventional Hubbard model, which is believed to contain the essential elements for high-temperature superconductivity in the cuprates and for this reason has been the subject of intense research activity. Despite its deceptive simplicity, the model remains unsolved except in one dimension 27 and in the limit of infinite dimensions 13 . But, when we have J ij↓ = J and J ij↑ = 0, i. e., if the electrons with spin up are held static while those with spin down are allowed to hop between nearest-neighbor sites, we obtain the Falicov-Kimball model, which can also be viewed as describing a mixture of heavy and light fermions on a lattice. This model has the advantage of displaying a Mott transition while also being more amenable to numerical methods than the Hubbard model. We will study this model at half-filling when there are as many electrons with spin up as spin down (or, equivalently, as many light, ↓, as heavy, ↑, fermions) and, in total, as many electrons as there are lattice sites.

III. NONEQUILIBRIUM DMFT
Dynamical mean-field theory was introduced to address strongly correlated systems in equilibrium and has been used successfully to describe key aspects of strongly correlated systems such as the metal-to-Mott insulator transition [13][14][15][16] . The method has successfully been adapted to studies of nonequilibrium systems 12,17,18 . It maps the lattice model onto an impurity embedded in a self-consistently determined bath as illustrated in figure 2. The method then relies on an accurate solution of the impurity problem and has been implemented with a multitude of impurity solvers with varied levels of success. These solvers include diagrammatic perturbative approaches 28,29 , Quantum Monte Carlo 30 and exact diagonalization 31,32 . The Falicov-Kimball model considered in this case has the advantage of being able to be solved exactly, as shown below.

A. Transient Nonequilibrium DMFT
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 is the initial Hamiltonian and Z eq = Tr e −βHeq is the corresponding partition function.
The contour-ordered Green's function obeys the Dyson equation given by where G 0 kσ (t, t ′ ) is the noninteracting (U =0) Green's function [also defined by Eqs. (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 noninteracting 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 * /2 √ d 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 34 : 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 noninteracting 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) in the Hamiltonian. The noninteracting Greens function then depends on the band energy ǫ k and on the band velocityǭ k = − lim d→∞ sin(k i a) in the direction of the electric field. Settingh = 1 and c = 1, this leads to the following expression 33 : which becomes The dressed lattice Green's function is constructed from the noninteracting Green's function and the self-energy following Eq. (10). The local Green's function G σ (t, t ′ ) is found by summing the momentum-dependent dressed Green's function G kσ (t, t ′ ) over all momentum, or, equivalently integrating over the joint density of states, ρ 2 (ǫ,ǭ). In the nonequilibrium 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 w i /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.

B. Nonequilibrium steady State DMFT
The main shortcoming of the above nonequilibrium 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 nonequilibrium 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 nonequilibrium DMFT formalism can be formulated directly to study the steady-state regime. In this situation, we can use a noninteracting initial state with its noninteracting contour-ordered Green's function given by: Switching to Wigner coordinates, as illustrated in figure 3, with the relative time t rel = t − t ′ and average time t ave = (t + t ′ )/2, we obtain: The dashed gray line represents the time at which the system is driven away from equilibrium by switching on the electric field (above and to the right of the dashed lines is where the field is turned on). The average time axis is illustrated by the red line across the diagonal and each average time has an associated set of relative times identified by points running along the blue lines.
in the long-time limit (t ≫ t on and t ′ ≫ t on , with t on the time the field is initially turned on). We denote the Green's function in Wigner coordinates by g 0 k (t rel , t ave ). From Eq. (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 Eqs. (24) and (25). For the choice of τ = (t + t ′ )/2, we obtain G k−Eτ (t, t ′ ) = G k ((t − t ′ )/2, (t ′ − t)/2). The local Green's function is G(t, t ′ ) = k G k ((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 G k (t+ 2π E , t ′ + 2π E ) = G k (t, t ′ ) or equivalently, g k (t rel , t ave − 2π E ) = g k (t rel , t ave ) 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 nonzero 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) with 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 Eq. (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 nonzero spectral weight [−ν L , ν L ]. This single variable is relabelled as ω.
To complete the self-consistency loop for this nonequilibrium 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 Eq. (19)], the dressed Green's function is related to the noninteracting Green's function via Combining this with the Dyson equation yields The self-consistency loop is now complete and proceeds as follows: 1. For a given Σ(ω), solve the Dyson equation for g ǫ,ǭ (ω, t ave → ∞); 2. Next, solve Eq. (29) for the local Green's function; 3. Then, solve Eq. (41) for Σ(ω);

Repeat steps 1 − 3 until convergence is achieved.
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(ω)]/π.

A. Steady-state density of states
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 nonequilibrium steady state DOS shows much richer behavior than in equilibrium. Figures 5 and 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 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?

B. Monotonic thermalization
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 of 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 towards 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 frequencydependent 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 Eq. (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 towards zero (its infinitetemperature 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 towards zero.
This decay of the current towards 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 38 . For the noninteracting 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 n k (t) = −iG < k+A(t) (t, t) 39 . 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. . 2π/E is the period of Bloch oscillations while 2π/U is the timescale for the collapse and revival of Bloch oscillations 1,42,43 . 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.

V. DISCUSSION
In general, the fluctuation-dissipation theorem is not satisfied once a system is driven away from equilibrium. This review describes nonequilibrium DMFT studies for both the transient and the steady-state of the Falicov-Kimball model, describing a Fermi-Fermi mixture of heavy-light particles, when it is driven away from equilibrium by a  constant electric field. We showed the complex range of relaxation scenarios exhibited by this nonequilibrium system. In particular, the density of states is fundamentally altered by the electric fields with the formation of Wannier-Stark ladders and a dielectric breakdown that arises with the presence of mid-gap states that are absent in equilibrium. This density of states however switches to its steady-state value rapidly after the field is turned on. For an isolated system, the relaxation to a steady state satisfying the fluctuation-dissipation theorem can then occur in several identified scenarios. We have particularly examined the monotonic thermalization scenario where the system monotonically goes to an infinite-temperature thermal state (satisfying the fluctuation-dissipation theorem) and evolves through consecutive quasi-thermal states (approximately satisfying the fluctuation-dissipation theorem) en-route. Ongoing work will take advantage of this property to develop an extrapolation scheme that bridges the gap between the transient and the steady state at minimal computational cost. We have also described how the key timescales that appear in the current in the form of Bloch oscillations and their collapse and revival, or beats, are manifested in the Wigner distribution function and its evolution towards infinite temperature where all states are equally occupied [n k = 0.5]. These results illustrate the rich physics of field-driven correlated quantum systems and the role the