Spin-caloric resistance of Dirac plasma in a graphene Corbino device

Thermoelectric (TE) and thermomagnetic (TM) transport effects in metals reflect particle-hole asymmetry and arise from the coupling of charge and heat currents Abrikosov (1988). A different class of thermomagnetoelectric (TME) effects emerges when one considers coupling of charge and heat flow with the spin current Aronov (1976); Johnson and Silsbee (1987). The distinction between TM and TME phenomena lies in the former being driven by the Lorentz force of the applied external magnetic field exerted on the charge flow, while the latter arises from the flow of spins driven by a magnetic field.

In metallic magnetic heterostructures, such as spin valves, spin-dependent thermoelectric effects can be understood through two main models. The first model involves two parallel spin-transport channels with distinct thermoelectric properties. Alternatively, these effects can be attributed to collective phenomena, influenced by spin waves, which are also present in insulating ferromagnets. For further references, refer to the collection of papers and reviews on spin caloritronics Bauer et al. (2010, 2012). The two-channel model, incorporating spin-dependent conductivities, can be derived from the kinetic equation. This model assumes different drift velocities for two spin subsystems and is particularly suitable in the weak interaction limit. Essentially, it extends the single-particle Drude-Sommerfeld theory to account for spin-polarized transport. Consequently, it adheres to Mott’s law for the Seebeck coefficient and the Wiedemann-Franz law for spin-dependent thermal conductivity. In contrast, in the collision-dominated regime, the drift velocity is the same for all the spin components. The reason is that frequent collisions between electrons with different spins form a common drift of the electron system. This scenario can be realized in high-mobility and low-density semiconductor and semimetal devices, where strong correlations lead to the electron liquid reaching a hydrodynamic limit Gurzhi et al. (2006); Andreev et al. (2011). In this situation, the coupling of spin and heat currents exhibits a distinct character, leading to the specific mechanism of thermal resistance explored in this paper.

The hydrodynamic behavior of the electron liquid has been demonstrated in a number of recent experiments, see reviews Lucas and Fong (2018); Levchenko and Schmalian (2020); Narozhny (2022); Fritz and Scaffidi (2024) on this topic and references therein. Of particular relevance to this work are strong thermoelectric anomalies observed in graphene devices tuned to the proximity of charge neutrality Crossno et al. (2016); Ghahari et al. (2016). Physically, these profound features, namely giant enhancement of the Lorenz ratio and violation of the Mott relation by inelastic scattering, stem from the decoupling of charge and heat flows.

The hydrodynamic equations in the presence of finite spin-polarization of an electron liquid induced by an in-plane magnetic field should be supplemented by the conservation law for spin. This is particularly relevant to graphene multilayers, as the spin relaxation time is known to be long, even at room temperature, in the range of nanoseconds, with corresponding spin diffusion lengths reaching 10 $\mu$m Drögeler et al. (2014, 2016). In the linear response, the spin current consists of a spin diffusion current and a spin thermocurrent induced by temperature gradients. According to Onsager’s reciprocity principle Onsager (1931), the heat flux acquires an additional spin-dependent contribution related to the Peltier effect. These processes modify the viscous flow of the electron fluid, thus influencing dissipative friction and consequently strongly increasing thermal resistance.

To explore this physics, a particularly useful device geometry is a graphene Corbino membrane, as illustrated in Fig. 1. In this setup, the continuity of currents enforces the electronic flow profile to be spatially nonuniform, even at uniform particle and entropy density. This amplifies the interplay of dissipative viscous, spin-diffusion, and spin-thermal effects. It should be noted that charge and thermal transport measurements in gated graphene Corbino devices have been reported, including observations of viscous and magneto-thermal effects Yan and Fuhrer (2010); Faugeras et al. (2010); Waissman et al. (2024); Talanov et al. (2024); Zeng et al. (2024). Additionally, powerful thermal and electrochemical potential imaging techniques can be implemented Waissman et al. (2022); Kumar et al. (2022); Krebs et al. (2023); Völkl et al. (2024). Thus far, theoretical works have focused on the orbital effects of the magnetic field in the magnetotransport Shavit et al. (2019); Kavokin et al. (2020); Levchenko et al. (2022); Gall et al. (2023), while hydrodynamic effects dominated by spin thermotransport remain unexplored, which motivates this work.

To this end, consider a two-dimensional electron system subjected to an in-plane magnetic field. In the regime where the rate of momentum-conserving electron-electron collisions exceeds the rates of momentum, energy, and spin relaxation, the corresponding macroscopic hydrodynamic equations of the electron fluid can be written as follows:

Here $P$ is the pressure, $\bm{u}$ is the hydrodynamic velocity, $k$ denotes the disorder-induced friction coefficient, and $\hat{\bm{\Sigma}}\equiv\Sigma_{ij}$ denotes the viscous stress tensor Landau and Lifshitz (1987)

with $\eta$ and $\zeta$ being, respectively, shear (first) and bulk (second) viscosity of the electron liquid. The electric potential $\phi$ is related to the electron charge density $en$ by the Poisson equation. The densities of particles, entropy, momentum, and spin are denoted by $n$, $s$, $\bm{p}$, and $\sigma$ respectively. The corresponding current densities of particles, entropy, and spin are labeled respectively as $\bm{j}$, $\bm{j}_{s}$, and $\bm{j}_{\sigma}$. The local rate of entropy production due to electron-electron collisions is denoted by $\dot{s}$ in Eq. (1b) Likewise, $\dot{\sigma}$ in Eq. (1c) corresponds to the spin relaxation. The continuity equation for the spin current is written in analogy to the hydrodynamic theory of spin waves Halperin and Hohenberg (1969). It should be noted that only spin component along the field appears in this framework as the other spin components are not conserved due to spin precession.

We proceed to determine the entropy production rate. For this purpose, we recall that the entropy density $s$ is a function of the conserved quantities, and its differential is given by the general thermodynamic relation Landau and Lifshitz (1980)

where $\mu$ and $\mu_{\sigma}$ are chemical potentials for particle and spin. Noting that all these thermodynamic variables correspond to conserved quantities, we express their time derivatives in terms of the divergences of the corresponding currents from Eq. (1). This yields the following expression for $\dot{S}=\int\dot{s}d\bm{r}$,

where we used $\partial_{t}\varepsilon=-\bm{\nabla}\cdot\bm{j}_{\varepsilon}$ for the energy current density $\bm{j}_{\varepsilon}$, and similarly for $\partial_{t}n$, $\partial_{t}\sigma$, and $\partial_{t}\bm{p}$ from Eqs. (1a)-(1c). The tensor $\hat{\bm{\Pi}}\equiv\Pi_{ij}=P\delta_{ij}-\Sigma_{ij}$ denotes the local part of the momentum flux. Next, using the thermodynamic identity $\bm{\nabla}P=n\bm{\nabla}\mu+s\bm{\nabla}T+\sigma\bm{\nabla}\mu_{\sigma}$ and integrating by parts we get

where we introduced dissipative fluxes $\bm{j}^{\prime}_{s}=\bm{j}_{s}-s\bm{u}$, $\bm{j}^{\prime}=\bm{j}-n\bm{u}$, and $\bm{j}^{\prime}_{\sigma}=\bm{j}_{\sigma}-\sigma\bm{u}$ that are proportional to gradients of the equilibrium parameters. The entropy flux is related to the energy flux as follows $\bm{j}_{s}=\frac{1}{T}[\bm{j}_{\varepsilon}-(\mu+e\phi)\bm{j}-\mu_{\sigma}\bm{j}_{\sigma}]$.

In a situation with a vanishing electromotive force, $e\bm{\mathcal{E}}=-\bm{\nabla}(\mu+e\phi)$, the second term in the above expression is absent since $\bm{j}^{\prime}=0$. Recalling the definition $\dot{S}=-\sum_{i}X_{i}\dot{x}_{i}$ from Ref. Landau and Lifshitz (1980) and identifying $\dot{x}_{i}$ with the entropy current $\bm{j}^{\prime}_{s}$ and the spin diffusion current $\bm{j}^{\prime}_{\sigma}$ we conclude that the role of the corresponding thermodynamic forces is played by $\bm{\nabla}T/T$ and $\bm{\nabla}\mu_{\sigma}/T$ respectively. Introducing the kinetic coefficients following the standard convention $\dot{x}_{i}=-\gamma_{ij}X_{j}$ we write

We can express $\gamma_{11}=\kappa$ in terms of the thermal conductivity $\kappa$, and $\gamma_{22}$ in terms the spin diffusion $D_{\sigma}$ as $\gamma_{22}=TD_{\sigma}$. Since the spin density changes sign with respect to time reversal symmetry and energy does not, the Onsager symmetry principle gives $\gamma_{12}(H)=-\gamma_{21}(-H)$. On the other hand both $\gamma_{12}(H)$ and $\gamma_{21}(H)$ are odd functions of the magnetic field $H$. Therefore, we conclude that $\gamma_{12}(H)=\gamma_{21}(H)$. Thus the matrix of kinetic coefficients is symmetric. In the following we denote $\gamma_{12}\equiv\gamma_{\sigma}$ and introduce

As a reminder, for temperatures $T<E_{F}$ in the Fermi liquid theory $\kappa\sim E^{2}_{F}/T$ Abrikosov and Khalatnikov (1959) and $D_{\sigma}\sim A/T^{2}$ Hone (1961), modulo logarithmic renormalizations specific to the 2D case, see e.g. Ref. Lyakhov and Mishchenko (2003). A detailed analytical calculations of these coefficients in the limit of weak interactions, $r_{s}=V/E_{F}\lesssim 1$, where $V$ is the characteristic energy of Coulomb interaction and $E_{F}$ is the bare Fermi energy, can be found in Ref. Sykes and Brooker (1970).

For compactness of further expressions, it is convenient to introduce the column-vector notations

Then it can be readily checked that the sum of the first and third term in the right hand side of Eq. (II) can be rewritten as $\vec{\bm{X}}^{\mathbb{T}}\hat{\Xi}\vec{\bm{X}}$, with the superscript $\mathbb{T}$ denoting the transposition. As a result,

where summation over the repeated indices is implicit. The current densities take the final form in explicit notations

where the first term in the right hand side represents the equilibrium components of the currents, while the remaining two terms in each expression represent the dissipative components. With the appropriate boundary conditions for a given device geometry, the system of equations Eqs. (1), (9), and (10) enable determination of the linear response transport coefficients (see Appendix A for additional discussion).

We apply the above formalism to consider heat transport in a graphene flake shaped like a Corbino disk, see Fig. 1. The aspect ratio of the disk is defined by the ratio of its radii, $p=r_{2}/r_{1}>1$. In the linear response regime, a small heat current is induced in the device by a temperature difference $\Delta T=T_{1}-T_{2}$ between the electrodes. The calculations presented below will focus on the regime of charge neutrality, where $n\to 0$. For simplicity, we also ignore the effect of bulk viscosity, which is known to vanish for both parabolic and linear spectra of quasiparticle excitations Pitaevskii and Lifshitz (1981). The latter corresponds to the case of graphene monolayer.

In a steady state, the continuity equation of entropy and spin currents can be presented in the form specific to the cylindrical geometry Li et al. (2022); Levchenko et al. (2022)

where $u(r)$ is the radial hydrodynamic velocity and $\vec{X}(r)$ is the column-vector of thermodynamically conjugate forces, with $r\in[r_{1},r_{2}]$ being the radial coordinate. The force balance condition exerted on an element of the fluid takes the form of the Navier-Stokes equation, which can be found with the help of Eqs. (1d) and (2):

Equations (11) and (12) determine the spatial dependence of the flow velocity, temperature gradients, and gradients of the spin chemical potential in the interior of the Corbino disk in terms of the entropy and spin currents $\vec{I}$. The crucial point to realize here is that even though thermal resistance is defined under the global condition of vanishing spin current, locally dissipative parts of these currents are nonzero due to finite $\vec{X}(r)$ in the bulk of the flow. This can be seen explicitly.

For this purpose, we use Eq. (11) to exclude $\vec{X}(r)$ and arrive at the equation for $u(r)$ in the form

Here we introduced the characteristic length scale $l$ defined by

The general solution of Eq. (13) consists of the linear superposition of solutions to the homogeneous equation and a particular solution due to the right hand side. The former is given by a linear combination of modified Bessel functions of the first and second kinds. These functions describe deviations of the hydrodynamic flow from that in the bulk of the disk. Due to properties of Bessel functions they are exponentially localized near the inner and outer boundaries on the length scale of the order of $l$. For the graphene monolayer this length scale is of the order of thermal de Broglie length $l\sim l_{T}=v/T$, which in itself is of the order of electron mean free path $l_{\text{ee}}$. Therefore, these solutions contribute to the resistance of the contacts. Consideration of these effects goes beyond the limit of hydrodynamic theory. We are interested in the solution that describes the hydrodynamic mode in the bulk of the device. This solution can be easily found as a particular solution of the inhomogeneous equation, yielding the result

The factor defined by $f$ can be safely set to unity. Indeed, viscosity of the Dirac fluid is $\eta\sim(T/v)^{2}\sim l^{-2}_{T}$ Müller et al. (2009). In addition, friction coefficient can be related to the spatial average of the local density fluctuations $\delta n(\bm{r})$ of the globally charge neutral system, namely $k\sim\frac{e^{2}}{\varsigma}\langle\delta n^{2}\rangle$ Li et al. (2020), where $\varsigma$ is the intrinsic conductivity of graphene Mishchenko (2008); Fritz et al. (2008); Kashuba (2008). Therefore, $kl^{2}/\eta\sim\langle\delta n^{2}\rangle l^{4}_{T}\ll 1$. Based on this argument, we set $f\to 1$ in what follows. Knowledge of $u(r)$ enables us to compute dissipative friction term and also the viscous term in Eq. (9). For a purely radial flow, there are only two nonvanishing components of the stress tensor that contribute to the dissipated power. In cylindrical coordinates, these are

As the next step, we need to determine $\vec{X}(r)$, which governs contribution of the relative transport mode to the dissipation rate in Eq. (9). The required solution can be obtained by multiplying Eq. (11) by a row vector $\vec{\varkappa}^{\mathbb{T}}=(\sigma,-s)$ from the left. The result is

which completes solution of the hydrodynamic equations.

We introduce spin caloric resistance matrix $\hat{\mathcal{R}}$ that can be determined by equating the Joule heat,

to the rate of energy dissipation in the bulk flow as defined by Eq. (9). Here $I_{q}=TI_{s}$ is the heat current. The thermal resistance $R_{\text{th}}$ can be found by setting the spin current to zero $I_{\sigma}\to 0$. The computed dissipation in Eq. (18) is then simply $\mathcal{P}=\mathcal{R}_{11}I^{2}_{q}$. On the other hand, employing general thermodynamic relations, one finds the entropy production rate $\dot{S}=I_{q}\Delta T/T^{2}=R_{\text{th}}I^{2}_{q}/T^{2}$. As a consequence, from the condition $\mathcal{P}=T\dot{S}$ we determine that

We insert Eqs. (15), (29), and (17) into Eq. (9), set $I_{\sigma}\to 0$, compute each of the radial integrals individually, and bring the resulting expression into Eq. (19) that will cancel $I^{2}_{s}$ in the denominator. The resulting formula can be written as a sum of respective contributions (see Appendix B for technical details)

The first term describes viscous dissipation modified by the spin diffusion and spin thermocurrent, which reads

Eqs. (46) and (21) represent the main results of this work. Next we analyze several limiting cases.

In pristine systems, $k\to 0$, and without magnetic field, $\sigma\to 0$, the result simplifies to $R_{\text{th}}=\eta(p^{2}-1)/\pi Ts^{2}r^{2}_{2}$ (see Appendix C for an alternative derivation). For a graphene monolayer $\eta\sim s\sim(T/v)^{2}$, so that $R_{\text{th}}\propto 1/T^{3}$. In contrast, in bilayer graphene (BLG) the entropy density has different temperature dependence, $s\sim m^{*}T$, where $m^{*}$ is the effective mass of the band structure. From the bound on the viscosity to entropy ratio, known from the strongly interacting quantum field theories Kovtun et al. (2005), one can deduce that $\eta\propto T$, therefore $R_{\text{th}}\propto 1/T^{2}$ in BLG. These predictions can be tested in experiments.

The leading field dependent correction to $R_{\text{th}}$ comes from $R_{\Xi}$. At low magnetic fields, we write the spin density in the form $\sigma=\chi H$, where $\chi$ denotes the spin susceptibility. As a result, for the relative thermal magnetoresistance, one obtains

which is controlled by the spin diffusion coefficient. We conclude that the spin caloric effect enables the extraction of $D_{\sigma}$ and, thus, the spin conductance via the corresponding Einstein relation, in the transport regime dominated by electron collisions. This approach could offer valuable complementary experimental tools in addition to earlier methods based on spin Coulomb drag measurements Weber et al. (2005).

These results can be readily generalized to the magnetotransport away from charge neutrality. There are three key modifications that need to be incorporated. (i) At finite particle density $n$, the thermoelectric effects couple charge and heat flows. Therefore, we have to augment Eq. (10) by the respective expression for the electric current density

where $\gamma$ is the intrinsic thermoelectric coefficient and the electromotive force (EMF) $e\bm{E}$ is the gradient of the local electrochemical potential. It is important to realize that finite $\bm{E}$ is generated even in a purely thermal bias setup. Since thermal resistance is defined at vanishing spin and electric currents, vanishing of $\bm{j_{e}}$ implies a local relationship $\bm{E}=-(en/\varsigma)\bm{u}$. Note that at $n\ll s$ the thermoelectric contribution to $\bm{j}_{e}$ may be neglected because $\gamma/T\sim n/s\ll 1$. (ii) The Coulomb force $-en\bm{E}$ must be added to the right-hand-side of the force balance condition, Eq. (12). Inserting $\bm{E}$ into the force-balance condition gives a force term $(en)^{2}\bm{u}/\varsigma$, which amounts to the following redefinition of the effective friction coefficient in Eq. (12), $k\to k(n)=k+\frac{e^{2}}{\varsigma}n^{2}$. The form of the equation remains intact thus its solution for $\bm{u}$, which is still given by Eq. (15). (iii) An additional term must be added to the entropy production rate in Eq. (9), which takes a simple form

with the neglect of thermoelectric contribution $\propto\gamma/T\ll 1$. With these ingredients, we use Eq. (15) to calculate Eq. (24) and thus determine additional contribution to thermal resistance due to intrinsic electric conductivity

In this limit $f$ defined in Eq. (15) has noticeable density dependence. Therefore $R_{\Sigma}$ and $R_{k}$ should be also multiplied by a factor $f^{2}$, since $u\propto f$, but not $R_{\Xi}$, which remains the same. The limit of charge neutrality corresponds to $n\ll\sqrt{\langle\delta n^{2}\rangle}$. In modern encapsulated devices density inhomogeneity of monolayer graphene is in the range $\sqrt{\langle\delta n^{2}\rangle}\sim(5\div 10)\times 10^{9}$ cm^-2. This is still well within the domain of Dirac plasma $\{n,\sqrt{\langle\delta n^{2}\rangle}\}<s$ at $T\sim 100$ K. In this limit one may approximate $f\approx s^{2}/(s^{2}+\frac{e^{2}}{\varsigma}n^{2})$. Therefore, $R_{\Xi}$ remains the leading spin-dependent contribution to thermal resistance at all relevant densities $n<s$.

It is instructive to estimate the order of magnitude of computed thermal resistances. For this purpose consider Eq. (46) in the limit $\sigma\to 0$ as an example. Recall that for a monolayer graphene $\eta\sim s\sim(k_{B}T/\hbar v)^{2}$, where Boltzmann and Planck constants were restored to have proper units in SI system. Let us take $T\sim 100$ K and $v\sim 10^{6}$ m/s. The typical radius of the Corbino disk in experiments is $r_{1,2}\gtrsim 1\mu$m. For this values, let us convert the conventional K/W units of $R_{\text{th}}$ into $\Omega$ units of electrical resistance by using the Lorenz number $L$ as a conversion factor through the Wiedemann-Franz law $L=R_{\text{el}}/TR_{\text{th}}$. This gives $R_{\text{th}}\sim(\hbar/e^{2})1/(sr^{2})\sim 100\,\Omega$. The measured values of thermal resistance reported in Ref. Talanov et al. (2024) are in the range $\sim 50\div 600\,\Omega$ for temperatures in the range $T\sim 25\div 200$ K. The estimate matches well.

The kinematic viscosity measured in Corbino geometry as reported in Refs. Talanov et al. (2024); Zeng et al. (2024) is in the range $\nu\sim 0.1\div 0.5$ m²/s. The kinematic viscosity is given by the shear viscosity per mass density $\nu=\eta/\rho$. The mass density of the Dirac thermal cloud can be estimated as $\rho\sim T^{3}/v^{4}$. Putting this together and restoring proper conversion units one finds $\nu\sim\hbar v^{2}/k_{B}T$, which gives $\nu\sim 0.1$ m²/s for $T=100$ K. Again this gives good agreement.

It is harder to estimate the magnitude of the relative MR in Eq. (22) because of unknown values of spin diffusion and spin susceptibility in the hydrodynamic regime, but form the structure of the expression it is expected that the result is significant even at relatively weak field.

I am grateful to Anton Andreev and Philip Kim for insightful discussions. This work was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES) under Award No. DE-SC0020313 and H. I. Romnes Faculty Fellowship provided by the University of Wisconsin-Madison Office of the Vice Chancellor for Research and Graduate Education with funding from the Wisconsin Alumni Research Foundation.