Microscopic formulas for thermoelectric transport coefficients in lattice systems

A macroscopic description of thermoelectric phenomena involves several tensorial transport coefficients. Textbook microscopic Kubo formulas for them are plagued with ambiguities in the definitions of the current operators and the magnetization. We derive a version of these formulas for lattice systems which is free from ambiguities but contains additional terms compared to the textbook results. For symmetric components of thermoelectric tensors, we identify a large class of lattice systems for which the additional terms vanish with a natural choice of the energy current. To eliminate ambiguities in the skew-symmetric components, one needs to interpret them as relative quantities: only their differences for pairs of materials are well-defined.


I. INTRODUCTION
Thermoelectric effects have many scientific and technological applications [1]. They can also serve as probes of novel materials. Thus it is important to develop a theoretical framework for computing thermoelectric coefficients in the most general setting, including strongly interacting materials without well-defined quasi-particles.
Traditionally, the starting point for microscopic transport theory is provided by Kubo formulas. These formulas express transport coefficients in terms of correlators of volumeaveraged current densities of conserved quantities. But although Kubo formulas go back to [2,3] and can be found in many textbooks and monographs [4,5], there are a number of subtleties in their derivation. It is well appreciated by now that naive Kubo formulas for skew-symmetric parts of the transport tensors must be supplemented with additional terms involving magnetization and "energy magnetization" [6]. Such terms affect thermal Hall conductivity and the skew-symmetric parts of thermoelectric coefficients. Since magnetizations are intrinsically ambiguous, it is not obvious how to evaluate such terms, see [7] for a thorough discussion of magnetizations in general, [8,9] for the semiclassical case, and [10,11] for geometric approaches to defining energy magnetization. Another rarely discussed issue is the ambiguity in the definition of the energy density. One might expect that transport coefficients, being measurable quantities, are not affected by this ambiguity, but as far as we know this has been demonstrated only for the thermal conductivity and only for a special class of systems [12].
The theory of transport coefficients developed in [3,6,7] applies to continuum systems.
It cannot be directly applied to lattice systems because it assumes certain scaling relations for electric and energy currents which do not hold on a lattice. (In fact, they do not hold for interacting continuum systems either, except after some spatial averaging [6]). An even more basic issue is the lack of an accepted definition of charge and energy current densities on a lattice. Many expositions of linearized transport theory (see e.g. [4,5]) derive only the expressions for the volume-averaged current densities. But in order to define transport coefficients one needs to separate currents into transport and magnetization contributions [6,7]. Such a separation does not make sense for volume-averaged currents.
Since tight-binding models and other lattice Hamiltonians are ubiquitous in theoretical condensed matter physics, it is important to develop a formalism for describing currents of conserved quantities in such systems. In fact, such a formalism has been described by A. Kitaev many years ago [13], but it is rarely applied to transport theory. In our recent work we used it to prove a Bloch theorem for energy currents [14] and to derive Kubo-type formulas for the electric Hall conductivity and thermal Hall conductivity of general lattice systems [15]. In this paper use the same approach to derive microscopic formulas for thermoelectric coefficients of general lattice systems.
The main results of the paper are as follows. Our formulas for the symmetric parts of conductivity and thermal conductivity tensors [15] are completely analogous to continuum formulas. Surprisingly, this not the case for the symmetric parts of thermoelectric tensors. In general, microscopic formulas for them contain local terms as well as the expected Kubo term. We show that these extra terms are in fact required to ensure that transport coefficients are unaffected by the ambiguities in the definition of the microscopic energy density. We also show that in special cases, such as systems of free particles or systems with only density-dependent interactions, the additional terms vanish with a natural definition of currents.
The skew-symmetric parts of all transport tensors except conductivity contain contributions from magnetizations. Since magnetizations are defined only up to additive constants, this leads to ambiguities. In the case of thermal conductivity, a way to resolve the ambiguities on a lattice was described in [15] (building on the results of [6,7]), and the same approach works for thermoelectric tensors. Namely, although skew-symmetric parts of these tensors are "contaminated" with edge effects, ambiguities cancel when one considers differences of transport tensors for two materials. We express this by saying that skew-symmetric tensors are relative transport coefficients. Microscopic formulas for relative transport coefficients take a more complicated form: they are integrals of differential 1-forms along a path in the space of parameters. These issues do not affect the skew-symmetric part of the conductivity tensor because one can, in principle, measure it in a torus geometry, where no boundaries are present. This is not possible to do even in principle for the skew-symmetric parts of other transport coefficients.
The content of the paper is as follows. In Section II, we explain how ambiguities in definition of transport and magnetization currents leads to a natural separation of transport coefficients into absolute and relative ones. In Section III, we derive microscopic formulas for thermoelectric transport coefficients. We end with a discussion of possible generalizations of our results in Section IV. In one of the appendices, we specialize our formulas to the case of non-interacting fermions and express thermoelectric coefficients in terms of zero-temperature 1-particle Green's functions in coordinate space.
The work was supported in part by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award Number DE-SC0011632. A. K. was also supported by the Simons Investigator Award.

II. RELATIVE AND ABSOLUTE TRANSPORT COEFFICIENTS
The total current densities are usually divided into two parts: where magnetization currents are by definition divergence-free vector fields which do not contribute to net currents across any section of the system. Therefore they must have the where M N,E are defined by these equations and are usually called magnetization density and energy magnetization density, respectively. In the following we will use the same term magnetization for both magnetization and magnetization density which should not lead to a confusion.
The magnetization currents can be present even in an equilibrium state. The transport currents j N,E , on the other hand, can be present only in a non-equilibrium steady state created by slowly varying gradients of external electric potential and temperature (for simplicity of presentation we assume the chemical potential to be constant). This follows from the Bloch theorem [16,17] and its energy analogue [14]. This constrains the form of transport currents.
Further constraints arise from gauge-invariance. It requires the transport electric current j N and the transport heat current 1 j E − (ϕ + µ)j N to be invariant under shifts of the 1 The correction terms in the heat current originate from dN contributions to the heat T dS = dE − µdN − ϕdN , where the last term represents the work done by electromagnetic field. electrical potential ϕ by a constant. This follows from the equations defining the current and fact that under a constant gauge transformation ϕ(r) → ϕ(r) + c the energy density also transforms as h(r) → h(r) + cρ(r), where ρ is the electric charge density.
Taking all these considerations into account, one finds that to leading order in the derivative expansion the transport electric current is given by where the conductivity tensor σ km and the thermoelectric tensor ν km are functions of temperature only. For the energy current the expansion is A crucial point for this paper is that the separation of the current densities in (1) is ambiguous. One can always remove a curl of a vector field from j N,E and add it to j N,E mag without affecting the conservation equation and the form of the transport equations (5,6).
While this should have no effect on physically observable quantities, it can affect the transport coefficients. Let us specialize to the 2d case and decompose all two-index tensors into symmetric and anti-symmetric parts: σ km = σ S km + km σ A and similarly for the tensors ν, η, and κ. Taking into account the requirement of gauge-invariance, the allowed redefinitions of the transport currents have the form Here σ 0 is a constant and f (T ), g(T ) are arbitrary functions of T . Simultaneously magnetizations are redefined as follows: After the redefinition transport coefficients change: Using such a redefinition we can always set σ A (T ) to vanish at T = 0 and make κ A and η A vanish identically for any homogeneous material. Instead of setting η A (T ) = 0, one can choose to set ν A (T ) = 0 and use the remaining freedom to set η A (0) = 0. Note also that dσ A dT and ν A − dη A dT are invariant under such redefinitions.
So far we have ignored the vector potential, or equivalently gauge transformations which depend on the spatial coordinates. Allowing such gauge transformations changes the analysis as follows. Transport electric current j N and transport heat current j E − (ϕ + µ)j N are now required to depend on ϕ only through the electric field E k = −∂ k ϕ− ∂A k ∂t . The only difference this makes is that only transformations (7) with σ 0 = 0 are allowed. As a result, the Hall conductivity σ A is now free from ambiguities.
There is a natural way to fix ambiguities in M N and M E and therefore also in ν A , η A and κ A [6,15]. If we consider a material with a boundary, the magnetizations as well as all transport coefficients can be set to zero outside. This removes all ambiguities from transport coefficients, but obscures the fact that some transport coefficients are defined relative to vacuum, while others do no depend on any choices and can be measured in the bulk. We will call them relative and absolute transport coefficients, respectively. According to the above analysis, all symmetric transport coefficients as well as σ A are absolute, while This distinction has consequences for the microscopic formulas that can be derived for transport coefficients which are usually called Kubo formulas [15]. As we just explained, determination of relative transport coefficients require considering a system with boundaries.
On the other hand, as we show in the paper, derivatives of relative transport coefficients with respect to temperature or the parameters of the Hamiltonian involves only correlation functions of a system without boundary. The values of relative transport coefficients for any particular material can be found by integrating this differential over the parameters and/or temperature. The non-uniqueness in the choice of the base point of the integral reflects the ambiguity in the definition of the magnetization currents and can be fixed by choosing the base point to be a trivial insulator. The resulting microscopic formula for a relative transport coefficient is manifestly independent of the choice of boundary conditions at the cost of depending on the correlation functions of a whole family of systems which interpolates between the system of interest and a trivial insulator. On the other hand, microscopic formulas for absolute transport coefficients depend only on the correlation function of the system at a fixed temperature and values of all parameters.
As a consistency check, let us verify that physical bulk quantities depend only on absolute transport coefficients. For the time derivatives of charge and energy densities we get where the external fields are assumed to be time-independent and thus ∇ × E = 0. As expected, these time derivatives are unaffected by the transformations (11).
The time-derivative of entropy density is usually written as follows [18]: where the entropy current density is The r.h.s. of eq. (14) seems to depend on some relative transport coefficients. However, if one redefines the entropy current as follows: then one can write It is manifest now that both the entropy production rate and the entropy current depend only on the absolute transport coefficients.
Since only absolute transport coefficients enter the expressions for the divergences of currents, measuring net currents through closed curves (or surfaces, if we are discussing a 3d material) does not allow to determine relative transport coefficients. This applies even to infinite curves with a boundary at infinity, provided ϕ and T tend to fixed values at infinity. The latter condition must be imposed to eliminate the contribution of magnetization currents. For example, if we compute the electric current I N x through a vertical line x = 0, then the contribution of ν A drops out because The above considerations apply to a homogeneous material whose transport coefficients are constants. If one considers a heterogeneous material, such as an interface between two homogeneous ones, then the expressions for net currents will involve differences between relative transport coefficients. For example, consider a sample such that ν A interpolates between ν A 1 for y 0 and ν A 2 for y 0. Suppose that the temperature is a function of y only which is equal to T b throughout the interface region and approaches T ∞ at y → ±∞.
Then the ν A -dependent contribution to the net electric current in the x direction is Similarly, the contribution of κ A to the net heat current is ( . By creating an electric potential ϕ which is equal to ϕ b in the interface region and approaches ϕ ∞ at y → ±∞, one can also measure σ A 2 − σ A 1 and η A 2 − η A 1 . In the case of the electric Hall conductivity one can do better by utilizing a time-dependent vector potential rather than a scalar potential and working in a cylinder geometry or a torus geometry. Then in principle one can determine σ A for a single material by measuring the net flow of electric charge across a section of a cylinder or a torus as one inserts a unit of magnetic flux through this section. This does not work for η A because the physical quantity that needs to be measured is the net amount of heat transferred to the heat bath as one inserts a unit of magnetic flux (see Fig. 1). Therefore heat transport will receive a contribution from the work of the electromotive force E on the net electric edge currents I N edge . The edge currents are proportional to the jump in the magnetization along the boundary and they make η A relative even in the cylinder geometry.

A. Currents on a lattice
We follow the conventions of [15]. We consider a lattice system with a Hamiltonian H = p∈Λ H p , where Λ ⊂ R 2 is a not necessarily regular lattice. The operators H p have a finite range, i.e. there exist R such that H p acts trivially on site q if |p − q| > R. The space of states at each site of the lattice is assumed to be finite-dimensional. The electric charge operator Q has the form Q = p∈Λ Q p , where Q p has integral eigenvalues (we set the electric charge of electron to be 1) and acts only on site p. This means that the U (1) symmetry is on-site. In particular, [Q p , Q q ] = 0 for all p, q ∈ Λ.
The electric current from site q to site p is defined as The energy current from site q to site p is defined as These currents enter the charge and energy conservation equations The net current from a subset B ⊂ Λ to its complement A = Λ\B is given by This observable measures the net current across the boundary of A and B. In this paper we will need its mild generalization. For a given skew-symmetric function η(p, q) : This equation is a lattice analog of the continuum equation In the continuum the general solution to this equation defines the magnetization M N and energy magnetization M E . Analogously, on a lattice the solution of equation (27) is where M N,E pqr are skew-symmetric functions of the lattice points p, q, r ∈ Λ. These are lattice analogs of the magnetization and the energy magnetization. Physically, in continuum case M N,E (r) represent the circulating currents of the system in equilibrium. Similarly, M N,E pqr physically can be thought as quantity which measures the circulating current around a triangle formed by p, q, r.
Unfortunately, M N,E pqr is not unique: one can always redefine where N pqrs is skewsymmetric function of its subscripts which decays whenever any two of them are far apart. This corresponds to ambiguity in splitting of the circulating currents into contributions of magnetization from the different triangles and it is absent in continuum case.
There is an additional ambiguity corresponding to existence of solutions to the equation r∈Λ M N,E pqr = 0 which are not of the form s∈Λ N pqrs . It corresponds to an ambiguity of addition of a constant to the magnetization in continuum case. A standard method to deal with the later is to consider a system with a boundary and fix the magnetization to be zero outside of the system. In this paper, we want to think about all transport coefficients as manifestly bulk quantities and avoid dealing with boundaries. While magnetization itself suffers from ambiguities and depends non-locally on the boundary conditions, the variation of magnetization with respect to parameters of the Hamiltonian is local. Indeed, consider the variation of the equation (30) with respect to a parameter λ of the Hamiltonian where µ N,E pqr, = ∂M N,E pqr ∂λ and is given by [13] where A; B denotes the Kubo canonical pairing [2]. Using the properties of the Kubo pairing (see Appendix A), one can easily verify the identity (32). In the following we will combine derivatives of magnetizations with respect to different parameters into 1-forms on the parameter space µ N,E pqr = µ N,E pqr, dλ .

B. Equilibrium conditions and driving forces
In the following sections we will follow Luttinger [3] and study the behavior of the system coupled to external potentials where ϕ(p) is external electric potential and ψ(p) can be thought of as gravitational potential.
The potentials are assumed to infinitesimally small slowly-varying functions of p which vanish at infinity. After coupling to external potentials the system will eventually relax into a state with density matrix where T 0 and µ 0 are the temperature and local chemical potential of the system at infinity.
On the other hand, on physical grounds we expect local observables supported in some small but macroscopic region around site p to be described by a thermal density matrix where the local temperature T (p) and chemical potential µ(p) are slowly varying functions of p. The equilibrium conditions can be found to be [3,6,7] (1 + ψ(p))(µ(p) + ϕ(p)) = µ 0 , These relations together with the absence of transport currents in equilibrium can be used to derive Einstein relations between transport coefficients [3]. The latter also leads to transport currents being proportional to the gradients of the left hand sides of eqs. (37) and (38). The fact that the driving forces depend only on specific combinations of ψ, ϕ, T, µ will be used to relate the response to variations of the thermodynamic parameter T to the response to variations of the external field ψ.

C. Nernst effect
In order to find the thermoelectric coefficient coefficient ν xy we deform the Hamiltonian density by where g (p) is a hat-shaped function as in Fig. 2b, is an infinitesimal parameter, and s is a small positive number which controls how fast the perturbation is turned on. We will consider the so-called fast regime [3] in which the characteristic time 1/s is large but not large enough in order for the two slopes of the hat g (p) to come into equilibrium. The change of the state of the system can be found as follows. The density matrix satisfies the quantum Liouville equation where ρ 0 is the equilibrium density matrix at t = −∞ and dots represent higher order terms in . The solution to this equation is where the dot denotes the time derivative. The change of the observable A can be found to where we used the Kubo pairing notation (see Appendix A), β = 1 T is inverse temperature, and ∆A is the variation of the operator arising from explicit dependence of A on the Hamiltonian. Using explicit form of the perturbation (39), energy conservation law and properties of the Kubo pairing we can rewrite this formula as where we neglected term proportional to small s.
The change in the electric current across a vertical line x = a is where f (p) = θ(a − x(p)) is a step function. The explicit variation of the current is where in the last line we separated the result into two contribution formally skew-symmetric and symmetric in f, g . The second term depends only on difference of values of g at different sites as expected for a transport current. On the other hand, the first term depends on the value of g and does not seem to be of the form expected for a transport current. As was explained in [6] this contribution is related to magnetization currents. Indeed we can rewrite where we used skew-symmetry of M N pqr .
Written in this form, the response is proportional to the differences of g at different points and therefore receives appreciable contributions only from regions I and II in Fig. 2d.
However, the transformation (47) contains an important subtlety. The right-hand side hand side contains a magnetization contribution which is ambiguously defined while left-hand side is unambiguous. There is no contradiction because the ambiguity in the region I will compensate the one in the region II. Moreover, in a homogeneous system the response will be zero, because these two regions compensate each other exactly. One way to deal with it is to introduce a boundary somewhere in between the two regions and enforce the magnetization M to be 0 outside of the sample. This approach is used in [6] in the continuum case. However, an explicit boundary introduces additional computational challenges and makes the bulk nature of the Nernst effect obscure.
In this paper, we will use an alternative approach proposed in [15]. Instead of introducing a sharp boundary, we will make one parameter of the Hamiltonian λ to have slightly different values in regions I and II (see Fig. 2d). We can write the hat-shaped function g as a difference of two functions g I and g II , g (p) = g II (p) − g I (p), each of which is a translate of the function g(p) which depends only on y(p) and is shown in Fig. 2c. The functions g I,II are non-constant only in regions I and II respectively. Then the magnetization contribution can be rewritten as where µ N pqr,λ = ∂M N pqr ∂λ , and we introduced a notation Combining this with other contributions we find where The function g as in Fig. 2c is not compactly supported and thus it takes an infinite time for the system to equilibrate. Therefore, one can take the limit s → 0 while staying in the "fast" regime. Using the Einstein relation following from eqs. (37,38), one finds that electric current after perturbation by gravitational potential g(p) is equal to the current generated by From continuum phenomenological equation (5) one finds the current across x = a line to be Comparing it to (50) we find where we introduced the notation J Q = J E − µJ N for the heat current and we dropped the subscript 0 from T 0 and µ 0 since this formula contains correlation functions of the unperturbed system in equilibrium. We combined differential with respect to parameter into the differential form µ N (δf ∪δg) = µ N λ (δf ∪δg)dλ and the derivation can be straightforwardly extended to involve several parameters. The exterior derivative d = dλ ∂ ∂λ acts on the parameter space.
Since rescaling the temperature is equivalent to rescaling the Hamiltonian, we can extend this 1-form to the enlarged parameter space which includes T . Then we can define the difference of coefficients η xy for any two 2d materials, regardless of the temperature. Explicitly, let us define the rescaled Hamiltonian H λ 0 = λ 0 H, where we introduced an additional scaling parameter λ 0 . Then Therefore we can define the T -component of the 1-form on the enlarged parameter space as follows: where τ N (δf ∪ δg) is given by eq. (49) with µ pqr replaced with which is obtained from µ N by replacing dH p with −H p .

D. Ettingshausen effect
One can derive a formula for the coefficient η xy in a similar way. In this section, we will display only the key steps, since all arguments are the same.
In order to find the coefficient we deform the Hamiltonian density by where g (p) is a hat-shaped function of y(p) as in Fig. 2b.
The change in the energy current across a vertical line x = a is where f = θ(a − x(p)). The explicit variation of the current is The first term can be expressed in terms of magnetization as in (47). We write g (p) as a difference g (p) = g II (p) − g I (p), where g I,II are translates of a smeared step-function g(p).
Then we rewrite the response as a difference of conductivities of different materials: On the other hand, from the continuum phenomenological equation (6) one finds the current across the line x = a to be where the second term originates from the contribution of ϕj N to j E . Comparing it to (62) we find This 1-form can be extended to include temperature as a parameter using the scaling Therefore we can define the T -component of the 1-form on the enlarged parameter space as follows: d dT where τ N is given by (58).

E. Symmetric parts of transport coefficients
Note that the 1-forms µ N (δf ∪ δg) and τ N (δf ∪ δg) are formally skew-symmetric under the exchange of f and g. To make use of this symmetry, we need to argue that f can be replaced with a smeared step-function in the x direction. This can be argued using matching between microscopic theory and hydrodynamics. Namely, we expect that the microscopic linear response can be used to compute the properties of a Non-Equilibrium Steady State and f (p) is a hat-shaped function which depends only on x(p). On the other hand, if the microscopic linear response is to reproduce the expected properties of a NESS, the expectation value of these observables must be zero. Thus we can replace f with a smeared step-function in the x-direction without affecting ν xy or η xy .
After this has been done, exchanging f and g is equivalent to exchanging x and y. Thus As we show in Appendix C, the two terms on the right-hand side of these equation are not separately invariant under Hamiltonian density redefinition, but the full transport coefficients are invariant. The correction term U (δf, δg) is zero for fermionic systems with only density-dependent interactions (see Appendix D for more details).

F. Skew-symmetric parts of transport coefficients
In a similar way one can find skew-symmetric parts of transport coefficients These formulas give only derivatives of transport coefficients with respect to parameters.
Integration of these formulas over parameters or temperature gives the difference of relative transport coefficients at different values of parameters. It is natural to define relative transport coefficients of a trivial insulator to be zero. Determination of the relative transport coefficient in this case would correspond to integration over a path in the parameter space from a trivial insulator to the material of interest.
There are many paths which one can use to deform a system into a trivial one. Consistency requires the integral of dν A or dη A to depend only on the endpoints of the path. This means that the 1-form µ N (δf ∪ δg) must be exact. This can be proved using the techniques of [15] where µ E (δf ∪ δg) was shown to be exact.

IV. DISCUSSION
In this paper we have derived microscopic formulas for "transverse" thermoelectric coefficients of general 2d lattice systems. It was convenient to decompose them into symmetric and anti-symmetric parts, since they have qualitatively different behavior: the former are absolute transport coefficients, while the latter are relative. Similar formulas for electric Hall conductivity and thermal Hall conductivity have already been derived in [15].
The usual Kubo formulas for transport coefficients require averaging the correlators of currents over the whole space. In contrast, our microscopic formulas involve net currents through two perpendicular lines. This is because a current on a 2d lattice is a function of a pair of points, and the natural observable associated to it is localized on a line rather than at a point. Despite this, after the limit s → 0 has been taken, our formula computes the same quantity as the usual continuum Kubo formula.
It is natural to ask whether longitudinal components of conductivity, thermal conductivity, and thermoelectric tensors of a 2d lattice system can be computed in a similar manner. This is easily achieved: one simply replaces two perpendicular lines with two lines making a nonzero angle θ. It is easy to see determine from hydrodynamics which linear combination of components of the transport tensors describes the corresponding linear response. For example, if both currents involved are electric currents, one of the lines is given by y = 0, and the other one is y = x · tan θ, then the correlator measures σ xy − 1 tan θ σ yy . Thus by changing the functions f, g one can extract all four components of the conductivity tensor. The same is true about other transport coefficients.
In this paper we focused on the case of 2d materials, but the 3d case can be accommodated as well. One can simply replace a lattice in R 2 with a lattice in R 2 × [0, L], impose periodic boundary conditions in the third direction, divide all formulas by L, and take the limit L → ∞. The functions f, g remain independent of the third coordinate. It should not matter whether the limit L → ∞ is taken before or after the limit s → 0, since the problem is translationally invariant in the third direction.
In the 2d case, the quantity ν A (T ) (normalized relative to the vacuum) is dimensionless, and Bloch's theorem implies that for gapped systems ν A (0) does not change under the variations of the Hamiltonian which do not close the gap. Thus if ν A (0) were nonzero, it would represent a new topological invariant of gapped 2d phases of matter. However, one can show that on very general grounds ν A (0) vanishes for all gapped systems [19]. By Onsager reciprocity, the T → 0 limit of η A (T )/T also vanishes. Thus topological invariants of gapped 2d systems arise only from the Hall conductivity and the thermal Hall conductivity. of time-reversal-invariant 2d systems. For the skew-symmetric thermal conductivity κ A , where λ is a parameter of the Hamiltonian. Thus κ A can be a function of temperature only.
Further, if we treat T as a parameter, then scaling analysis gives Hence κ A (T ) = aT , where a does not depend on parameters. The parameter a has no physical significance, but it is natural to set it to zero, so that the vacuum has zero thermal Hall conductivity. Thus we reach the standard conclusion that for a system with time- The case of thermoelectric coefficients is slightly different. Usually one says that Onsager reciprocity requires ν km = T −1 η mk , which implies ν A + T −1 η A = 0 [18]. Since both ν A and η A are relative transport coefficients, one should interpret this as Hence ν A + T −1 η A can depend only on the temperature. If we treat temperature as a parameter, then the scaling analysis gives Hence ν A + T −1 η A = a, where a is a constant which is independent of any parameters or temperature and has no physical significance. One can choose it to be zero. Then So for a time-reversal-invariant 2d system there is only one independent skew-symmetric thermoelectric transport coefficient, namely ν A .

Appendix C: Invariance under Hamiltonian density redefinition
For a given Hamiltonian, there are many ways to define a Hamiltonian density. A typical example of this is the ambiguity in splitting an interaction term between two sites p and q into H p and/or H q . In this appendix, we will show that our microscopic formulas for physically observable transport coefficients are independent of the choice of the Hamiltonian density, even though individual terms in the microscopic formulas are not invariant. For some systems this can be used to simplify the microscopic formulas.

Invariance of the electric current
Consider the following change of the Hamiltonian density where A rp is skew-symmetric in r, p. We want the final Hamiltonian to be U (1)-invariant.
Therefore, we have to impose For a general choice of A pq a stronger condition will not hold. However, one can always redefine A pq (by subtracting the U (1)-non-invariant part) in such a way that (C3) holds without affecting H p . In the following we will assume this was done and (C3) is true.
Under the transformation (C1) the electric current changes as Even though the current density changes, the net current through any section is invariant. Indeed, and the last term is zero since where we have used (C3) and the skew-symmetry of [A rq , Q p ] + [A pr , Q q ] + [A qp , Q r ].

Covariance of the energy current
Let us now consider the effect of the redefinition of the Hamiltonian density on the energy current. Imposing an energy analog of (C2) or (C3) is far too restrictive, since it would only allow changes of the Hamiltoniain density by conserved quantities. For example, the difference between putting the interaction term between the two sites p and q either into H p or into H q corresponds to A pq equal to the interaction term. Obviously, interaction terms are not integrals of motion in general. Because of this we will not impose either of the equations in (C7).
Under the redefinition of the Hamiltonian density (C1) the energy current changes as while the net current transforms as The last term can be rewritten as where we have defined We find that the net energy current transforms as follows under a redefinition of the Hamiltonian density: But this should be expected since a redefinition of the energy density changes how we define the energy of sub-regions and therefore should affect the net energy current. Indeed, one can see that (C9) is exactly the transformation needed in order to satisfy the energy conservation for the new energy density H p + q∈Λ A qp . By summing this transformation law over p weighted by a function f (p) with a compact support we find thaṫ which reproduces (C11). Here we used an identity which is true for any f with a compact support.
From the above discussion, one can see that energy current is not invariant but covariant under energy density redefinitions. If we choose f (p) to be 1 when p is in some compact set B and zero otherwise, the physical meaning of (C11) is very clear. It corresponds to ambiguities in the energy currents due to interaction terms along the boundary of B. Depending on how we distribute the interaction terms among H p we can change the energy stored in the region B as well as energy current through its boundary.

Invariance of the microscopic formulas for thermoelectic coefficients
In this section we will show that the coefficients ν xy and η xy are invariant under a redefinition of the Hamiltonian density. We will start with skew-symmetric coefficients Here we defined the Kubo parts as Under Hamiltonian density redefinition the Kubo parts transform as where we used properties of the Kubo pairing.
Before finding the variation of the magnetization term it is useful to rewrite it slightly: Note that one cannot expand the square brackets, since the two resulting sums over p, q will not converge separately.
Let us find the variation of 1 2 J N pq (g(p) + g(q)) under a Hamiltonian density redefinition. It reads 1 2 J N pq (g(p) + g(q)) → 1 2 J N pq (g(p) + g(q)) The last term can be rewritten as follows: where we used the properties of the Kubo pairing. The first term in this expression can be rewritten as = − 1 2 s,r∈Λ J N sp (g(s) + g(p)) + J N ps (g(s) − g(p)); A rq − (p ↔ q) = J N (δg); A qp − 1 2 r,s∈Λ J N rp (g(r) + g(p)); A sq + J N ps (g(s) − g(p)); A rq + 2 perms , (C24) where "2 perms" means the two cyclic permutations in p, q, r. Note that the term in square brackets is skew-symmetric in p, q, r. The second term can be rewritten as Note that term in square brackets is skew-symmetric in p, q, r By combining equations (C21-C25) we find that the magnetization contribution changes under a redefinition of the Hamiltonian density as follows: where C pqr is a skew-symmetric function of p, q, r which is combination of skew-symmetric parts (and their derivatives) in the right-hand sides of Eqs. (C21-C25). Due to its skewsymmetry we find that 1 2 p,q,r∈Λ C pqr (f (q)−f (p)) = 1 6 p,q,r∈Λ The variation of Kubo parts were already determined before, so we focus on the transformation of U . Under (C1) it transforms as follows: We can rewrite this equation by noticing that Then using eqs. (C24, C25) we find We see that the variation of this term cancels the varitions of the Kubo parts.
One can do the same checks for the thermal Hall conductivity and verify that the microscopic formula derived in [15] is in invariant under a redefinition of the Hamiltonian density.
To linear order in A pq all the manipulations are almost the same except for the replacement The electric current can be found from the conservation equation: The net current through a section defined by δf (p, where a bounded function f ∈ 2 (Λ) is understood as an operator acting on the one-particle Hilbert space 2 (Λ) by multiplication. Summation over sites is implicit.
The energy current operator is The net energy current is The state of the system at a temperature T = 1/β is defined via Wick's theorem and Gibbs distribution a p (t)a † q = p e −iht 1 + e −βh q , (D7) where a p (t) are operators in the Heisenberg picture.
Using these formulas we find where the trace is over the 1-particle Hilbert space 2 (Λ), and the functions f : Λ → R and where ε n are 1-particle Hamiltonian energy eigenvalues.
Multiplying this by e −st and integrating over t, we arrive at where f(ε) = 1 1+e β(ε) is the Fermi-Dirac distribution. We absorb the chemical potential into a shift of the Hamiltonian.
The above expressions assume a discrete energy spectrum and thus can only be used for finite-volume systems. To get an expression applicable to infinite-volume systems, let us rewrite it in terms of the one-particle Green's functions G ± (z) = 1/(z − h ± i0). Some of the useful formulas are where we have dropped z for G ± (z). Here A and B are operators acting on the one-particle Hilbert space, and in the second formula we assumed in addition that their average is zero: a † Aa = a † Ba = 0. Note also that hG ± = G ± h = zG ± − 1, [G ± , A] = G ± [h, A]G ± .
Using this notation the formula for the electric conductivity takes the form where the integration is over the real axis in the z-plane.

Magnetization term
The magnetization differential for an arbitrary deformation dh of the 1-particle Hamiltonian is given by For temperature variations this expression can be simplified to These expressions are needed only for the evaluation of skew-symmetric parts of the thermoelectric coefficients.

U -term
Let us study the term (51) for free fermionic system. In this case the relevant many-body operators become where δ p is a Kronecker delta function equal 1 on site p and 0 on all other sites. A product of two delta functions enforces q = p in the summation over p and q. Since U also involves a factor of (g(p) − g(q))(f (p) − f (q)), U (δf, δg) vanishes for systems of free fermions.
More generally, one can consider a system of fermions with only density-dependent interactions. Namely, suppose we allow the following interaction term in the Hamitonian where V (p 1 , . . . , p n ) is a function of n sites which describes the potential energy of manybody interaction and decays rapidly when the points p 1 , . . . , p n are far from each other. One can see that this term will leave eq. (D13) unaffected since Q p commute with each other.
We conclude that for fermionic system with only density-dependent interactions there is no correction originating from U to the symmetric thermoelectric coefficients provided H p is chosen in the manner explained above.

Skew-symmetric part
Consider the variation of the Kubo parts (C17,C18) of the skew-symmetric thermoelectic coefficients under a rescaling of the Hamiltonian: dh = h dλ 0 . We get Summing up this contributions with the magnetization contribution gives Integrating over the temperature and using the formula where gives Here we normalized the thermoelectric coefficients to be 0 in the infinite-temperature state.
Note that since in the limit T → 0 the Fermi-Dirac distribution f(z) becomes a step-function,