\section{Introduction}

Hydrodynamics in a general sense is a coarse-grained description of matter in local thermal equilibrium which incorporates all local conservation laws. It depends on the microscopic details only through transport coefficients.
Dissipative transport coefficients such as viscosity, thermal conductivity and electric conductivity are particularly important for applications and their theory is better known. More recently a lot of attention has been devoted to non-dissipative  time-reversal-odd  transport coefficients such as electric Hall conductivity and Hall viscosity. This interest is explained in part by the contribution that gapless  edge modes make to some of these coefficients at low temperatures.

Theoretical studies of transport phenomena often take Kubo formulas as their starting point. These are microscopic formulas for transport coefficients in terms of correlators of currents of conserved quantities. Although Kubo formulas go back to \cite{Kubo,Luttinger} and can be found in many  textbooks, 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" \cite{CHR}.  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 \cite{Niuetal} for a thorough discussion.

The theory of transport coefficients developed in \cite{Luttinger,CHR,Niuetal} 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 \cite{CHR}). Since tight-binding models and other lattice Hamiltonians are ubiquitous in theoretical condensed matter physics, it is important to develop an analogous formalism for them. In \cite{thermalHallpaper} we described such a formalism and used it to derive microscopic formulas for the Hall conductivity and thermal Hall conductivity of lattice systems. In this paper we extend it to other transport  coefficients. For simplicity, we focus on the case of 2d lattice systems and only briefly discuss the 3d case.

The main results of the paper are as follows. Our formulas for the symmetric parts of conductivity and thermal conductivity tensors \cite{thermalHallpaper} 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 was described in \cite{thermalHallpaper} (building on the results of \cite{CHR}), 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 \ref{sec: hydro}, 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 \ref{sec: themoelectric Kubo}, we derive  microscopic formulas for thermoelectric transport coefficients. We end with a discussion of possible generalizations of our results in Section \ref{sec: thermoel discussion}. 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.

\section{Relative and absolute transport coefficients}\label{sec: hydro}
The total current densities are usually divided into two parts:
\begin{align}\label{eq: tot into mag + tr}
\begin{split}
    {\bf j}^N_{\rm tot} &=     {\bf j}^N +     {\bf j}^N_{\rm mag},\\
      {\bf j}^E_{\rm tot} &=     {\bf j}^E +     {\bf j}^E_{\rm mag},
\end{split}
\end{align}
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 form
\begin{align}
    {\bf j}^N_{\rm mag} &=     \nabla \times    {\bf M}^N,\\
      {\bf j}^E_{\rm mag} &=      \nabla \times    {\bf M}^E,
\end{align}
where ${\bf 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 ${\bf 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 \cite{Watanabe,Bloch} and its energy analogue \cite{energyBloch}. This constrains the form of transport currents.

Further constraints arise from gauge-invariance. It requires the transport  electric current ${\bf j}^N$ and the transport ``heat current'' ${\bf j}^E-\varphi {\bf j}^N$ to be invariant under shifts of the electrical potential $\varphi$ by a constant. This follows from the equations defining the current operators~${\bf J}^{N,E}({\bf r})$
\begin{equation}\label{eq: hydro currents def}
i[ H, h({\bf r})]=-\nabla\cdot {\bf J}^E({\bf r}),\quad  i[ H, \rho({\bf r})]=-\nabla\cdot {\bf J}^N({\bf r}),
\end{equation}
and fact that under a constant gauge transformation $\varphi({\bf r})\mapsto\varphi({\bf r})+c$ the energy density also transforms as $h({\bf r})\mapsto h({\bf r})+c \rho({\bf r})$, where $\rho$ 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
\begin{equation}\label{eq: jn hydro}
{ j}^N_k=-\sigma_{km}\partial_m\varphi -\nu_{km} \partial_m T,
\end{equation}
where the conductivity tensor $\sigma_{km}$ and the thermoelectric tensor $\nu_{km}$ are functions of temperature only. For the energy current the expansion is
\begin{equation}\label{eq: je hydro}
 { j}^E_k=\varphi { j}^N_k-\eta_{km}\partial_m\varphi -\kappa_{km} \partial_m T.
\end{equation}

A crucial point for this paper is that the separation of the current densities in (\ref{eq: tot into mag + tr}) is ambiguous. One can always remove a  curl of a vector field from ${\bf j}^{N,E}$ and add it to ${\bf j}^{N,E}_{\rm mag}$ without affecting the conservation equation  and the form of the transport equations (\ref{eq: jn hydro},\ref{eq: je hydro}).  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: $\sigma_{km}=\sigma^S_{km}+\epsilon_{km}\sigma^A$ and similarly for the tensors $\nu,\eta,$ and $\kappa$. Taking into account the requirement of gauge-invariance, the allowed redefinitions of the transport currents have the form
\begin{align}\label{currenttransf}
   j^N_k &\mapsto j^N_k + \epsilon_{km}\partial_m  \left(\sigma_0\varphi +  f(T)\right), \\
    j^E_k &\mapsto j^E_k + \epsilon_{km}\partial_m\left(\frac12\sigma_0 \varphi^2+ f(T) \varphi+ g(T)\right).
\end{align}
Here $\sigma_0$ is a constant and $f(T), g(T)$ are arbitrary functions of $T$. Simultaneously magnetizations are redefined as follows:
\begin{align}
M^N &\mapsto M^N-\sigma_0\varphi-f(T),\\
M^E &\mapsto M^E-\frac12\sigma_0\varphi^2-f(T)\varphi-g(T).
\end{align}
After the redefinition transport coefficients change:
\begin{align}\label{eq: antisym part shifts}
\begin{split}
\sigma^A&\mapsto\sigma^A-\sigma_0,\\
\nu^A&\mapsto\nu^A-\frac{df(T)}{dT},\\
\eta^A&\mapsto\eta^A-f(T), \\
\kappa^A&\mapsto \kappa^A - \frac{dg(T)}{dT}.
\end{split}
\end{align}
Using such a redefinition we can always set $\sigma^A(T)$ to vanish at $T=0$ and make $\kappa^A$ and $\eta^A$ vanish identically for any homogeneous material. Instead of setting  $\eta^A(T)=0$, one can choose to set $\nu^A(T)=0$ and use the remaining freedom to set $\eta^A(0)=0$. Note also that $\dfrac{d{\sigma}^A}{dT}$ and  $\nu^A-\dfrac{d\eta^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 ${\bf j}^N$ and transport heat current ${\bf j}^E-\varphi {\bf j}^N$ are now required to depend on $\varphi$ only through the electric field $E_k=-\partial_k\varphi-\frac{\partial A_k}{\partial t}$. The only difference this makes is that only transformations (\ref{currenttransf}) with $\sigma_0=0$ are allowed. As a result, the Hall conductivity $\sigma^A$ is now free from ambiguities.

There is a natural way to fix ambiguities in $M^N$ and $M^E$ and therefore also in $\nu^A,\eta^A$ and $\kappa^A$  \cite{CHR,thermalHallpaper}. 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 $\sigma^A$ are absolute, while $\nu^A,\eta^A,$ and $\kappa^A$ are relative. The combination $\nu^A-\dfrac{d\eta^A}{dT}$ is also absolute.

This distinction has consequences for the microscopic formulas that can be derived for transport coefficients which are usually called Kubo formulas \cite{thermalHallpaper}. 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
\begin{align}\label{eq: div jn}
\frac{\partial\rho^N}{\partial t}=-\nabla\cdot {\bf j}^N &=-\sigma_{km}^S \partial_k {\bf E}_m-\frac{{d\sigma}_{km}} {dT} {\bf E}_m \partial_k T-\nu^S_{km}\partial_k\partial_m T-\frac{{d\nu}^S_{km}}{dT} \partial_k T\partial_m T,\\
\begin{split}\label{eq: div je}
\frac{\partial\rho^E}{\partial t}=-\nabla\cdot {\bf j}^E &=-\varphi \nabla\cdot {\bf j}^N+\sigma^S_{km}{\bf E}_k{\bf E}_m-\left(\nu_{km}+\frac{d{\eta}_{mk}}{dT}\right) {\bf E}_k \partial_m T-\eta^S_{km} \partial_k{\bf E}_m\\
&+\kappa^S_{km}\partial_k\partial_m T+\frac{d{\kappa}^S_{km}}{dT}\partial_k T\partial_m T,
\end{split}
\end{align}
where the external fields are assumed to be time-independent and thus  $\nabla\times {\bf E}=0$. As expected, these time derivatives are unaffected by the transformations (\ref{eq: antisym part shifts}).

The time-derivative of entropy density is usually written as follows \cite{LandauLifshits}:
\begin{equation}\label{eq:entropy prod}
\frac{\partial s}{\partial t}=\frac{1}{T}\sigma^S_{km} {\bf E}_k {\bf  E}_m+\frac{1}{T^2}\kappa^S_{km} \partial_k T\partial_m T-\frac{1}{T^2}(T\nu_{km}+ \eta_{mk}) {\bf E}_k \partial_m T-\nabla\cdot {\tilde {\bf j}}^S_k,
\end{equation}
where the entropy current density is
\begin{equation}\label{entropycurrent}
{\tilde {\bf j}}^S_k=\frac{1}{T}\eta_{km} {\bf E}_m-\frac{1}{T}\kappa^S_{km} \partial_m T .
\end{equation}
The r.h.s. of eq. (\ref{eq:entropy prod}) seems to depend on some relative transport coefficients. However, if one redefines the entropy current as follows:
\begin{equation}
{\bf j}^S_k=\frac{1}{T}\eta^S_{km}{\bf E}_m-\frac{1}{T}\kappa^S_{km}\partial_m T,
\end{equation}
then one can write
\begin{equation}
\frac{\partial s}{\partial t}=\frac{1}{T}\sigma^S_{km} {\bf E}_k {\bf  E}_m+\frac{1}{T^2}\kappa^S_{km} \partial_k T\partial_m T-\frac{1}{T^2}\left(T\nu^S_{km}+ \eta^S_{mk}+\epsilon_{km}T\left(\nu^A-\frac{d{\eta}^A}{dT}\right)\right) {\bf E}_k \partial_m T-\nabla\cdot {\bf j}^S_k.
\end{equation}
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 $\varphi$ 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 $\nu^A$ drops out because
\begin{equation}
\int_{-\infty}^\infty dy\,  \nu^A \partial_y T dy=0.
\end{equation}

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 $\nu^A$ interpolates between $\nu^A_1$ for $y\ll 0$ and $\nu^A_2$ for $y\gg 0$. Suppose that the temperature is a function of $y$ only which is equal to $T_b$ throughout the interface region and approaches $T_\infty$ at $y\rightarrow\pm\infty$. Then the $\nu^A$-dependent contribution to the net electric current in the $x$ direction is
\begin{align}
 -\int_{-\infty}^\infty dy\, \nu^A \partial_y T  &  = (T_b-T_\infty) (\nu_2^A - \nu_1^A),
\end{align}
Similarly, the contribution of $\kappa^A$ to the net heat current is
$(T_b-T_\infty) (\kappa_2^A - \kappa_1^A).$
By creating an electric potential $\varphi$ which is equal to $\varphi_b$ in the interface region and approaches $\varphi_\infty$ at $y\rightarrow\pm\infty$, one can also measure $\sigma^A_2-\sigma^A_1$ and $\eta_2^A-\eta^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 $\sigma^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 $\eta^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. \ref{fig:flux}). Therefore heat transport will receive  a contribution from the work of the electromotive force $\mathcal E$ on the net electric edge currents $I^N_{\rm edge}$. The edge currents are proportional to the jump in the magnetization along the boundary and they make $\eta^A$ relative even in the  cylinder geometry.

\begin{figure}
    \centering
    \includegraphics{Fig1.pdf}
    \caption{Insertion of a flux $\Phi$ into a cylinder creates electric field around it. If $\eta_{xy}\neq 0$, this electric field $\bf E$ drives a heat current $I^{\mathcal Q}$ along the axis of the cylinder. The jump of the magnetization gives rise to an edge current $I^N_{\rm edge}$. Work done on the edge current by the electric field contributes to the net heat $\mathcal Q$ transferred to the heat bath.
    }
    \label{fig:flux}
\end{figure}

\section{Microscopic formulas for thermoelectric coefficients}\label{sec: themoelectric Kubo}
\subsection{Currents on a lattice}

We follow the conventions of \cite{thermalHallpaper}. We consider a lattice system with a Hamiltonian $H=\sum_{p\in\Lambda} H_p$, where $\Lambda\subset{\mathbb 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=\sum_{p\in\Lambda} 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\in\Lambda$.

The electric current from site $q$ to site $p$ is defined as
\begin{equation}
J^N_{pq}=i[H_q,Q_p]-i[H_p,Q_q].
\end{equation}
The energy current from site $q$ to site $p$ is defined as
\begin{equation}\label{energycurrent}
J^E_{pq}=-i[H_p,H_q].
\end{equation}
These currents enter the charge and energy conservation equations
\begin{align}
    \frac{d Q_q}{dt} &= i\sum_{p\in \Lambda} [H_p,Q_q] = -\sum_{p\in \Lambda} J^N_{pq},\\
    \frac{d H_q}{dt} &= i\sum_{p\in \Lambda} [H_p,H_q] = -\sum_{p\in \Lambda} J^E_{pq}.
\end{align}

The net current from a subset $B\subset \Lambda$ to its complement $A=\Lambda \backslash B$ is given by
\begin{align}
    J^{N,E}(A,B) = \sum_{p\in A}\sum_{q\in B}  J^{N,E}_{pq}.
\end{align}
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  $\eta(p,q): \Lambda \times \Lambda \rightarrow \mathbb R$ satisfying
\begin{equation}
 \eta(p,q)+\eta(q,r) + \eta(r,p) = 0, \quad \forall p,q,r \in \Lambda,
\end{equation}
define
\begin{align}
    J^{N,E}(\eta) = \frac 1 2 \sum_{p,q\in \Lambda} \eta(p,q)J^{N,E}_{pq}.
\end{align}
Current $J(A,B)$ from $B$ to $A=\Lambda\backslash B$ corresponds to the case $\eta(p,q)= \chi_B(q)-\chi_B(p)$, where the function $\chi_B(p)$ is equal 1 on $B$ and 0 otherwise.

The equilibrium expectation value of the currents satisfy
\begin{equation} \label{eq: div j = 0 lat}
\sum_{p\in \Lambda}\langle J^{N,E}_{pq} \rangle=0.
\end{equation}
This equation is a lattice analog of the continuum equation
\begin{align}
    \nabla \cdot \langle {\bf J}^{N,E}({\bf r})\rangle = 0.
\end{align}
In the continuum the general solution to this equation
\begin{align}\label{eq: cont magnetization}
    \langle J_k^{N,E}({\bf r})\rangle = -\epsilon_{kj}\partial_j M^{N,E}({\bf r})
\end{align}
defines the magnetization $M^N$ and energy magnetization $M^E$. Analogously, on a lattice the solution of equation (\ref{eq: div j = 0 lat}) is
\begin{equation}\label{magnetization}
\langle J^{N,E}_{qr}\rangle=\sum_{p\in \Lambda} M^{N,E}_{pqr},
\end{equation}
where $M^{N,E}_{pqr}$ are skew-symmetric functions of the lattice points $p,q,r\in \Lambda$. These are lattice analogs of the magnetization and the energy magnetisation. 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
\begin{equation}
M^{N,E}_{pqr}\mapsto M^{N,E}_{pqr}+\sum_{s\in \Lambda} N_{pqrs},
\end{equation}
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 $\sum_{r\in \Lambda}M^{N,E}_{pqr}=0$ which are not of the form $\sum_{s\in \Lambda} 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 (\ref{magnetization}) with respect to a parameter $\lambda^\ell$ of the Hamiltonian
\begin{equation}\label{magnetizationderivative}
\frac{\partial}{\partial \lambda^\ell}\langle J^{N,E}_{pq}\rangle=\sum_{r\in\Lambda} \mu^{N,E}_{pqr,
\ell},
\end{equation}
    where $\mu^{N,E}_{pqr,\ell}=\dfrac{\partial M^{N,E}_{pqr}}{\partial \lambda^\ell}$ and is given by \cite{kitaev}

\begin{equation}\label{mue}
\mu^{N,E}_{pqr,\ell}=-\beta\langle\langle \frac{\partial H_p}{\partial \lambda^\ell}; J^{N,E}_{qr}\rangle\rangle-\beta\langle\langle \frac{\partial H_r}{\partial \lambda^\ell}; J^{N,E}_{pq}\rangle\rangle-\beta\langle\langle \frac{\partial H_q}{\partial \lambda^\ell}; J^{N,E}_{rp}\rangle\rangle,
\end{equation}
where $\langle\langle A;B\rangle\rangle$ denotes the Kubo canonical pairing \cite{Kubo}. Using the properties of the Kubo pairing (see Appendix A),
one can easily verify the identity (\ref{magnetizationderivative}). In the following we will combine derivatives of magnetizations with respect to different parameters into 1-forms on the parameter space $\mu_{pqr}^{N,E}=\sum_{\ell} \mu^{N,E}_{pqr,\ell} d\lambda^{\ell}$.

\subsection{Equilibrium conditions and driving forces}
In the following sections we will follow Luttinger \cite{Luttinger} and study the behavior of the system coupled to external potentials
\begin{equation}
    H_p^{\psi,\varphi} =(1+\psi(p)) (H_p + \varphi(p) Q_p),
\end{equation}
where $\varphi(p)$ is external electric potential and $\psi(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
\begin{equation}
    \rho \sim \exp\left( -\frac{H^{\psi,\varphi}-\mu_0 Q}{T_0}\right),
\end{equation}
where $T_0$ and $\mu_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
\begin{equation}
    \rho(p) \sim \left( -\frac{H-\mu(p) Q}{T(p)}\right),
\end{equation}
where the local temperature $T(p)$ and chemical potential $\mu(p)$ are slowly varying functions of $p$. The equilibrium conditions can be found to be \cite{Luttinger,CHR,Niuetal}
\begin{align} \label{eq: elec eq}
    (1+\psi(p))(\mu(p)+\varphi(p)) &= \mu_0,\\ \label{eq: therm eq}
    (1+\psi(p)) T(p) &= T_0.
\end{align}
These relations together with the absence of transport currents in equilibrium can be used to derive Einstein relations between transport coefficients \cite{Luttinger}. The latter also leads to  transport currents being proportional to the gradients of the left hand sides of eqs. (\ref{eq: elec eq}) and (\ref{eq: therm eq}). The fact that the driving forces depend only on specific combinations of $\psi,\varphi,T, \mu $ will be used to relate the response to variations of the thermodynamic parameter $T$  to the  response to variations of the external field $\psi$.

Note that, in section \ref{sec: hydro} for simplicity we assumed that chemical potential is constant. We will continue doing so in the rest of the paper and we will assume $\mu(p)$ to be constant. Generalization to a more general case of non-constant local chemical potential is straightforward, but clumsy due to ambiguities in the separation of electrochemical potential $\mu+\varphi$ into two parts.

\subsection{Nernst effect}
In order to find the thermoelectric coefficient coefficient $\nu_{xy}$ we deform the Hamiltonian density by
\begin{align}\label{eq: thermal pert}
    \Delta H_p = \epsilon e^{st} {\overset{\frown}{g}}(p) H_p,
\end{align}
where ${\overset{\frown}{g}}(p)$ is a hat-shaped function as in Fig. \ref{fig: g hat fig 2d}, $\epsilon$ 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 \cite{Luttinger} in which the characteristic time $1/s$ is large but not large enough in order for the two slopes of the hat ${\overset{\frown}{g}}(p)$ to come into equilibrium.

\begin{figure}[ht]

     \adjustbox{minipage=1.3em,valign=t}{\subcaption{}\label{fig: g kink fig 2d}}
      \begin{subfigure}[b]{0.45\textwidth}
    \centering
    \includegraphics[valign=t]{Fig2a.pdf}
\end{subfigure}
 \adjustbox{minipage=1.3em,valign=t}{\subcaption{}\label{fig: g hat fig 2d}}
  \begin{subfigure}[b]{0.45\textwidth}
    \centering
    \includegraphics[valign=t]{Fig2b.pdf}
\end{subfigure}
 \adjustbox{minipage=1.3em,valign=t}{\subcaption{}\label{fig: g kink fig}}
  \begin{subfigure}[b]{0.3\textwidth}
    \centering
    \includegraphics[valign=t]{Fig2c.pdf}
\end{subfigure}
\hfill
 \adjustbox{minipage=1.3em,valign=t}{\subcaption{}\label{fig: g hat fig}}
\begin{subfigure}[b]{0.3\textwidth}
  \centering
  \includegraphics[valign=t]{Fig2d.pdf}
\end{subfigure}
\hspace*{5cm}
    \caption{ (a) Heat map of the function $g(p)$ corresponding to a temperature gradient in a horizontal strip. (b) Heat map of the function ${\overset{\frown}{g}}(p)$ corresponding to two horizontal strips with the opposite signs of the temperature gradient. (c) The function $g(p)$ restricted to the dashed line in (a). (d) The function ${\overset{\frown}{g}}(p)$ restricted to the dashed line in (b). Dotted line in (d) represents the dependence of some parameter $\lambda$ of the Hamiltonian on $y(p)$. }
    \label{fig:g}
\end{figure}

The change of the state of the system can be found as follows. The density matrix
\begin{align}
    \rho(t) = \rho_0 + \Delta \rho (t)
\end{align}
 satisfies the quantum Liouville equation
 \begin{equation}
    \frac{d \Delta \rho }{dt}= - i [\Delta H,\rho_0] + \dots,
 \end{equation}
where $\rho_0$ is the equilibrium density matrix at $t=-\infty$ and dots represent higher order terms in $\epsilon$. The solution to this equation is
\begin{align}
    \Delta \rho (t=0) = - \rho_0 \int_0^\infty dt \int_0^\beta d\tau \Delta \dot H (-t-i\tau) ,
\end{align}
where the dot denotes the time derivative.
The change of the observable $A$ can be found to be
\begin{align}
    \Delta \langle A \rangle = \langle \Delta A\rangle -\beta \int_0^\infty dt   \langle \langle A;  \Delta \dot H (-t) \rangle \rangle,
\end{align}
where we used the Kubo pairing notation (see Appendix \ref{app: Kubo}) and $\Delta A$ is the  variation of the operator arising from explicit dependence of $A$ on the Hamiltonian.
Using explicit form of the perturbation (\ref{eq: thermal pert}), energy conservation law and properties of the Kubo pairing we can rewrite this formula as
\begin{align}
    \Delta \langle A \rangle = \langle \Delta A\rangle +\epsilon\beta  \int_0^\infty dt e^{-st} \langle \langle A(t);  J^E(\delta {\overset{\frown}{g}}) \rangle \rangle,
\end{align}
where we neglected term proportional to small $s$.

The change in the electric current across a vertical line $x=a$ is
\begin{equation}
       \Delta \langle J^N(\delta f) \rangle = \langle \Delta J^N(\delta f) \rangle +\epsilon\beta  \int_0^\infty dt e^{-st} \langle \langle J^N(\delta f,t) ;  J^E(\delta {\overset{\frown}{g}}) \rangle \rangle,
\end{equation}
where $f(p)=\theta(a-x(p))$ is a step function. The explicit variation of the current is
\begin{align}
    \begin{split}
    &\langle \Delta J^N(\delta f) \rangle= \frac {i \epsilon} 2 \sum_{p,q\in \Lambda}\big\langle {\overset{\frown}{g}}(q)[H_q,Q_p] - {\overset{\frown}{g}}(p)[H_p,Q_q]\big\rangle (f(q)-f(p))\\&=\frac \epsilon 4 \sum_{p,q\in \Lambda} \langle J^N_{pq}\rangle({\overset{\frown}{g}}(p)+{\overset{\frown}{g}}(q)) (f(q)-f(p))\\&\qquad + \frac {i\epsilon} 4 \sum_{p,q\in \Lambda} \langle [H_q,Q_p]+[H_p,Q_q](f(q)-f(p))({\overset{\frown}{g}}(q)-{\overset{\frown}{g}}(p))\rangle,
    \end{split}
\end{align}
where in the last line we separated the result into two contribution formally skew-symmetric and symmetric in $f,{\overset{\frown}{g}}$. The second term depends only on difference of values of ${\overset{\frown}{g}}$ at different sites as expected for a   transport current. On the other hand, the first term depends on the value of ${\overset{\frown}{g}}$ and does not seem to be of the form expected for a transport current. As was explained in \cite{CHR} this contribution is related to magnetization currents. Indeed we can rewrite
\begin{align}\label{eq: j to mag}
    \frac \epsilon 4 \sum_{p,q\in \Lambda} \langle J^N_{pq}\rangle({\overset{\frown}{g}}(p)+{\overset{\frown}{g}}(q)) (f(q)-f(p)) = \sum_{p,q,r\in \Lambda} M^{N}_{pqr} (f(q)-f(p))({\overset{\frown}{g}}(q)-{\overset{\frown}{g}}(r)),
\end{align}
where we used skew-symmetry of $M_{pqr}^N$.

Written in this form, the response is proportional to the differences of ${\overset{\frown}{g}}$ at different points and therefore receives appreciable contributions only from regions $I$ and $II$ in Fig.~\ref{fig: g hat fig}. However, the transformation (\ref{eq: j to mag}) 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 \cite{CHR} 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  \cite{thermalHallpaper}. Instead of introducing a sharp boundary, we will make one parameter of the Hamiltonian $\lambda$ to have slightly different values in regions $I$ and $II$ (see Fig.~\ref{fig: g hat fig}). We can write the hat-shaped function ${\overset{\frown}{g}}$ as a difference of two functions $g_I$ and $g_{II}$, ${\overset{\frown}{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.~\ref{fig: g kink fig}. The functions $g_{I,II}$  are non-constant only in regions $I$ and $II$ respectively. Then the magnetization contribution can be rewritten as
\begin{equation}
     \epsilon (\lambda_{II}-\lambda_{I}) \mu^N_\lambda(\delta g\cup \delta f) + O((\lambda_{II}-\lambda_{I}) ^2),
\end{equation}
where $\mu^N_{pqr,\lambda} = \frac{\partial M^N_{pqr}}{\partial \lambda}$, and we introduced a notation
\begin{align}\label{eq: mag contraction}
    \mu^{N}_\lambda(\delta f_1 \cup \delta f_2) = \frac 1 6 \sum_{p,q,r\in\Lambda} \mu^N_{pqr,\lambda} (f_1(q)-f_2(p))(f_2(r)-f_2(q)).
\end{align}
Combining this with other contributions we find
\begin{align}\label{eq: nerst current}
\begin{split}
    \delta \langle J^N(\delta f)\rangle &\approx \epsilon (\lambda_{II}-\lambda_{I}) \Bigg\{\frac{\partial} {\partial \lambda}\Big[\beta  \int_0^\infty dt e^{-st} \langle \langle J^N(\delta f,t) ;  J^E(\delta g) \rangle \rangle +U(\delta f,\delta g)\Big]+ \mu^N_\lambda(\delta g\cup \delta f)\Bigg\},
\end{split}
\end{align}
where
\begin{align}\label{eq: U def}
    U(\delta f,\delta g) = \frac {i} 4 \sum_{p,q\in \Lambda} \langle [H_q,Q_p]+[H_p,Q_q](f(q)-f(p))(g(q)-g(p))\rangle.
\end{align}

The function $g$ as in Fig. \ref{fig: g kink fig} is not compactly supported and thus it takes an  infinite time for the system to equilibrate. Therefore, one can take the limit $s\rightarrow0$ while staying in the ``fast'' regime.
From eqs. (\ref{eq: elec eq},\ref{eq: therm eq}) one finds that electric current after perturbation by gravitational potential $\epsilon g(p)$ is equal to the current generated by
\begin{align}
    T(p) &= \epsilon g(p) T_0,\\
    \varphi(p) &= \epsilon g(p) \mu_0.
\end{align}
From continuum phenomenological equation (\ref{eq: jn hydro}) one finds the current across $x=a$ line to be
\begin{equation}
    -\epsilon T_0 \int_{-\infty}^\infty \nu_{xy}\partial_y g dy  -\epsilon \mu_0 \int_{-\infty}^\infty \sigma_{xy} \partial_y g dy = \epsilon T_0 \int_{\lambda_I}^{\lambda_{II}} \frac{\partial\nu_{xy}}{\partial \lambda}d\lambda  + \epsilon \mu_0 \int_{\lambda_I}^{\lambda_{II}} \frac{\partial\sigma_{xy}}{\partial \lambda}d\lambda .
\end{equation}
Comparing it to (\ref{eq: nerst current}) we find
\begin{align}
\begin{split}
    d \nu_{xy} &=d\Big[\beta^2 \lim_{s\rightarrow0} \int_0^\infty dt e^{-st} \langle \langle J^N(\delta f,t) ;  J^{\mathcal Q}(\delta g) \rangle \rangle +\beta U(\delta f,\delta g) \Big]- \beta \mu^N(\delta f \cup\delta g ),
\end{split}
\end{align}
where we introduced the notation $J^{\mathcal Q}=J^E-\mu J^N$ for the heat current, inverse temperature $\beta=\frac 1 T$, we dropped the subscript 0 from $T_0$ and $\mu_0$ since this formula contains correlation functions of the unperturbed system in equilibrium. We combined differential with respect to parameter into the differential form $\mu^N(\delta f \cup\delta g ) = \mu^N_\lambda(\delta f \cup\delta g )d\lambda$ and the derivation can be straightforwardly extended to involve several parameters. The exterior derivative $d= \sum_\ell d\lambda \frac {\partial}{\partial \lambda^\ell}$ 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 $\eta_{xy}$ for any two 2d materials, regardless of the temperature. Explicitly, let us define the rescaled Hamiltonian $H_{\lambda_0}=\lambda_0 H$, where we introduced an additional scaling parameter $\lambda_0$. Then
\begin{equation}
\left(T\frac{\partial}{\partial T}+\left.{\lambda_0}\frac{\partial}{\partial{\lambda_0}}\right|_{{\lambda_0}=1}\right)\nu_{xy}({\lambda_0},T)=0.
\end{equation}
Therefore we can define the $T$-component of the 1-form on the enlarged parameter space as follows:
\begin{equation} \label{eta Tderivative}
\frac{d\nu_{xy}}{dT}=\frac{\partial}{\partial T}\left[\beta^2  \int_0^\infty dt e^{-st} \langle \langle J^N(\delta f,t) ;  J^{\mathcal Q}(\delta g) \rangle \rangle +\beta U(\delta f,\delta g) \right]-\beta^2\tau^N(\delta f \cup \delta g),
\end{equation}
where $\tau^N(\delta f \cup \delta g)$ is given by eq. (\ref{eq: mag contraction}) with $\mu_{pqr}$ replaced with
\begin{equation}\label{eq: tau def}
\tau^N_{pqr}=\beta\langle\langle H_p; J^N_{qr}\rangle\rangle+\beta\langle\langle H_r; J^N_{pq}\rangle\rangle+\beta\langle\langle H_q; J^N_{rp}\rangle\rangle,
\end{equation}
which is obtained from $\mu^N$ by replacing $dH_p$ with $-H_p$.

\subsection{Ettingshausen effect}

One can derive a formula for the coefficient $\eta_{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
\begin{align}
    \Delta H_p = \epsilon e^{st} {\overset{\frown}{g}}(p) Q_p,
\end{align}
where ${\overset{\frown}{g}}(p)$ is a hat-shaped function of $y(p)$ as in Fig.~\ref{fig: g hat fig 2d}.

The change in the energy current across a vertical line $x=a$ is
\begin{equation}
       \Delta \langle J^E(\delta f) \rangle = \langle \Delta J^E(\delta f) \rangle +\epsilon\beta  \int_0^\infty dt e^{-st} \langle \langle J^E(\delta f,t) ;  J^N(\delta {\overset{\frown}{g}}) \rangle \rangle,
\end{equation}
where $f=\theta(a-x(p))$. The explicit variation of the current is
\begin{align}
    \begin{split}
    \langle \Delta J^E(\delta f) \rangle&= \frac {i \epsilon} 2 \sum_{p,q\in \Lambda}\big\langle {\overset{\frown}{g}}(p)[H_q,Q_p] - {\overset{\frown}{g}}(q)[H_p,Q_q]\big\rangle (f(q)-f(p))\\&\qquad=\frac \epsilon 4 \sum_{p,q\in \Lambda} \langle J^N_{pq}\rangle({\overset{\frown}{g}}(p)+{\overset{\frown}{g}}(q)) (f(q)-f(p)) - U(\delta f, \delta {\overset{\frown}{g}}),
    \end{split}
\end{align}

The first term can be expressed in terms of magnetization as in (\ref{eq: j to mag}). We write ${\overset{\frown}{g}}(p)$ as a difference  ${\overset{\frown}{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:

\begin{align}\label{eq: etting current}
\begin{split}
    \delta \langle J^E(\delta f)\rangle &\approx \epsilon (\lambda_{II}-\lambda_{I}) \Bigg\{\frac{\partial} {\partial \lambda}\Big[\beta  \int_0^\infty dt e^{-st} \langle \langle J^E(\delta f,t) ;  J^N(\delta g) \rangle \rangle -U(\delta f,\delta g)\Big]+ \mu^N_\lambda(\delta g\cup \delta f)\Bigg\}.
\end{split}
\end{align}

On the other hand, from the continuum phenomenological equation (\ref{eq: je hydro}) one finds the current across the line $x=a$ to be
\begin{equation}
     -\epsilon  \int_{-\infty}^\infty \eta_{xy} \partial_y g dy  -\epsilon \mu \int_{-\infty}^\infty \sigma_{xy} \partial_y g dy =  \epsilon \int_{\lambda_I}^{\lambda_{II}} \frac{\partial\eta_{xy}}{\partial \lambda}d\lambda +\epsilon\mu \int_{\lambda_I}^{\lambda_{II}} \frac{\partial\sigma_{xy}}{\partial \lambda}d\lambda,
\end{equation}
where the second term originates from the contribution of $\varphi {\bf j}^N$ to ${\bf j}^E$. Comparing it to (\ref{eq: etting current}) we find
\begin{align}
\begin{split}
    d \eta_{xy} &=d\Big[\beta  \lim_{s\rightarrow0}\int_0^\infty dt e^{-st} \langle \langle J^{\mathcal Q}(\delta f,t) ;  J^N(\delta g) \rangle \rangle - U(\delta f,\delta g) \Big]- \beta \mu^N(\delta f \cup\delta g ).
\end{split}
\end{align}

This 1-form can be extended to include temperature as a parameter using the scaling relation
\begin{equation}
\left(T\frac{\partial}{\partial T}+\left.\lambda_0\frac{\partial}{\partial\lambda_0}\right|_{\lambda_0=1}\right) \frac {\eta_{xy}(\lambda_0,T)}{T}=0.
\end{equation}
Therefore we can define the $T$-component of the 1-form on the enlarged parameter space as follows:
\begin{equation}
\frac{d}{dT} \frac {\eta_{xy}}{T}=\frac{\partial}{\partial T}\left[\beta^2  \int_0^\infty dt e^{-st} \langle \langle J^{\mathcal Q}(\delta f,t) ;  J^N(\delta g) \rangle \rangle -\beta U(\delta f,\delta g) \right]-\beta^2\tau^N(\delta f \cup \delta g),
\end{equation}
where $\tau^N$ is given by (\ref{eq: tau def}).

\subsection{Symmetric parts of transport coefficients}

Note that the 1-forms $\mu^N(\delta f\cup\delta g)$ and $\tau^N(\delta f\cup\delta 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 (NESS). Replacing $f$ with a smeared step-function changes the operators  $J^N(\delta f)$ and $J^E(\delta f)$ by $J^N(\delta {\overset{\frown}{f}})$ and  $J^E(\delta {\overset{\frown}{f}})$. The latter  operators can be equivalently written as $\frac{d Q({\overset{\frown}{f}})}{dt}$ and $\frac{d H({\overset{\frown}{f}})}{dt}$, where
\begin{equation}
Q({\overset{\frown}{f}})=\sum_{p\in \Lambda} {\overset{\frown}{f}}(p) Q_p,\quad H({\overset{\frown}{f}})=\sum_{p\in \Lambda} {\overset{\frown}{f}}(p) H_p,
\end{equation}
and ${\overset{\frown}{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 $\nu_{xy}$ or $\eta_{xy}$.

After this has been done, exchanging $f$ and $g$ is equivalent to exchanging $x$ and $y$. Thus $\mu^N(\delta f\cup\delta g)$ does not enter the microscopic formulas for the symmetrized thermoelectric coefficients. These formulas then can be integrated, giving
\begin{align}  \nu^S_{xy} &=\frac{\beta^2}2  \lim_{s\rightarrow0} \int_0^\infty dt e^{-st} \left[\langle \langle J^N(\delta f,t) ;  J^{\mathcal Q}(\delta g) \rangle \rangle +\langle \langle  J^N(\delta g,t) ;  J^{\mathcal Q}(\delta f) \rangle \rangle \right]+\beta U(\delta f,\delta g),\\
\eta^S_{xy} &=\frac \beta 2  \lim_{s\rightarrow0}\int_0^\infty dt e^{-st} \left[\langle \langle J^{\mathcal Q}(\delta f,t) ;  J^N(\delta g) \rangle \rangle+\langle \langle J^{\mathcal Q}(\delta g,t) ;  J^N(\delta f) \rangle \rangle\right] - U(\delta f,\delta g).
\end{align}
As we show in Appendix \ref{App: invariance under Hamiltonian density redefinition}, 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(\delta f,\delta g)$ is zero for fermionic systems with only density-dependent interactions (see Appendix \ref{app: thermoelectric free fermions} for more details).

\subsection{Skew-symmetric parts of transport coefficients}

In a similar way one can find skew-symmetric parts of transport coefficients
\begin{align}
        d \nu^A &=d\Big[\frac {\beta^2}2 \lim_{s\rightarrow0} \int_0^\infty dt e^{-st} \left( \langle \langle J^N(\delta f,t) ;  J^{\mathcal Q}(\delta g) - \langle \langle J^N(\delta g,t) ;  J^{\mathcal Q}(\delta f\rangle \rangle\right) \Big]- \beta \mu^N(\delta f \cup\delta g ),\\
            d \eta^A &=d\Big[\frac \beta 2  \lim_{s\rightarrow0}\int_0^\infty dt e^{-st} \left(\langle \langle J^{\mathcal Q}(\delta f,t) ;  J^N(\delta g) \rangle \rangle-\langle \langle J^{\mathcal Q}(\delta g,t) ;  J^N(\delta f) \rangle \rangle \right) \Big]- \beta \mu^N(\delta f \cup\delta g ).
\end{align}
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\nu^A$ or $d\eta^A$ to depend only on the endpoints of the path. This means that the 1-form $\mu^N(\delta f \cup\delta g )$ must be exact. This can be proved using the techniques of \cite{thermalHallpaper} where   $\mu^E(\delta f \cup\delta g )$ was  shown to be exact.

\section{Discussion}\label{sec: thermoel 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 \cite{thermalHallpaper}.

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\rightarrow 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 $\theta$. 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\cdot\tan \theta $, then the correlator
\begin{equation}
\beta \lim_{s\rightarrow 0+}\int_0^\infty e^{-st} \langle\langle J^N(\delta f,t); J^N(\delta g)\rangle\rangle dt
\end{equation}
measures $\sigma_{xy}-\frac{1}{\tan\theta}\sigma_{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 ${\mathbb R}^2$ with a lattice in ${\mathbb R}^2\times [0,L]$, impose periodic boundary conditions in the third direction, divide all formulas by $L$, and take the limit $L\rightarrow\infty$. The functions $f,g$ remain independent of the third coordinate. It should not matter whether the limit $L\rightarrow\infty$ is taken before or after the limit $s\rightarrow 0$, since the problem is translationally invariant in the third direction.

In the 2d case, the quantity $\nu^A(T)$ (normalized relative to the vacuum) is dimensionless, and Bloch's theorem implies that for gapped systems $\nu^A(0)$ does not change under the variations of the Hamiltonian which do not close the gap. Thus if $\nu^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  $\nu^A(0)$ vanishes for all gapped systems  \cite{Thirdlaw}. By Onsager reciprocity, the $T\rightarrow 0$ limit of $\eta^A(T)/T$ also vanishes. Thus topological invariants of gapped 2d systems arise only from the Hall conductivity and the thermal Hall conductivity.