Hydrodynamic charge and heat transport on inhomogeneous curved spaces

TL;DR

This paper develops a hydrodynamic theory for charge and heat transport in strongly interacting systems on curved spaces, highlighting the impact of spatial inhomogeneity and viscous effects, with relevance to materials like suspended graphene.

Contribution

It introduces a novel hydrodynamic framework for inhomogeneous curved spaces, analyzing the effects of spatial curvature on transport properties in condensed matter systems.

Findings

Thermal and electrical conductivities are dominated by viscous effects.

Thermal conductivity is highly sensitive to spatial inhomogeneity.

The theory is applicable to systems like suspended graphene at the Dirac point.

Abstract

We develop the theory of hydrodynamic charge and heat transport in strongly interacting quasi-relativistic systems on manifolds with inhomogeneous spatial curvature. In solid-state physics, this is analogous to strain disorder in the underlying lattice. In the hydrodynamic limit, we find that the thermal and electrical conductivities are dominated by viscous effects, and that the thermal conductivity is most sensitive to this disorder. We compare the effects of inhomogeneity in the spatial metric to inhomogeneity in the chemical potential, and discuss the extent to which our hydrodynamic theory is relevant for experimentally realizable condensed matter systems, including suspended graphene at the Dirac point.

Figures

Equations

$$\mathrm{d}{s}^2=\left({\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}}_{ij}+{\partial }_{i}h{\partial }_{j}h\right)\mathrm{d}{x}_i\mathrm{d}{x}_j.$$

$${\nabla }_{\mu }{T}^{\mu \nu }=F_{\mathrm{ext}}^{\nu \mu }{J}_{\mu },$$

$${\nabla }_{\mu }{J}^{\mu }=0.$$

$$A_{\mathrm{ext}}={\mu }_0(\mathbf{x})\mathrm{d}t.$$

$$\mathrm{d}{s}^2=g_{\mu \nu }\mathrm{d}{x}^{\mu }\mathrm{d}{x}^{\nu }=-\mathrm{d}{t}^2+g_{ij}(\mathbf{x})\mathrm{d}{x}_i\mathrm{d}{x}_j,$$

$$T^{\mu \nu }$$

$$J^{\mu }$$

$$\left(\begin{array}[]{c}J_{i}^{\mathrm{avg}}\\ Q_{i}^{\mathrm{avg}}\end{array}\right)\equiv\left(\begin{array}[]{cc}\sigma_{ij}&\ T\alpha_{ij}\\ T\bar{\alpha}_{ij}&\ T\bar{\kappa}_{ij}\end{array}\right)\left(\begin{array}[]{c}E_{j}\\ \zeta_{j}\end{array}\right),$$

$$Q^i\equiv T^{ti}-\mu (\mathbf{x})J^i,$$

$${\kappa }_{ij}\equiv {\bar{\kappa }}_{ij}-T{\bar{\alpha }}_{ik}{\sigma }_{kl}^{-1}{\alpha }_{lj}.$$

$${\mu }_{\mathrm{eq}}(\mathbf{x})$$

$$T_{\mathrm{eq}}(\mathbf{x})$$

$$u_{\mathrm{eq}}^{\mu }(\mathbf{x})$$

$$\mu (\mathbf{x})$$

$$T(\mathbf{x})$$

$$u^{\mu }(\mathbf{x})$$

$$J^t$$

$$J^j$$

$$T^{tt}$$

$$T^{ti}$$

$$T^{ij}$$

$$A$$

$$\mathrm{d}{s}^2=-\mathrm{d}{t}^2+g_{ij}(\mathbf{x})\mathrm{d}{x}^i\mathrm{d}{x}^j+2\frac{{\mathrm{e}}^{-\mathrm{i}\omega t}}{\mathrm{i}\omega }{\zeta }_i\mathrm{d}{x}^i\mathrm{d}t.$$

$$-{\bar{\nabla }}^i\left({\sigma }_{\text{textsc}q}{\partial }_i\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}\mu \right)+{\bar{\nabla }}^i\left({\sigma }_{\text{textsc}q}\frac{{\mu }_0}{T_0}{\partial }_i\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}T\right)+{\bar{\nabla }}_i\left(n\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}{v}^i\right)$$

$${\bar{\nabla }}^i\left({\sigma }_{\text{textsc}q}{\mu }_0{\partial }_i\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}\mu \right)-{\bar{\nabla }}^i\left({\sigma }_{\text{textsc}q}\frac{{\mu }_0^2}{T_0}{\partial }_i\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}T\right)+{\bar{\nabla }}_i\left(s{T}_0\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}{v}^i\right)$$

$$n{\partial }_j\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}\mu +s{\partial }_j\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}T-{\bar{\nabla }}_i\left[\eta ({\bar{\nabla }}^i\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}{v}_j+{\bar{\nabla }}_j\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}{v}^i)\right]-{\partial }_j\left[(\zeta -\eta ){\bar{\nabla }}_i\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}{v}^i\right]$$

$${\sigma }^{ij}\approx \frac{n^2{\tau }^{ij}}{ϵ+P},$$

$${\alpha }^{ij}\approx \frac{ns{\tau }^{ij}}{ϵ+P},$$

$${\bar{\kappa }}^{ij}\approx \frac{T{s}^2{\tau }^{ij}}{ϵ+P},$$

$${\tau }_{ij}^{-1}=({\tau }_{ij}^{-1}{)}^{(\mu \mu )}+({\tau }_{ij}^{-1}{)}^{(\mu h)}+({\tau }_{ij}^{-1}{)}^{(hh)}$$

$$({\tau }_{ij}^{-1}{)}^{(\mu \mu )}=$$

$$({\tau }_{ij}^{-1}{)}^{(\mu h)}=$$

$$({\tau }_{ij}^{-1}{)}^{(hh)}=$$

$${\mathbb{P}}_{ij}={\text{mbox1 0 -.18 1to0.0pthssleavevmodeset@colorδ}}_{ij}-\frac{k_i{k}_j}{k^2}.$$

$$n({\mu }_0)=C_2\mu T+\mathrm{O}\left(\frac{{\mu }^3}{T}\right),$$

$$s({\mu }_0)=C_0{T}^2+\frac{C_2}{2}{\mu }^2+\mathrm{O}\left(\frac{{\mu }^4}{T^2}\right),$$

$$\eta ({\mu }_0)=T^2{\eta }_0+\mathrm{O}\left({\mu }^2\right),$$

$$\zeta ({\mu }_0)=0,$$

$$\sigma ({\mu }_0)={\sigma }_0+\mathrm{O}\left(\frac{{\mu }^2}{T^2}\right),$$

$${ϵ}_{\text{textsc}c}^2$$

Full text

1. 1 Introduction

2. 2 Relativistic Hydrodynamics on Curved Space

3. 3 Hydrodynamic Transport

4. 3.1 General Solution

5. 3.2 Perturbative Analytic Solution

6. 3.3 Numerical Solution

7. 4 Application to Suspended Graphene

8. 4.1 Classical Flexural Phonon Dynamics

9. 4.2 Contribution to Momentum Relaxation Time

10. 5 Conclusion

11. A Onsager Reciprocity

12. B Numerical Simulations

13. C Weak Disorder

14. C.1 Phonon Contribution

[TABLE]

**

Introduction**

A theory of electrical and thermal transport necessarily relies on a precise description of how translation symmetry is broken. In conventional weakly coupled quasi-particle theories, most collisions of electrons are with impurities or phonons, and relax momentum. In recent years, rapid progress towards a theory of transport which also accounts for momentum-conserving electron-electron interactions has been made [1]. One of the most useful tools that has arisen for understanding transport in this limit is hydrodynamics. Hydrodynamics is the effective theory describing the relaxation of any interacting system to thermal equilibrium on long wavelengths. Such a theory is suitable for any interacting metal where the disorder which breaks translation invariance only varies on long wavelengths compared to the electron-electron scattering length [2, 3, 4, 5, 6, 7, 8]. Although this is a difficult regime to reach experimentally, this has now become possible [9, 10, 11, 12] (see also [13]). A thorough understanding of the hydrodynamic regime of transport is certainly necessary as a “solvable” limit of any more complete theory of transport [14]. Hence, it is worthwhile to have a systematic understanding of hydrodynamic transport in a broad variety of systems.

The purpose of this paper is to describe hydrodynamic transport on curved spaces.111Our formalism is relatively similar to the emergent “hydrodynamic” formalism used to describe transport in strongly correlated systems described via the AdS/CMT correspondence [4, 15, 16, 17, 18]. However, in most of these papers, the random spatial metric is an emergent phenomenon from the point of view of the bulk description of the field theory — the exception is [19]. We emphasize that we are interested in scenarios where the inhomogeneous spatial metric is a physical effect. In electronic materials, the presence of internal strain on a crystal lattice can be interpreted as an effective distortion to the induced spatial metric [20]. As the electronic charge-carrying degrees of freedom move in this inhomogeneous metric, our results will be relevant for strongly correlated systems in inhomogeneously strained crystals. Following [7], we will focus on the relativistic hydrodynamic equations as a model for transport in monolayer graphene in the hydrodynamic limit. The techniques which we develop straightforwardly generalize to other hydrodynamic models.

Recent experimental evidence [9] indicates that electrons behave hydrodynamically in charge-neutral graphene. Collectively they behave as a Dirac fluid: a plasma of thermally excited electron and holes which is likely to be strongly interacting at ‘reasonable’ temperatures $T\sim 100$ K [21, 22, 23]. Crucial to the observation of this Dirac fluid is the reduction and smoothing of “charge puddle” disorder, which corresponds to inhomogeneities in the local chemical potential. This was achieved by placing the graphene sheet in between layers of another material: boron nitride [24]. Another way to reduce charge puddle disorder in graphene is to “suspend” graphene, leaving it unattached to any substrate [25, 26]. For mechanical reasons, dealing with such suspended graphene can be challenging. The aspect we focus on here is that in principle a suspended sheet of graphene, as it consists of a single two dimensional ‘membrane’ of carbon atoms, is susceptible to out-of-plane flexural distortions. From the point of view of a two-dimensional effective theory for the Dirac fluid, flexural disorder can be interpreted as disorder in the spatial components of the spacetime metric. Letting the local height of the membrane be $h(x,y)$, the metric is [20]

[TABLE]

In reality, $h(x,y)$ need not be time-independent. However, such flexural motion is expected to be quite slow relative to electronic time scales, and we may approximate it as static disorder. Hence, a study of hydrodynamic electron transport in suspended graphene should naturally include flexural distortions to the metric.

The outline of this paper, and our main conclusions, are as follows:

• $\blacktriangleright$

In Section 2, we review the theory of linearized relativistic hydrodynamics on curved spaces, relevant for transport.

• $\blacktriangleright$

In Section 3 we use this curved space hydrodynamics to solve for thermoelectric transport coefficients for a fluid in a slowly varying chemical potential and spatial metric. When the inhomogeneity is small, we give analytic expressions for the thermal and charge conductivities as functions, expressed entirely in terms of the inhomogeneous chemical potential and metric, and thermodynamic and hydrodynamic coefficients. When the inhomogeneity cannot be treated analytically, we compute the transport coefficients numerically. Because transport is dissipative, the transport coefficients depend on hydrodynamic dissipation, via viscosity and a “quantum critical” conductivity. In the presence of inhomogeneous chemical potentials, both dissipative channels affect the conductivity significantly. However, for inhomogeneous strain viscous dissipation is far more relevant – in fact, perturbatively it is the only source of dissipation.

• $\blacktriangleright$

We discuss the application of our formalism to suspended graphene in Section 4. This discussion includes a justification of some of the statements in the introduction. Our hydrodynamic transport theory allows us to describe electronic scattering off of certain long-wavelength phonons non-perturbatively in the strength of electronic interactions. Although we will see that most phonons cannot be accounted for in this limit, our results may nonetheless be valuable for a more detailed study of electron-phonon coupling.

Technical results are found in appendices. We mostly work in units where $\hbar=k_{\mathrm{B}}=1$, and we also set the effective speed of light $v_{\mathrm{F}}=1$,222In materials such as graphene, the effective speed of light is set by the Fermi velocity $v_{\mathrm{F}}$. as well as the electron charge $e=1$. When we discuss the application of our formalism to suspended graphene, we will briefly restore these dimensionful quantities.

**

Relativistic Hydrodynamics on Curved Space**

In this section we review and generalize to curved spacetime the hydrodynamic framework developed in [7]. This framework describes the collective motion of the relativistic electronic plasma in a disordered metal, where the disorder is introduced via a spatially dependent chemical potential $\mu_{0}(\mathbf{x})$. When the chemical potential varies on a length scale larger than the electron mean-free path, a hydrodynamic description of transport is sensible: all other microscopic degrees of freedom have already reached local thermodynamic equilibrium. The only relevant degrees of freedom for transport are locally conserved quantities: energy, charge and momentum. All the spatial dependence of the parameters (such as local energy density $\epsilon$ or shear viscosity $\eta$) is encoded by the functional dependence of these quantities on the local $\mu_{0}(\mathbf{x})$: e.g. $\eta(\mathbf{x})=\eta(\mu_{0}(\mathbf{x}))$. Charge/chemical potential disorder is natural for many metals, including graphene [7]. For slowly varying disorder, this is also convenient because it is very naturally included within a hydrodynamic framework.

Another type of universal disorder that is natural to consider within a hydrodynamic framework is local inhomogeneity in the spacetime metric: as we described previously, this is a model for strain in the crystal lattice. This strain can also be natural in a broad variety of solids: occuring from either in-plane strain, or (in the case of suspended graphene) out-of-plane bending of the crystal lattice. In the limit where this strain is long wavelength, we can account for it by simply solving the hydrodynamic equations of motion, written in a coordinate-independent fashion, on a curved spacetime.

Let us note that strain can also open up a gap $\Delta$ in certain crystals, including graphene [27]. This will alter the effective microscopic dispersion relation and hence the equations of state. In the present work we have neglected this contribution, and our theory is not valid if the strain is so large that $\Delta\sim T$. For smaller strain, our theory remains valid, but there will be additional $\mathbf{x}$-dependence of the thermodynamic and hydrodynamic coefficients due to the local value of the gap. For simplicity we will not explicitly acocunt for this effect. Up to the opening of a gap, the effects of strain are universal in the hydrodynamic limit.

As we previously noted, the only quantities which are globally conserved (up to external sources) are charge, energy and momentum [28]. The natural degrees of freedom are the local number density $n(x)$, the energy density $\epsilon(x)$ and the momentum density $\Pi_{i}(x)$. A more convenient approach is to use their thermodynamic conjugates: the chemical potential $\mu(x)$, temperature $T(x)$, and velocity $v^{i}(x)$ respectively. These are the standard choice of hydrodynamic variables. In relativistic systems this velocity is commonly written covariantly as a four-velocity $u^{\mu}(x)=(1,v^{i})/\sqrt{1-v^{2}}$, constrained to equal $u^{\mu}u^{\nu}\eta_{\mu\nu}=-1$ with $\eta_{\mu\nu}$ the Minkowski metric: $\eta_{\mu\nu}=\mathrm{diag}(-1,1,\ldots,1)$.

The equations of motion are local conservation laws:

[TABLE]

In the absence of an external electric field or temperature gradient, there remains an external electromagnetic field due to an inhomogeneous chemical potential: $F^{\mu\nu}_{\mathrm{ext}}=\nabla^{\mu}A^{\nu}_{\mathrm{ext}}-\nabla^{\nu}A^{\mu}_{\mathrm{ext}}$, with

[TABLE]

The only non vanishing components of the Maxwell tensor are $F_{ti}^{\mathrm{ext}}=-F_{it}^{\mathrm{ext}}=\partial_{i}\mu_{0}(\mathbf{x})$. $T^{\mu\nu}(\mathbf{x},t)$ and $J^{\mu}(\mathbf{x},t)$ are the expectation values of the local relativistic stress-energy tensor and charge current, respectively. These conservation equations, understood in terms of the covariant derivative $\nabla_{\mu}$ with respect to the metric

[TABLE]

are valid in any (curved) spacetime, including those with an inhomogenous spatial curvature of interest to us.

In order for the equations (2) to be well-posed, we must express the expectation value of $T^{\mu\nu}$ and $J^{\mu}$ in terms of the hydrodynamic variables $\mu$, $T$ and $u^{\mu}$. We will expand $T^{\mu\nu}$ and $J^{\mu}$ in a gradient expansion in derivatives: more physically, the small parameter of the perturbative expansion is $\ell_{\mathrm{ee}}k$, with $k$ the wave number of our perturbation, and $\ell_{\mathrm{ee}}$ the electron-electron scattering length. In this paper, we will only include terms with zero or one derivatives of $\mathbf{x}$ and $t$. This expansion is well-known for a relativistic fluid [28, 29]:

[TABLE]

with $\eta$ and $\zeta$ the shear and bulk viscosity respectively, and $\sigma_{\textsc{q}}$ a microscopic dissipative coefficient. As emphasized in [29], $\sigma_{\textsc{q}}$ should be interpreted as the finite electrical conductivity of the charge neutral plasma (up to hydrodynamic long-time tails [1]), and for historical reasons it is sometimes called the “quantum critical” conductivity. Finally, $\mathcal{P}^{\mu\nu}$ is the projector orthogonal to the rest frame of the fluid, set by the velocity $u^{\mu}$: $\mathcal{P}^{\mu\nu}=g^{\mu\nu}+u^{\mu}u^{\nu}$.

**

Hydrodynamic Transport**

We now wish to compute the thermoelectric conductivity matrix of a fluid in such an inhomogeneous background. These coefficients are defined as follows:

[TABLE]

where $J_{i}^{\mathrm{avg}}$ is the spatial average of the charge current defined above, $Q_{i}^{\mathrm{avg}}$ is the spatial average of the heat current, defined as

[TABLE]

$E_{j}$ is an infinitesimal externally applied uniform electric field, and $\zeta_{j}$ is an infinitesimal “thermal drive” analogous to a homogeneous temperature gradient $-\partial_{j}\log T$. This more formal notation will prove useful for our purposes.

Our goal is to compute $\sigma_{ij}$, $\alpha_{ij}$, $\bar{\alpha}_{ij}$ and $\bar{\kappa}_{ij}$, using the hydrodynamic equations of motion. We will explicitly show how this is done. First, let us note a few formal results. Onsager reciprocity states that (with time-reversal symmetry) $\bar{\alpha}_{ij}=\alpha_{ji}$, and that $\sigma$ and $\bar{\kappa}$ are symmetric. In the hydrodynamic framework on a curved space, we prove this in Appendix A. Secondly, it is experimentally more common to measure a thermal conductivity defined by

[TABLE]

This can be interpreted as the ratio of the average heat current to a constant temperature gradient, subject to the constraint that no net charge current flows. We will show results for $\bar{\kappa}$ and for $\kappa$.

**

General Solution**

We now present the formal computation of the thermoelectric conductivity matrix. First, we note that in an inhomogeneous metric $g_{ij}(\mathbf{x})$ and chemical potential $\mu_{0}(\mathbf{x})$, there is an exact solution to the nonlinear equations of motion, encoding that the fluid is at rest in local thermal equilibrium:

[TABLE]

where $T_{0}$ is a constant. Then, because we are applying an infinitesimal electric field and thermal drive, we only look for the perturbations around equilibrium within linear response:

[TABLE]

Because the disorder explicitly picks out a preferred fluid rest frame, it is often helpful to decompose (2) into timelike and spacelike components. The hydrodynamic expansion of the electric current within linear response gives:

[TABLE]

while the stress-energy tensor reads

[TABLE]

where $\bar{\nabla}_{i}v^{j}\equiv\partial_{i}v^{j}+\Gamma^{j}_{~{}kl}v^{j}$ is the covariant derivative with respect to the spatial metric $g_{ij}$ and $\Gamma^{j}_{~{}kl}=\frac{1}{2}g^{jm}(\partial_{k}g_{ml}+\partial_{l}g_{mk}-\partial_{m}g_{kl})$ is the Christoffel symbol. For simplicity, we henceforth specialize to two spatial dimensions: $d=2$.

The external electric field $E_{i}$ and thermal drive $\zeta_{i}$ are added by modifying the background vector potential $A$ and spacetime metric $g$ [1]:

[TABLE]

We are interested in the thermoelectric conductivities within linear response, and so we need only calculate the perturbations \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}\mu, \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}T and \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}v_{i} to linear order in $E_{i}$ and $\zeta_{i}$. After some algebra, the linearized hydrodynamic equations can be found:

[TABLE]

where $\bar{\nabla}$ is the covariant derivative with respect to the spatial component of the metric $g_{ij}$. These are elliptic differential equations which can be straightforwardly solved numerically, as we describe in Appendix B.

**

Perturbative Analytic Solution**

In the limit where $g_{ij}$ is a perturbatively small deviation from flat space $g_{ij}=\delta_{ij}+\hat{g}_{ij}$, and the spatial variation of the chemical potential around the average $\mu_{0}(\mathbf{x})=\bar{\mu}_{0}+\hat{\mu}(\mathbf{x})$ is also perturbatively small, we can analytically compute the conductivity matrix to leading order. The calculation is rather tedious and is presented in Appendix C.333Note that $\hat{g}_{ij}$ is quadratic in the height function of out-of-plane distortions. For a chemical potential and induced metric with an explicit small parameter $u$, $\mu_{0}(\mathbf{x})=\bar{\mu}_{0}+u\hat{\mu}(\mathbf{x})$ and g_{ij}=\text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}_{ij}+u\partial_{i}\hat{h}\partial_{j}\hat{h}, with $\hat{\mu}$ and $\hat{h}$ O(1) functions. Hence fluctuations in the height function $\sqrt{u}\hat{h}$ must be parametrically larger amplitude than the fluctuations in the chemical potential, $u\hat{\mu}$, to have the same effect on transport. When $u\rightarrow 0$, the transport coefficients $\sigma_{ij}$, $\alpha_{ij}$ and $\bar{\kappa}_{ij}$ will be $\mathrm{O}(u^{-2})$. The transport coefficients can be expressed in terms of the relaxation rate $\tau^{-1}_{ij}$ for momentum. Assuming the density $n$ is finite, one expects on general grounds [1]:

[TABLE]

where

[TABLE]

with

[TABLE]

We have defined the projector

[TABLE]

The pure charge disorder scattering rate $(\tau_{ij}^{-1})^{(\mu\mu)}$ was found before in [7]. If the disorder in the chemical potential is uncorrelated with the strain disorder, then after disorder averaging we expect $(\tau^{-1})^{(\mu h)}\approx 0$.

As we have stated, $\tau^{-1}_{ij}$ is the rate at which the fluid can relax its momentum on the long wavelength disorder. As this is a dissipative process, it is necessarily the case that $\tau^{-1}_{ij}$ must depend on $\sigma_{\textsc{q}}$, $\eta$ and/or $\zeta$. As we show in Appendix C, in this perturbative regime the charge and heat currents are – at leading order – uniform. When there are inhomogeneities in the chemical potential, momentum relaxation can be non-negligible, even in the limit where inhomogeneity is very long-wavelength ($k\rightarrow 0$). This is due to the fact that a uniform fluid velocity, uniform charge and uniform heat current are not simultaneously consistent with both charge and heat conservation – the heat and charge currents must contain a purely dissipative component, which carries no momentum. The conductivity associated with this incoherent current is proportional to $\sigma_{\textsc{q}}$; this explains the $1/\sigma_{\textsc{q}}$ scaling of the first term in $(\tau^{-1})^{(\mu\mu)}$. A non-relativistic avatar of this effect was emphasized in [2]. However, in the presence of strain or metric disorder, there is no such impediment to a uniform flow. In this case, dissipation arises due to viscous effects. As viscosity vanishes, dissipation becomes weaker and hence the thermoelectric conductivity matrix is proportional to $1/\eta$.

**

Numerical Solution**

In order to make more precise predictions, the theory introduced in the previous section ought to be supplemented with specific equations of state. As noted in [7], the equations of state for a quasirelativistic fluid with gapless excitations, such as the Dirac fluid in graphene, are rather constrained. If we focus on the physics near the charge neutrality point for simplicity, we find:

[TABLE]

where the constants $\sigma_{0}$, $\eta_{0}$ and $C_{0,2}$ are dimensionless. For simplicity we have assumed that the bulk viscosity $\zeta=0$; we did not find that a finite $\zeta$ led to qualitatively different physics than a finite $\eta$.

Using the spectral methods of [7], described in Appendix B, we have numerically solved (14) with the equations of state (19), in inhomogeneous chemical potentials and metrics. We have always taken periodic boundary conditions, and assumed that the metric disorder and chemical potential disorder are uncorrelated, for simplicity.

Denoting spatial averages with $\mathbb{E}[\cdots]$, let us define

[TABLE]

These two parameters quantify the relative amount of charge vs. strain disorder. The overall prefactors of temperature $T$ are chosen so that $\epsilon_{\textsc{c,s}}$ are dimensionless numbers. In Figure 1, we demonstrate quite clearly the dramatic effects of viscosity on transport in the presence of strain disorder, as explained in the previous section.

The other dissipative channel is the one controlled by the microscopic conductivity $\sigma_{\textsc{q}}$. The presence of $\sigma_{\textsc{q}}$ is essential: for vanishing $\sigma_{\textsc{q}}=0$ there can be no heat current in the absence of an electric current, and so $\kappa=0$. So a clear way to observe the effects of $\sigma_{\textsc{q}}$ is in the Lorenz ratio

[TABLE]

where for simplicity we have assumed isotropic transport coefficients (this is the case for isotropic disorder). For perturbatively small disorder, we estimate the Lorenz ratio [7]

[TABLE]

as can be seen from the analytic results. For chemical potential disorder, we expect that (at small $\sigma_{\textsc{q}}$) $\tau\sim\sigma_{\textsc{q}}$ and so $\kappa$ does not depend strongly on $\sigma_{\textsc{q}}$. When all disorder is in the strain, $\tau$ does not depend on $\sigma_{\textsc{q}}$ and so $L$ has much stronger dependence on $\sigma_{\textsc{q}}$. This is shown in Figure 2.

The numerical results in Figure 1 and Figure 2 are still fully in the perturbative analytic regime. For larger disorder the analytic results are no longer quantitatively correct, even though the differences remain small and the qualitative features stay the same. This is shown in Figure 3.

A clear indication that one is outside the perturbative regime is that the results can no longer be described in terms of a sum of inverse scattering times. This is depicted in Figure 4. Beyond the perturbative regime, we find that the analytic expression overestimates the conductivity in the presence of strain disorder, and underestimates the conductivity in the presence of charge disorder.

**

Application to Suspended Graphene**

We now turn to the application of our formalism to hydrodynamics in the Dirac fluid in monolayer graphene [21, 22]. Graphene is a honeycomb lattice of carbon atoms in two spatial dimensions, with a low energy dispersion relation

[TABLE]

The $a$ label denotes spin and valley indices, and will mostly be ignored for the purposes of this paper – neither the interactions nor the disorder couples to spin here. These electrons interact with one another via long-range Coulomb interactions. Thus, strictly speaking, the hydrodynamics of graphene cannot be relativistic hydrodynamics.

However, as we have seen, transport is a linear response calculation. The key input from relativistic hydrodynamics was that the energy current and momentum density were identical – this reduced the number of hydrodynamic variables present. This follows trivially from the (weak-coupling) action for the Dirac fluid, and so we expect that the non-relativistic nature of the interactions will not play an important role in a transport calculation. Furthermore, as one can show following [30, 7], the effect of Coulomb interactions can be absorbed into a (nonlocal) redefinition of $\mu_{\mathrm{eq}}$ and \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}\mu, and so the final equations governing transport remain unchanged. Some of the literature also includes a long-lived (but not exactly conserved) imbalance mode into the hydrodynamic description [31, 32, 33]; for simplicity, we have not accounted for this effect. Indeed, the predictions of relativistic hydrodynamics have been confirmed experimentally in [7, 9]; see also [34, 35].

The key advance for the observation of hydrodynamic behavior was the growth of high quality graphene crystals, sandwiched between layers of hexagonal boron nitride. This dramatically reduced the size and the number of “charge puddles” – local inhomogeneity in the chemical potential [24]. As a consequence the disorder in graphene became weak enough that hydrodynamic effects were observable at $T\lesssim 100$ K. (When $T\gtrsim 100$ K, electron-phonon coupling appears to significantly degrade the electronic energy and momentum, and hence hydrodynamic behavior).

Another possibility for limiting the amount of disorder in graphene is to suspend it [25, 26]. The charge puddles in suspended graphene are also inherently quite weak. However, suspending graphene leads to a new source of disorder: flexural (out-of-plane) distortions of the graphene crystal. As we noted in the introduction, these distortions lead to an effective spatial metric $g_{ij}$ given by (1). In the limit where electron-electron interactions are negligible, these flexural modes are known to dominate the resistivity at low temperatures [36]. Our goal is to understand the implications of these flexural distortions on transport in suspended graphene, in the hydrodynamic limit. As we have already shown the consequences of (1) on transport, our goal here is simply to estimate the size of $h(x,y)$ in suspended graphene, and to comment on whether the hydrodynamic approximation is ever sensible.

In this paper, we will account for the flexural modes by considering motion on a curved space. In the limit where there are well-defined quasiparticles, it is common to interpret the strain not as the metric deformations (1) but as emergent magnetic fields [37, 38]. A priori, this is quite subtle – a magnetic field breaks time-reversal symmetry, while (14) preserves time-reversal symmetry. The resolution to this puzzle is that there are two Dirac points in the Brillouin zone in graphene, and the emergent magnetic field has opposite signs in each valley. The Dirac fluid of graphene, accounting for both valleys, will remain invariant under time-reversal in the presence of strain. Nonetheless, as we mentioned previously, strain can open a gap, and so it may be possible that in graphene the presence of strain leads to modifications of the effective hydrodynamics. These are questions worth considering more carefully in future work.

With these caveats, let us nonetheless estimate the hydrodynamic momentum relaxation rate due to long wavelength flexural fluctuations in graphene.

**

Classical Flexural Phonon Dynamics**

The classical action describing flexural phonons in graphene is [36]

[TABLE]

where $h$ is the height of the graphene membrane at position $(x,y)$. The parameters $\kappa\sim 1$ eV and $\rho\sim 7\times 10^{-7}\;\mathrm{kg/m}^{2}$ [39]. Assuming a square membrane of size $L$, and writing

[TABLE]

with allowed wave vectors \mathbf{k}=2\text{\mbox{\hskip 5.70027pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\pi}}}}/L\times(n_{x},n_{y}), we obtain

[TABLE]

As expected, we find a set of decoupled harmonic oscillators with

[TABLE]

In quantum mechanics, phonons are quantized, and so we should check the length scales at which this classical description will fail. This occurs when the occupation number of a given phonon mode is comparable to 1, which occurs when $\hbar\omega_{\mathbf{k}}\sim k_{\mathrm{B}}T$. This occurs when

[TABLE]

For the remainder of this section, we will restore factors of $\hbar$, $k_{\mathrm{B}}$ etc. The hydrodynamic description fails when

[TABLE]

At any reasonable experimental temperature, there are a classically large number of thermally excited flexural phonons at wave numbers in the hydrodynamic regime. We also learn that there is a large range of wave numbers where the phonons cannot be treated hydrodynamically. Hence, we expect a further contribution to momentum relaxation due to these higher wave number phonons, which must be computed using a more microscopic description, such as kinetic theory.

Next, we must ask whether or not the classical dynamics of flexural phonons is slow enough that the background metric may be treated as static. The fastest phonon dynamics in the hydrodynamic regime occurs for fluctuations $h_{\mathbf{k}}$ with $\mathbf{k}$ of order (29). Plugging into (28) we see that the fastest phonon dynamics in the hydrodynamic regime occurs at a rate

[TABLE]

Hence, the metric configuration $h(x,y)$, on hydrodynamic length scales, is essentially frozen in place on electronic time scales, justifying the assumption in our previous hydrodynamic analysis that the background geometry is time-independent.

**

Contribution to Momentum Relaxation Time**

We now compute the contribution of long wavelength fluctuations to the relaxation time for momentum. First, we must compute the typical size of thermal fluctuations in $h_{\mathbf{k}}$. Using the classical equipartition theorem, and recalling that $h_{\mathbf{k}}$ contains two independent harmonic oscillators (real and imaginary part):

[TABLE]

We have once again reverted to natural units. A straightforward computation, presented in Appendix C.1, reveals that

[TABLE]

The hydrodynamic result can only be trusted until for $\xi\lesssim 1/T$, and so we estimate that the contribution of (hydrodynamically) long wavelength flexural phonons to the momentum relaxation time is

[TABLE]

Near the charge neutrality point, the thermodynamic prefactors scale with known powers of temperature, and we obtain $\tau^{-1}\sim T^{3}$.

Of course, this must be compared with the other contributions to the momentum relaxation time, including the scattering off of short wavelength phonons. Using kinetic theory, this has been estimated to be $\tau^{-1}\sim T^{2}$ [39, 40]. Typically one would account for electron-phonon scattering using kinetic theory, treating each electron-phonon scattering event as a rare and independent process. However, we have just seen that in the hydrodynamic limit, a classical electron fluid with many electron-electron collisions moves in an approximately frozen phonon background. Thus, one electron can be correlated with the same phonon for over many collisions. These correlations suggest that the molecular chaos assumption underlying the kinetic description (that scattering events are uncorrelated with each other) could easily break down.

Additional mechanical strain induced by the contacts in a realistic sample of graphene changes the low frequency dispersion relation of flexural modes from quadratic to linear [39, 40]. Such a change would alter (31). But, from the form of (32) it is clear that the smallest wavelength phonons are most efficient at relaxing momentum. Hence, so long as the quadratic dispersion relation is restored by $k\sim\ell_{\mathrm{ee}}^{-1}$, we expect that (32) approximately accounts for the hydrodynamic contribution to the electron-phonon momentum relaxation rate.

Depending on the nonlinear properties of an elastic membrane, there can be significant renormalization of the effective $\kappa$ which should be used in (31), due to thermal fluctuations [41]. This effect has been seen recently in molecular dynamics simulations [42] and in experiment [43]. In a very simple approximation, one estimates that $\kappa_{\mathrm{eff}}\sim\sqrt{TK}/q$ as $q\rightarrow 0$; $K$ is a constant associated with certain nonlinearities in the elastic free energy. If this renormalization is significant in the hydrodynamic regime, then we expect that the temperature scaling in $1/\tau$ would be reduced by a factor of approximately $T^{3}$.

Finally, we note that there are other phonon modes which we could account for. In particular, there are also in-plane deformations that naturally arise, where the point $x_{i}$ is displaced to $x_{i}+d_{i}(x)$. In the presence of both a fluctuating height $h(x)$ and displacement $d_{i}(x)$, the general expression for $g_{ij}$ is [20]

[TABLE]

The in-plane phonons of graphene are linearly dispersing, and so $\langle|d_{\mathbf{k}}|^{2}\rangle\sim k^{-2}$, in contrast to (31). However, the metric itself depends on $d$, not on $d^{2}$, and contains one fewer spatial derivative. Putting this together and generalizing the discussion of Appendix C.1, we estimate that $\tau^{-1}\sim T^{4}$. Hence flexural phonons are more important than longitudinal phonons in the hydrodynamic limit.

**

Conclusion**

In this paper, we have described the effects of inhomogeneous slowly varying strain on hydrodynamic transport in strongly correlated electron fluids. We have demonstrated that for a (quasi)relativistic system with only strain disorder, the conductivities depend only on the viscosity of the electronic fluid (at least when inhomogeneity is small).

The conventional theory of electron-phonon scattering estimates the relaxation rate by simply computing low-order Feynman diagrams. Such an approach is sensible when the mean free path is much larger than the wavelength of both the electrons and the phonons. However, it is plausible that in charge-neutral graphene and other strongly correlated electron fluids, the electronic mean free path is short compared to the wavelength of some phonons. Our hydrodynamic description is the appropriate description of scattering off of these long wavelength phonons, though we must bear in mind that there will inevitably be a larger number of shorter wavelength phonons, which are not entirely captured by our hydrodynamic model.

A large open problem involves extending the theory of transport beyond the hydrodynamic limit. In the limit of weak interactions, this can be achieved using kinetic theory: while challenging, it is possible to completely characterize the ballistic-to-hydrodynamic crossover in this limit [14]. It would be interesting to understand how the hydrodynamic limit of electron-phonon coupling that we have demonstrated in this work can be understood from such a kinetic theory framework.

Acknowledgements

We are grateful to Vadim Cheianov and Philip Kim for useful discussions. VS and KS were supported in part by a VICI (KS) award of the Netherlands Organization for Scientific Research (NWO), by the Netherlands Organization for Scientific Research/Ministry of Science and Education (NWO/OCW), and by the Foundation for Research into Fundamental Matter (FOM). AL was supported by the Gordon and Betty Moore Foundation’s EPiQS Initiative through Grant GBMF4302.

Appendix A Onsager Reciprocity

In this appendix we show that Onsager reciprocity is satisfied on a curved background. This is a non trivial consistency check of our formalism, as it has to be satisfied for any time-reversal symmetric theory of transport.

We begin by introducing some shorthand notation for our proof, following [4]. We denote a uniform spatial average with $\mathbb{E}[X]=\int\frac{\mathrm{d}^{d}x}{L^{d}}\sqrt{g}X$ , where $g$ is the determinant of the spatial metric $g_{ij}$. Define the vectors

[TABLE]

It is straightforward to see that Eqs. (14) are equal to

[TABLE]

We have denoted

[TABLE]

To prove Onsager reciprocity we must prove that

[TABLE]

By linearity, we may write the solutions to these equations of motion as

[TABLE]

We have denoted the index $J$ in capital letters to emphasize that $\Phi^{J}$ is not a contravariant vector, and that $v^{iJ}$ is a contravariant vector, not a tensor. (36) then becomes

[TABLE]

Now, by definition

[TABLE]

We now use the left hand side of (40) to re-write the coefficients of the first two terms of (41), after integrating by parts the second term of (41) (note that $\mathbb{E}[X\nabla_{i}Y^{i}]=-\mathbb{E}[Y^{i}\nabla_{i}X]$ – this can be easily seen using $\nabla_{i}Y^{i}=g^{-1/2}\partial_{i}\left(g^{1/2}Y^{i}\right)$):

[TABLE]

We next integrate by parts the viscous part of the first term, as well as the second term. Using the fact that $\mathbb{E}[Y_{i}\nabla_{j}Z^{ij}]=-\mathbb{E}[Z^{ij}\nabla_{j}Y_{i}]$, this can easily be done. We find

[TABLE]

It is straightforward to determine from this expression, as well as the symmetry properties of $\eta_{ijkl}$, that $\sigma^{\alpha\beta IJ}=\sigma^{\beta\alpha JI}$, proving Onsager reciprocity.

Appendix B Numerical Simulations

The hydrodynamic equations can be solved on a periodic domain of $L\times L$ grid points with pseudospectral methods [44]. The domain is reduced to a set of uniformly distributed grid points $\mathbf{x}_{I}$ and the continuous partial differential equation (14) is approximated by the discrete form

[TABLE]

where

[TABLE]

is the data we wish to solve for,

[TABLE]

is the source for the linearized hydrodynamic equations due to the external fields $E_{i}$ and $\zeta_{i}$ (henceforth, the dependence on $\mathbf{x}_{I}$ of coefficients is implicit), and $\mathsf{L}$ is the discretized differential operator

[TABLE]

The derivative operators in this equation are understood to act on everything to their right, unless contained within parentheses. They are approximated on the discrete grid using spectral methods. By inverting444The operator $\mathsf{L}$ has two zero eigenvectors, corresponding to a constant shift of $\mu$ and $T$. This can be fixed eliminating an appropriate pair of rows/columns from the matrix $L$, so as to enforce the constraints $\mu(\mathbf{0})=T(\mathbf{0})=0$. the system (44), we find expressions for \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}\mu, \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}T and \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}v_{i}, from which it is simple to compute the spatially averaged charge and heat currents, and hence the thermoelectric conductivity matrix.

The disorder consists of random sine waves both for the chemical potential and for the strain:

[TABLE]

where $\bar{\mu}_{0}$ is constant, $\phi_{x,y}$ is uniformly distributed on [0,2\text{\mbox{\hskip 5.70027pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\pi}}}}) and $\hat{\mu_{0}}(n_{x},n_{y})$, $\hat{h}(n_{x},n_{y})$ are uniformly distributed on $[-c,c]$, with c=\sqrt{(2-\text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}_{n_{x},0}-\text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}_{n_{y},0})/2}. We have chosen this expression for $c$ such that both charge disorder and strain are not altered by a random homogeneous offset. $k$ denotes the number of sine waves in a given direction. In all simulations we have used $L=8k+3$ grid points in the $x$ and $y$ directions.

The numerical simulations were performed using a finite number of disorder modes $k$, both for the charge disorder and for the strain disorder. This introduces another form of finite size effect to the results. It can be partially reduced by averaging over multiple disorder realizations. We also caution that due to the nonlinear way that strain couples to the hydrodynamic equations, an $h(\mathbf{x})$ which varies on the wavelength $\xi$ leads to a varying metric on the scale $\xi/2$, which can source the equations of motion at even shorter wavelengths due to the nonlinear metric dependence in covariant derivatives.

Despite these caveats, we expect that our numerical methods converge exponentially quickly upon increasing $L$ [44]. Simple tests confirm this [7]. We expect therefore that the results presented in the main text are a reasonable quantitative characterization of hydrodynamic transport in the presence of disorder. We have also found that the deviations in transport coefficients from one realization of disorder to the next is relatively mild: see Figure 5 for an example of the size of deviations from one realization of disorder to the next.

Appendix C Weak Disorder

In this appendix, we perturbatively calculate the transport coefficients assuming that the inhomogeneity is small. The calculation is analogous to [4, 7]. As we discussed in the main text, it is necessary for $\mathrm{\Delta}\mu\sim u$ and $\mathrm{\Delta}h\sim\sqrt{u}$, with $u\rightarrow 0$; this will lead to transport coefficients which are perturbatively large: $\sigma_{ij}\sim u^{-2}$. We will denote with $\hat{n}(\mathbf{k})$ the charge density at wave number $\mathbf{k}$, and denote the metric as

[TABLE]

with $\hat{g}\sim u$ perturbatively small. To leading non-trivial order, g^{ij}\approx\text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}^{ij}-\hat{g}^{ij}, with the indices raised trivially using the flat space metric.

Upon Fourier transforming \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}\mu, \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}T and \text{\mbox{\hskip 4.44444pt{1 0 -.18 1}\hbox to0.0pt{\hss\leavevmode\hbox{\set@color\delta}}}}v_{i}, we write

[TABLE]

We will show that to leading order in $u$, $\mu(\mathbf{k})\sim T(\mathbf{k})\sim v_{i}(\mathbf{k})\sim u^{-1}$ while $\bar{v}_{i}\sim u^{-2}$. Let us begin by writing down the $\mathrm{O}(u)$ equations for $\mu(\mathbf{k})$, $T(\mathbf{k})$ and $v^{i}(\mathbf{k})$. In these equations, indices may be raised and lowered freely, since all corrections are higher order in $u$. We obtain

[TABLE]

Upon solving for $T(\mathbf{k})$, $\mu(\mathbf{k})$ and $k_{i}v_{i}(\mathbf{k})$ at leading order, we obtain:

[TABLE]

It is also helpful to combine these equations to obtain

[TABLE]

We now spatially average the momentum conservation equation at $\mathrm{O}(u^{0})$, and find

[TABLE]

It is helpful to spatially average the last column of (47), without an overall metric factor, to obtain this equation, as many total derivative terms are shown explicitly. Next, we note that in all terms in this equation, we may raise and lower spatial indices with the flat space metric, so long as we keep in mind

[TABLE]

Our goal is to eliminate $\mu(\mathbf{k})$, $T(\mathbf{k})$ and $v_{i}(\mathbf{k})$ in favor of $\bar{v}_{j}$. Given (52) and (53) this is straightforward but tedious. One finds that

[TABLE]

where

[TABLE]

From (56) we obtain an expression for $\bar{v}_{j}$. Using $J_{i}\approx n_{0}\bar{v}_{i}$ and $Q_{i}\approx T_{0}s_{0}\bar{v}_{i}$, we obtain (15) and (17).

**

Phonon Contribution**

Let us now evaluate the contribution of flexural phonons to the momentum relaxation time. In order to evaluate the momentum relaxation time, we need to evaluate thermal averages which can be written in the generic form

[TABLE]

In the second line, we have used (31) to remove thermal averages of $h_{\mathbf{k}}$, and in the third line we have used the fact that $k_{i}\mathbb{P}_{ij}=0$. Now, we need to evaluate the sum over $\mathbf{q}$. In general, we will find

[TABLE]

but due to the projectors coming from the formula for $\tau^{-1}$, we need only compute the coefficient $A$. We do so by assuming that $\mathbf{k}=k\hat{\mathbf{x}}$ and then set all indices $k=l=m=n=y$ in (59). In order to approximate the sum over $\mathbf{q}$ with an integral, we employ the standard substitution

[TABLE]

Hence, we find

[TABLE]

Hence, we find a finite value for $A$. In contrast to $A$, one can show that other coefficients in (59) do have logarithmic divergences. Given the value for $A$, we find

[TABLE]

Employing

[TABLE]

we find (32).