Nonlinear Peltier Effect and the Nonequilibrium Jonson-mahan Theorem

We generalize the many-body formalism for the Peltier effect to the nonlinear/nonequilibrium regime corresponding to large amplitude (spatially uniform but time-dependent) electric fields. We find a relationship between the expectation values for the charge current and for the part of the heat current that reduces to the Jonson-Mahan theorem in the linear-response regime. The nonlinear-response Peltier effect has an extra term in the heat current that is related to Joule heating (we are unable to fully analyze this term). The formalism holds in all dimensions and for arbitrary many-body systems that have local interactions. We illustrate it for the Falicov-Kimball, Hubbard, and periodic Anderson models.


Introduction to heat transport and the Jonson-Mahan theorem
In 1834, Peltier introduced the notion that an electrical current also carries heat current with it [1].The Peltier effect is now viewed as the (linear-response) heat current that is generated when a small electric field is applied to a material held at a constant temperature; the proportionality between the d.c.heat current and the d.c.charge current is called the Peltier coefficient.In recent years, there has been much interest in the thermal transport of strongly correlated materials, because correlations, like in the Kondo effect, could lead to commercially viable thermal transport effects at low temperatures, perhaps resulting in low-temperature solid-state coolers (the Peltier coefficient is large near the Kondo temperature, which could be useful to employ for cooling applications).
In a series of two papers, in 1980 and 1990, Jonson and Mahan derived an exact relationship between the Kubo response formulas for charge and heat transport, which is called the Jonson-Mahan theorem [2,3].The Jonson-Mahan theorem can be summarized simply by saying that if one writes the linear-response conductivity in the form with f (ω) = 1/[1 + exp(βω)] the Fermi-Dirac distribution function, β = 1/T the inverse temperature, and τ (ω) a generalized many-body relaxation time (that includes all effects of vertex corrections, if appropriate), then the (linear-response) Peltier effect coefficient (determined from the ratio of the heat to charge currents) satisfies note that our convention for the Peltier coefficient is negative for electrons and positive for holes.The Jonson-Mahan theorem follows from the fact that one can find exact operator relationships between the charge-and heat-current correlation functions, which involve one derivative with respect to imaginary time.When analytically continued to the real axis, the time derivative becomes a power of frequency.The Jonson-Mahan theorem is important for two reasons: (i) first, it shows that all the manybody effects for charge transport are essentially the same for heat transport, so even though the operators are quite different, the correlation functions (including all effects of vertex corrections) retain the relationship they have for noninteracting systems, and (ii) second, the theorem allows for a consistent approximation scheme for charge and thermal transport -if one has an approximate result for the charge transport, then the use of the Jonson-Mahan theorem will produce the most reasonable approximation for the thermal transport.
In this contribution, we are interested in nonlinear-response effects.Hence we apply a large electric field to our material and watch how it evolves over time.If we keep the temperature fixed throughout the material, so there are no thermal gradients present, then we can study the charge current and the heat current induced by the electric field.The ratio of these two currents can be thought of as a nonlinear Peltier coefficient.Surprisingly, one can generalize the Jonson-Mahan theorem to the nonlinear response regime when a large spatially uniform electric field (with arbitrary time dependence) is applied.The generalization allows for a simple calculation of one piece of the heat current, which we carefully derive here.We find that the nonequilibrium heat current also contains an additional term related to the Joule heat production, which is not easily amenable to calculation.
Studying the nonlinear response of a strongly correlated material to the presence of a large electric field requires the use of the nonequilibrium many-body formalism of Kadanoff and Baym [4] and of Keldysh [5].This is because one can exactly solve the problem of a noninteracting material in the presence of an arbitrarily large electric field (for example, see [6] for the noninteracting Green's functions in an arbitrary uniform electric field in the infinite-dimensional limit), and then the many-body interactions can be added.Since the formalism does not use correlation functions (as in the Kubo-Greenwood formula [7,8]), one directly evaluates expectation values for the currents as a function of time; indeed, the solution in a field usually has a current flowing through it.
This nonequilibrium formalism requires physical fields to be added to the Hamiltonian in order to determine the response to those fields.Hence, it is difficult to determine the response to a temperature gradient, because a temperature gradient cannot be described by such a driving field.Instead, one needs to consider the variation of temperature throughout the material to be slow enough that the system can be divided into small blocks, with each block having enough atoms in it so that the concept of a temperature applies to the block.Then one can try to study the change in the temperature from block to block through the system and the currents that arise from those changes in temperature.One can even make contact with Kubo-like formulas by adding a gravitational field to the Hamiltonian, which couples to the energy current, and then showing the equivalence of the response to a gravitational field with the response to a temperature gradient [9].Such a procedure is difficult to generalize to the nonequilibrium/nonlinear-response regime, so we focus here on the Peltier effect, which occurs due to the presence of an external electric field at a constant temperature.

Many-body Hamiltonians and their charge and heat current operators
The derivation of the nonequilibrium Jonson-Mahan theorem should hold for any many-body Hamiltonian that has local interactions.We explicitly demonstrate the results for three popular models of strong electron correlations: (i) the Falicov-Kimball model [10]; (ii) the Hubbard model [11]; and (iii) the periodic Anderson model [12].We need to make explicit choices for the many-body Hamiltonians because the heat-current operator depends on the specific form of the potential energy, and is different for the three different models.
The Hamiltonians for the Falicov-Kimball and Hubbard models can be summarized by the following form: In equation ( 3), we have Fermionic creation (annihilation) operators c † i (c i ) for conduction electrons and f † i (f i ) for localized electrons at site i.The terms with µ and µ f are the chemical potential terms for each species of electron, and U is the on-site Coulomb repulsion.In the Falicov-Kimball model, we have t f = 0 and the f -electrons are localized on the lattice.In the Hubbard model, we have t = t f and we identify the c electrons as the up-spin electrons and the f electrons as the down-spin electrons.Then µ = µ f for the paramagnetic phase (no external magnetic field).Note that the generalized Hamiltonian in equation (3) interpolates continuously between the Falicov-Kimball and Hubbard models, and if we choose different chemical potentials, then it also includes the Zeeman term in a magnetic field.Our results will also hold for this generalized model.
The periodic Anderson model is written in a standard form, where the hybridization is local and independent of the lattice site i (this is not the most general form, as the hybridization can be local, but depend on the lattice site i, or the hybridization can be between orbitals on neighboring sites-we do not include such generalizations here) with the same conduction and localized electrons as in the Falicov-Kimball model, except now they have spin (E f is the localized electron site energy).The parameter V is the hybridization between conduction and localized electrons on the same site, and U is the on-site Coulomb repulsion between f -electrons with different spin.
We introduce an external electric field E(r, t) by employing a vector potential A(r, t) in the Hamiltonian or temporal gauge (where the scalar potential vanishes): Here we use units such that the speed of light satisfies c = 1 and Planck's constant satisfies = 1.We shall consider spatially uniform fields only, hence A is independent of r.Technically speaking this is only possible for a field that has no time dependence, because a time-varying electric field always yields a time-varying magnetic field via Maxwell's equations (the magnetic field is proportional to the curl of A in three dimensions).But if the time-variation is slow enough, then we can neglect the magnetic-field effects, which we do here.If desired, the small magnetic-field effects can be added back later via a gradient-based approach but we do not discuss such ideas further here.The electric field is coupled to the electrons via the Peierls substitution [13,14], which involves modifying the hopping matrix elements by a phase that depends on the line integral of the vector potential: and similarly for the t f ij term (where appropriate).Note that the Hamiltonian in the presence of an external field, H(A), is identical in form to that shown in equations ( 3) and (4), except it uses the hopping matrices in equation (6).Note also that t f = 0 for the Falicov-Kimball model.When the Hamiltonian gauge is used, translational symmetry is never broken, so the vector potential, electric field, and local charge density are all spatially uniform.
The "Peierls substituted" Hamiltonians become quite simple in momentum space: for the Falicov-Kimball-Hubbard model, where the Fermionic creation and annihilation operators now create or annihilate electrons with well-defined momentum.The periodic Anderson model becomes Since these Hamiltonians have essentially the same form as the case when there is no field, the derivation of the current operators is similar to the derivation for the linear-response operators; the only modification is that there is an additional time derivative that needs to be taken into account for the heat-current operator.We follow the procedure outlined in Mahan's textbook [15], which involves determining the charge-current operator by taking the commutator of the polarization operator with the Hamiltonian and the heat-current operator by taking the commutator of the energy-polarization operator with the Hamiltonian (plus the aforementioned partial derivative with respect to time).Since the commutator (plus partial derivative) can be identified with the total time derivative, the quantum statistical average of these current operators will satisfy the continuity equation.
We begin with the charge-current operator, which has an identical derivation as in the linearresponse case.We define the polarization operator to be Π pol = e i R i c † i c i , with the obvious generalization if we need to include spin (e is the charge of an electron and is less than zero).Since the potential energy for all Hamiltonians we consider is a function of the number operators for the various Fermions in the model, the polarization operator always commutes with the potential energy.So we need evaluate only the commutator with the kinetic-energy piece of the Hamiltonian (this does not hold if the hybridization depends on the lattice site i or if it occurs between different lattice sites in the periodic Anderson model, but we do not consider such situations here).The definition of the charge-current operator is which is the time derivative of the polarization operator.Evaluating the commutator is straightforward.If we take the result, evaluated in real space, and convert to a momentum-space representation, we find with the velocity operator defined by v(k) = ∇ (k).Equation ( 10) is generalized to include spin for the Hubbard and periodic Anderson model (for the generalized Falicov-Kimball-Hubbard model, we add a similar term with the f dispersion relation and the f operators).The heat-current operator is more complicated to construct.One first needs to define the energy-polarization operator, which requires assigning an energy operator h i to each lattice site i (H = i h i ).For a local potential energy, such an assignment is completely straightforward, but the hopping terms in the kinetic energy are more complicated -the standard procedure is to associate one half of the operator that connects sites i and j to the local energy operator at site i and at site j (this approach also automatically avoids double counting the kinetic energy contributions).Note that in the nonequilibrium case, the local energy operator h i is time-dependent, due to its dependence on the vector potential.The energy-polarization operator is then Π energy = i R i h i .The energy-current operator is defined to be and the heat-current operator follows as j Q = j E − µj/e.This is because heat is measured relative to the chemical potential of the material.Note that we need to include the partial time derivative because the energy polarization operator is time dependent.Since all of the time dependence is in the kinetic-energy piece of the energy polarization operator, that term will be the same for all models.We denote the two pieces of the heat-current operator as j Q = jQ + j Joule Q because the term coming from the time derivative is related to the Joule heat and is model-independent.The Jonson-Mahan theorem applies to the first term of the heat current.The expectation value of the second term (the Joule-heat-current operator) vanishes in linear response, because it is an order E 2 correction.
Unlike the charge-current operator, which was essentially independent of the model considered, the heat-current operator is different for different Hamiltonians.So we summarize the results for each case considered here.For the Falicov-Kimball model, the first term of the heat-current operator becomes And for the periodic Anderson model we have is the same as before and we only have to calculate This term becomes and in the Bloch representation it can be written as the symbol δ refers to the nearest neighbors of site i, so i + δ is a nearest neighbor site, and the symbol δ refers to the nearest neighbor translation vector.Putting these results together yields the final result for the first term of the heat-current operator of the periodic Anderson model jpam The Joule-heat-current operator can be written in the two equivalent forms for all of the Hamiltonians considered here where spin can be added trivially when necessary.Note that the Joule heat current is similar to the standard form described by a scalar potential multiplied by the local charge current and summed over all space for the scalar-potential gauge, which has the form of the top line in equation (19), It is possible that both forms yield identical results, but we have not established this.In order to evaluate the operator averages for the charge and heat current operators, we need to examine a formalism for the nonequilibrium many-body problem to complete the calculation.This is done in the next section.Before we move onto the next section, though, it is worthwhile to comment on the issue of Joule heat and its relation to the nonlinear heat current.Joule heating, determined by the I 2 R losses in a wire when current flows, is a heat created per unit volume in the presence of charge current flow.It arises from the fact that if we measure the heat current flowing through an area perpendicular to the direction of the flow, then the expectation value of the heat-current operator flowing through the area changes as we move in the direction of the heat-current flow.The change in the heat-current flow gives rise to the heat generated per unit volume in the system, which is equal to j • E. If we perform such an analysis for the Joule-heatcurrent term defined above, we find that analyzing the difference of the heat currents, between neighboring planes, yields precisely this Joule heating within the volume.This is why we gave it the name Joule heat current.

Nonequilibrium formalism for the many-body problem
Since the nonlinear Peltier effect and the nonequilibrium Jonson-Mahan theorem can be examined without explicitly solving the full many-body problem, we do not need to develop all of the nonequilibrium many-body formalism.Instead, we focus on the main issues needed to understand the relation between the charge and heat currents for nonlinear response.We imagine our system to be prepared in an equilibrium state, initially, with a temperature T = 1/β.Hence, the quantum-mechanical statistical averages will be with respect to the equilibrium Hamiltonian H(0) (the vector potential A(t) = 0 for early times, so H(A) = H(0) for those times).We use the symbols to denote the quantum statistical average [weighted trace over the states of H(0)] of an operator denoted by the dots with Z the partition function.In order to evaluate the current averages, we will need to use the so-called momentum-dependent lesser Green's function G < k , which is defined to be with the time-dependent Fermionic operators expressed in the Heisenberg picture, where all time dependence is in the operators, and the states are time independent (note that the lesser Green's function has no time ordering in its definition).These time-dependent Fermionic operators can also be represented in terms of Schroedinger operators, which only have explicit time-dependence in the parameters that multiply operators, since the quantum-mechanical states now have time dependence, via the connection for a (possibly time-dependent) operator O S (t) in the Schroedinger representation.The evolution operator U (t, t 0 ) is a time-ordered product, which satisfies the differential equation with t 0 being a reference time, which can be any time occurring before the external field is turned on, i. e., we must have A(t 0 ) = 0.The Heisenberg representation for the Hamiltonian is note that all Hamiltonians and current operators derived in the previous section are in the Schroedinger representation.
The lesser Green's function depends on two times.There is no time-translation invariance for the system anymore, because the field is turned on at some definite time.Rather than use the two times t and t in the definition of the lesser Green's function, it turns out to be much more useful to express results in terms of two other times, the average time T ave = (t + t )/2 and the relative time t rel = t − t , which are called the Wigner coordinates [16].We shall use the notation G < k (T ave , t rel ) to denote the Green's function expressed in Wigner coordinates, which is equal to G < k (T ave + t rel /2, T ave − t rel /2) when expressed in the original variables as in equation (21).
Using this new notation allows us to immediately learn that Hence, the expectation value for the charge-current operator can be expressed as in terms of the lesser Green's function.In the case where two spin species can conduct electricity, this result is generalized to include a summation over spin.Our goal, for the nonequilibrium Jonson-Mahan theorem, is to determine a similar formula for the expectation value of the first term of the heat-current operator.At first glance this looks to be hopeless, because the operator for the heat current contains two-particle operators (from the potential-energy terms), and it does not seem possible to relate it to the single-particle operators of the lesser Green's function.But the situation is not hopeless, as we quickly learn when we recall that a derivative of the Green's function with respect to t rel , brings down commutators of the creation and annihilation operators with the Hamiltonian.This analysis was carried out for both the first and second derivatives in order to determine the first two moments of the spectral function in [17].It is also related to the fact that one can determine the expectation value of the Hamiltonian from the Green's function, because of the special form for the equation of motion of the Green's function and how it relates to the Hamiltonian.
Nevertheless, the derivation of a similar formula for the heat-current operator requires a number of nontrivial steps, which we now describe in detail.We start by taking the derivative of the lesser Green's function with respect to t rel in the limit (27) The time derivatives are evaluated by using the Heisenberg equation of motion for any time-independent Schroedinger operator O S , with H(t) the total Hamiltonian of the system; one has a choice for which picture makes the most sense to use for a given calculation.The commutators are straightforward to calculate in either representation, and become for the Falicov-Kimball-Hubbard model in the Heisenberg representation and for the periodic Anderson model in the Heisenberg representation.Note that [c † k (t), H U (t)] = 0 for the periodic Anderson model.The commutator for the c operators is similar to evaluate.Putting these results together for the derivative of the lesser Green's function yields for the Falicov-Kimball-Hubbard model.For the periodic Anderson model, the U -term in equation (31) is replaced by a hybridization term, and we have a summation over spins, which yields Now, we multiply equations ( 31) and (32) by −v(k − eA(T ave )) and sum over k to find the expectation value for the first term of the heat-current operator.Hence we discover jQ for all models considered here (with the generalization to a sum over spin where needed).We believe this form will hold for all models with local interactions, but we offer no such general proof.The Joule-heat-current term is much more difficult to try to calculate.We have not succeeded in finding a way to fully evaluate it yet.The reason why one needs to be careful in evaluating the polarization-like objects is because the elements in the sum will diverge in the thermodynamic limit (in the real-space representation), due to the R i terms.So the evaluation of the expressions requires a thoughtful analysis.To begin, consider the charge-polarization operator, which can be expressed as follows: The expression in momentum space appears to be well-defined and feasible to calculate, but we are not aware of any simple way of evaluating the expectation values of gradients of creation or annihilation operators.If we examine the operator average in the Heisenberg picture, we find its time derivative satisfies where the velocity term arises from the kinetic-energy piece of the Hamiltonian, which commutes to the bandstructure times the creation or annihilation operator, and then produces the velocity when the gradient is taken (V is the potential-energy and hybridization pieces of the Hamiltonian).
To be concrete, we write down the result for the Falicov-Kimball model, which yields Evaluating the term proportional to the velocity is straightforward, but the other term appears to be quite complex, and may not be capable of being determined from just the Green's functions and the self-energies.This expression needs to be integrated over time to yield the desired expectation value.Naively, we expect the expectation value to vanish before the time when the field was turned on, but even that is not so obvious.The potential-energy terms are complicated in general, but if they are summed over k they vanish for the Falicov-Kimball-Hubbard model and for the periodic Anderson model.Hence, if we sum equation (36) over momentum, we immediately can derive the relation that d Π pol (t) /dt = j(t) , as we must.
If we try to proceed similarly with the derivation of a formula for the Joule heat current, we find that the procedure does not appear to work, but it is not obvious where it fails.To begin, if we take the expectation value of equation ( 19), we find If we take equation ( 35), integrate it over t, and substitute into equation (37), it seems like we should be able to derive an explicit formula for the Joule heat current, which will depend explicitly on the different models due to the commutators with the potential-energy piece of the Hamiltonian.But we do not trust this procedure because it gives a nonzero Peltier coefficient for the noninteracting particle-hole symmetric case (see below), and it also yields a contribution that is linear in the field, when we expect the Joule heat current to be at least quadratic in the field (otherwise this term would need to be included in the linear-response formalism, and it is not).We have not been able to resolve where the possible error comes from.Furthermore, in the interacting case, it appears to be quite a challenge to try to evaluate the terms that arise from commutators with the potential energy, because they still have gradients of the creation and annihilation operators in them, which we do not know how to evaluate.An alternative way to proceed is to find a way to clearly articulate how the schematic equation ∇ holds, and then integrate the Joule heat production per unit volume to yield the Joule heat current, but such a procedure does not appear simple to carry out either.

Relation to the equilibrium Jonson-Mahan theorem
The results in equations ( 26) and (33) are the nonequilibrium extension of the Jonson-Mahan theorem, and they allow one to determine part of the nonlinear Peltier effect for the response of a many-body system to a large electric field.Before we discuss such issues further though, it makes sense to examine how the nonequilibrium Jonson-Mahan theorem relates to the well-known (equilibrium) linear-response result.In this case, we explicitly drop the Joule heat current piece, because it should be quadratic in the field.
The first step is to perform a Fourier transform with respect to the t rel coordinates, so that Expressed in terms of the Fourier transform, we find the two current expectation values to satisfy because the Joule heat current term vanishes in linear response.This is then in the form of the Jonson-Mahan theorem, because the Peltier coefficient satisfies the form of equation ( 2), but one cannot immediately extract the result for τ (ω) without using the quantum Boltzmann equation to determine the change in the lesser Green's function to first order in the field.This change is proportional to the derivative of the Fermi-Dirac distribution, which then allows for τ (ω) to be extracted.Interested readers can consult [18] for a general derivation of the quantum Boltzmann equation in the continuum and [19] for a derivation on the lattice in the dynamical mean-field theory limit.Of course, the quantum Boltzmann equation approach is completely equivalent to the Kubo linear-response approach in the limit of small fields.

Nonlinear Peltier effect
Before we move onto a specific many-body example, it is worthwhile discussing these results in general first.In the noninteracting limit, the nonlinear current response has Bloch oscillations, which have a period inversely proportional to the strength of the electric field, and with an amplitude that varies with the temperature [20,21].Although we are unaware of any studies on nonequilibrium heat transport on a lattice, the heat current should also have Bloch oscillations in such a circumstance, because of the nonequilibrium extension to the Jonson-Mahan theorem.What is less clear is whether these oscillations are exactly in phase or not.If they are not in phase, then the Peltier coefficient will diverge when the charge current vanishes, but the Peltier coefficient could be a constant (or have more complex behavior) if the oscillations are in phase.
We can directly examine the nonlinear Peltier effect for the noninteracting system, which is known to have Bloch oscillations for the charge current.Instead of performing a completely general analysis, we shall focus here on the case of a hypercubic lattice in d-dimensions which has the vector potential oriented along the main (1, 1, 1, . ..) diagonal.Then the Peierls-substituted band structure becomes for the vector potential A(t) = (A(t), A(t), A(t), . ..).We define the ordinary band energy as = −2t i cos k i and the associated band energy as ¯ = −2t i sin k i .Then we study the joint density of states ρ( , ¯ ) for these two energies with δ(x) the Dirac delta function.Note that by shifting the wavevector the summation over the Brillouin zone, we can immediately prove the following symmetries for the joint DOS: The joint density of states is a completely even and symmetric function of and ¯ .This allows us to determine the general form for the nonlinear Peltier effect of a noninteracting system on a hypercubic lattice in d-dimensions.From the definition of j i (t) and G < k , equations ( 25) and (26), we obtain for each component.Since the joint density of states is symmetric in ¯ , the first term vanishes, and we find for each component that The first term of the heat-current operator is also easy to evaluate.Its components satisfy Using the symmetry of the joint density of states again, allows us to simplify this result to such that the first contribution to the nonlinear Peltier coefficient, defined by the ratio jQ (t) / j(t) , becomes the constant C satisfies This has the Heikes form [22,23] for the first term, plus an additional temperature-dependent correction.In the limit T → ∞, the numerator and denominator in equation ( 52) vanish, so we need to evaluate C via l'Hôpital's rule, which gives so that the high-temperature limit Heikes formula holds for the first contribution to the Peltier coefficient, as it must.For finite T , if the correction term is large enough, then the Peltier effect can change sign as a function of time, but if the magnitude of the constant is smaller than the magnitude of µ, it will simply oscillate, but not change sign.Note that at the particle-hole symmetric point of half-filling, where µ = 0, one can directly show that C = 0, so that the first contribution to the Peltier coefficient vanishes, as it must.
We can try the same type of analysis for the Joule heat current term.We work with a constant field, described by a linear vector potential turned on a time t = 0 and given by A(t) = −Et.Then, using equations (36) and (37), yields Performing the integral over t and using the symmetry of the joint density of states gives The term with a cos and sin adds to the first term of the heat current.But there are two problems with this result: (i) the expectation value has a term that is linear in E (the ¯ 2 term) and (ii) the expectation value does not vanish when µ = 0 (the sin term), which is the particle-hole symmetric case.This implies there is some error in the procedure used to evaluate the Joule heat current term, but we are unable to find it.Hence, we concentrate just on the first term of the heat current operator for the remainder of the manuscript.It is interesting to examine numerical values for the nonlinear Peltier effect, and to understand the time dependence for different cases.The simplest case to examine is the infinite-dimensional case with dynamical mean-field theory.Here, the joint density of states can be easily calculated to satisfy and then one finds that the coefficient C is similar in size to, but always slightly smaller in magnitude than µ, so that the first contribution to the nonlinear Peltier coefficient never vanishes when off of half filling, but the heat current becomes quite small in the vicinity of the zeros of the charge current Bloch oscillations [where cos eA(t) = sgn(µ)].We illustrate this behavior in figure 1, where we show the charge current, first term of the heat current, and first contribution to the nonlinear Peltier coefficient as functions of time, for a constant electric field turned on in the distant past [A(t) = −Et].The time axis is plotted in dimensionless units eEt, the filling is ρ e ≈ 0.251, and the temperatures are T = 0.1 and T = 1.As the temperature is increased, the amplitude of the Bloch oscillations decreases, and the first contribution to the Peltier coefficient approaches a more constant value, as expected due to the Heikes limit.When scattering is turned on, the general belief is that the Bloch oscillations will disappear due to the scattering in the system.Even if this is the case, there will be transient oscillations if the field is large enough, and so the question about whether oscillations are in phase or not becomes important again, at least in the short-time regime.
We test these ideas by using nonequilibrium dynamical mean-field theory to solve for the Green's functions, and then use the Jonson-Mahan theorem to calculate the expectation values for both the charge current and first term of the heat-current for the system.Details of the formalism and of the numerical calculations will be given elsewhere (preliminary discussions can be found in [17,24]).But we can briefly describe the numerical procedure.Like in the noninteracting case above, we take the vector potential to point along the hypercube diagonal in the d → ∞ limit.Then the summation over the Brillouin zone becomes a two-dimensional integral over the joint density of states.The numerical challenge arises from the fact that the momentum-dependent Green's function is represented by a general complex matrix in a real-time formalism, with time values on the Kadanoff-Baym-Keldysh contour.To generate the momentum-dependent Green's function for each quadrature point in the two-dimensional integration, we need to perform one matrix inversion and one matrix multiplication.Even if we use a Gaussian integration scheme for each dimension with, on the order of N = 100 points per dimension, we still need to evaluate on the order of 10,000 matrix inversions per iteration; typically 10-30 iterations are needed for reasonable convergence.
For our calculations, we take a time range of 40 units (taken in the limit were the reduced hopping t * = 1), and we turn the field on after 5 units of time.The electric field has a magnitude of 1 for each component, and the temperature (before the field is turned on) is 0.1.We use time steps of 0.067, 0.05, and 0.04, and extrapolate to the limit ∆t → 0 by using a pointwise quadratic Lagrange interpolation formula for the extrapolation.We take an average of a N = 54 and N = 55 Gaussian integration scheme for the integration in each dimension over the joint density of states.We find that the results for the charge current are much more accurate than for the first term of the heat current -the extrapolated charge current is truly zero for times before the time the field is turned on, but the extrapolated heat current is a nonzero constant, that heads towards zero as the number of points in the Gaussian integration scheme is increased, but is still nonzero even for N = 100 and N = 101.
The accuracy of the calculations can be tested by comparing with moment sum rules for the lowest few moments of the Green's function, as detailed in [17].The zeroth, first, and second moments of the local retarded Green's function can be found from the interaction, the chemical potential, and the fillings.In our case, the chemical potential is µ = −0.150868, the zeroth moment is 1, the first moment is 0.5259, and the second moment is 0.8234.Only the zeroth moment can be found from the local lesser Green's function without knowledge of more complex operator expectation values, which are beyond the scope of this work.This zeroth moment is equal to 0.25.When checking the retarded moments versus the scaled numerical data, we find the following agreement: all scaled moments are accurate to better than 0.5% once the field has been turned on.The accuracy is significantly worse in the equilibrium case, before the field is turned on (as large as 2-3%).The accuracy is less than 0.2% for times well after the field has been turned on.In general, the moments are less accurate for the times before the field is turned on; it turns out that the real-time formalism is generically more accurate for the nonequilibrium Green's functions rather than the equilibrium ones, but it isn't clear exactly why this is so.(c) the first contribution to the nonlinear Peltier effect for an interacting system on an infinitedimensional hypercubic lattice.The filling is ρe ≈ 0.25 and the temperature is T = 0.1.The Falicov-Kimball model interaction strength is U = 0.5, and the density of localized particles is ρ f = 0.75.Note the similarities of the currents with the noninteracting case, but that the amplitude of the oscillations decreases as time increases.The large values of the first term of the Peltier coefficient near the points where the charge current goes to zero are probably an artifact of the numerics not having enough accuracy to show identical zeros for the charge and heat transport.But we cannot confirm this with our limited computational resources.
Our results are shown in figure 2 for the case ρ e = 0.25 and ρ f = 0.75 and the spinless Falicov-Kimball model.The screened Coulomb interaction is U = 0.5 which corresponds to a strongly scattering metal (the metal-insulator transition occurs at U ≈ 1.8).Note how the oscillations in the currents have the same shape and a similar magnitude as the noninteracting case (the initial charge current is almost identical to the noninteracting charge current, the initial heat current is about a factor of 2 larger, which may be related to the fact that the chemical potential is changed dramatically due to the interactions, which can have a serious effect on the heat current), but the amplitude of the oscillations now decay as time increases; we cannot determine whether the steady state has oscillations or becomes a constant with the computer power that we have available to us.The nonlinear Peltier coefficient has strong additional oscillations in it, that are most likely a numerical artifact of the charge and heat currents not vanishing at the same times (the charge current vanishes at approximately integer multiples of π for E = 1).But even ignoring those "glitches" the data shows much richer behavior than the noninteracting result.These results for the first contribution to the nonlinear Peltier coefficient can have much more severe error in them due to the numerics.
We are unable to yet determine the Joule heat current contribution, which will certainly modify the overall heat current carried by the system.There does not appear to be any clear way of calculating those results at this point.
Note that both currents appear to oscillate around zero, rather than approach a steady-state limit.Because our calculations are limited to a finite time window, we cannot extend them beyond the current times without needing significantly more computational time.From these results, it is unclear to us exactly what the steady-state limit will look like.We cannot rule out the possibility of there being oscillations in both the charge and heat currents at long times.If the system instead settles down to a constant steady-state response, we cannot determine whether it goes to a finite limiting value or whether it goes to zero.The point is that once the system has built up a large number of nonlinear oscillations, intuition from the linear-response Kubo results need not apply any more for the steady state, so there are a number of different possible scenarios that can occur.We hope this work will interest future researchers into trying to answer just such questions.

Conclusions
In this work we have shown that there is a simple nonequilibrium generalization of the Jonson-Mahan theorem for Hamiltonians that have local potential energies.The heat current involves the Jonson-Mahan contribution plus an additional term related to the heat current; the latter piece does not appear to be capable of being evaluated in any simple way.Hence, there is a simple relationship between the charge and heat transport of the (nonlinear) Peltier effect, which involves the response of a material, kept at a fixed temperature, to a (possibly large) electric field.This relationship provides a deep connection between the charge and heat transport of strongly correlated materials.It holds in any dimension, and we conjecture for any model with a local potential energy; we explicitly demonstrated the result for the Falicov-Kimball, Hubbard, and periodic Anderson models.In addition, we found an extra nonlinear contribution to the heat current which we call the Joule heat current.This term is difficult to evaluate explicitly for any system.
One might ask whether these results can be generalized to include nonlinear transport effects in systems that have large temperature differences between different locations on the sample.This may be feasible, as it certainly holds in the linear-response limit, but much more theoretical work needs to be invested on developing the formalism for a nonequilibrium response in the presence of large temperature gradients, and this has not yet been attempted by anyone (including ourselves).However, the simplicity of the final result for the Peltier effect leads us to believe that there should be similar results for the Seebeck effect and for the electronic contribution to the thermal conductivity (but likely modified by other additional nonlinear terms as well).We do not try to provide such generalizations here.

Figure 1 .
Figure 1.Bloch oscillations of the (a) charge current, (b) first term of the heat current, and(c) first contribution to the nonlinear Peltier effect for a noninteracting system on an infinitedimensional hypercubic lattice.The filling is ρe ≈ 0.251 and the temperature is T = 0.1 (solid lines) and T = 1.0 (dotted lines).Note how the Peltier coefficient becomes more constant as the temperature is increased, ultimately assuming the Heikes form, and how the first term of the heat current does not appear to be modified as much by the temperature as the charge current.

Figure 2 .
Figure 2. Bloch oscillations of the (a) charge current, (b) first term of the heat current, and(c) the first contribution to the nonlinear Peltier effect for an interacting system on an infinitedimensional hypercubic lattice.The filling is ρe ≈ 0.25 and the temperature is T = 0.1.The Falicov-Kimball model interaction strength is U = 0.5, and the density of localized particles is ρ f = 0.75.Note the similarities of the currents with the noninteracting case, but that the amplitude of the oscillations decreases as time increases.The large values of the first term of the Peltier coefficient near the points where the charge current goes to zero are probably an artifact of the numerics not having enough accuracy to show identical zeros for the charge and heat transport.But we cannot confirm this with our limited computational resources.