Spin accumulation in the spin Nernst effect

Spin caloritronics is a subfield of spintronics and aims to convert waste heat into spin that carries information [1, 2]. The door to this new research field was opened by the celebrated discovery of the spin Seebeck effect [3, 4], in which the spin current is induced parallel to an applied temperature gradient. On the other hand, the spin Nernst (SN) effect [5] is a heat analog of the spin Hall (SH) effect [6], namely, a phenomenon in which the spin current flows perpendicular to the temperature gradient and then turns into spin at the boundaries, as depicted in Fig. 1(a). This phenomenon was theoretically proposed [7, 8, 9] and experimentally observed in heavy metals such as platinum and tungsten using SN magnetoresistance [10, 11, 12]. It will enrich candidate materials and device design for spincaloritronics.

The Landauer-Büttiker formalism is a powerful theoretical approach to the SH [13, 14] and SN effects [7, 8]. In this formalism, the spin current is measured only in the leads, where a spin-orbit coupling is assumed to be absent, and hence is well defined. The generalized Mott relation between the SH conductance [13, 14] and SN coefficient [7, 8] is obvious although not mentioned.

A complementary theory is based on the linear-response theory. However, there are two difficulties in the field-theoretical formulation of the SN effect. One is the ambiguity regarding the definition of the spin current. In the presence of a spin-orbit coupling, spin is not conserved, and the spin current is not well defined. We may use the conventional spin current, $\hat{J}_{sa}^{\phantom{sa}i}(\bm{k})=\{\hat{s}_{a},\hat{v}^{i}(\bm{k})\}/2$, but this current is unphysical because its uniform equilibrium expectation value is nonzero in noncentrosymmetric systems such as the Rashba and Dresselhaus models [15]. Here, $\hat{s}_{a}$ and $\hat{v}^{i}(\bm{k})$ are the spin and velocity operators, respectively. We may also use the conserved spin current [16, 17]. If the spin torque density $\tau_{a}(t,\bm{x})$ vanishes on average over the whole system, we can define the spin torque dipole density as $\tau_{a}(t,\bm{x})=-\partial_{x^{i}}P_{\tau a}^{\phantom{\tau a}i}(t,\bm{x})$, and $\tilde{J}_{sa}^{\phantom{sa}i}(t,\bm{x})=J_{sa}^{\phantom{sa}i}(t,\bm{x})+P_{\tau a}^{\phantom{\tau a}i}(t,\bm{x})$ is conserved on average. This current has a desirable property of Onsager’s reciprocity relation. However, there is no convincing way to define the spin current, since spin nonconservation is inevitable.

The other difficulty is magnetization correction. According to one of the first theoretical proposals [9], the Kubo formula of the SN conductivity of the conventional spin current diverges towards zero temperature. Such unphysical divergence is well known in the context of the Nernst and thermal Hall effects and can be corrected by orbital and heat magnetizations, respectively [18, 19, 20, 21, 22]. These corrections arise because Luttinger’s gravitational potential [23] perturbs the charge and heat current densities as well as the density matrix, leading to the Kubo formulas. In the absence of inelastic scattering, the correct Nernst and thermal Hall conductivities are related to the Hall conductivity via the generalized Mott relation and Wiedemann-Franz law, respectively. Similarly, the spin magnetic quadrupole moment is necessary for the gravitomagnetoelectric effect [24, 25], in which the magnetization is induced by a temperature gradient. Since the SN effect is a spin analog of the Nernst effect, that of orbital magnetization will be necessary for the SN effect.

So far, the SN conductivity of the conventional spin current has been computed assuming the generalized Mott relation in first-principles calculations [26, 27, 28, 29, 30, 31], but its justification remains unsettled. In Refs. [32, 33], using the quantum-mechanical linear-response theory, the authors introduced the spin-resolved orbital magnetization and obtained the SN conductivity that vanishes at zero temperature. More strongly, it was recently shown that the nonequilibrium part satisfies the generalized Mott relation using the semiclassical theory of wave packet dynamics [25] and the non-Abelian Keldysh formalism [34]. However, the equilibrium part does not always take the form of rotation, and hence its integral over a cross section does not vanish. Regarding the conserved spin current, the equilibrium part takes the form of rotation, and the nonequilibrium part satisfies the generalized Mott relation [35, 36]. These properties are desirable but not sufficient to define the spin current.

In the context of the SH effect, it is well understood that the essence is the spin accumulation, which is well defined theoretically and directly observable experimentally. This idea was pointed out already in the first theoretical proposal [37], and the spin accumulation has been calculated via the spin diffusion equation and Landauer-Keldysh formalism [38, 39, 40, 41, 42, 43, 44, 45, 46]. However, it is difficult to deal with finite geometries for realistic multiband systems with complicated spin-orbit couplings. Recently, a complementary theory was proposed considering the spin response to an electric-field gradient [47, 48, 49]. The electric field and resulting charge current must vanish outside the system and vary near the boundaries. The gradient captures the lowest order of such a spatial variation and is considered to induce spin accumulation in this theory. Since a nonuniform field can be described by the nonzero wave number under periodic boundary conditions, this theory can be implemented to first-principles calculations.

In this Letter, we study the thermal spin accumulation at the boundaries in the SN effect by calculating the linear response of spin to a gradient of a temperature gradient, as depicted in Fig. 1(b). We prove the generalized Mott relation assuming the presence of time-reversal symmetry and the absence of inelastic scattering only. The aforementioned difficulties can be avoided in this formulation; we do not rely on the ambiguous spin current, and no correction term is necessary because the spin accumulation is a dissipative Fermi-surface term. We apply this relation to the three-dimensional isotropic Luttinger model [50] for $p$-type semiconductors and the two-dimensional Rashba model for $n$-type semiconductor heterostructures. We find that thermal spin accumulation vanishes for the former but is nonzero for the latter. Onsager’s reciprocity relations are discussed.

In order to implement a temperature gradient $-\bm{\nabla}_{x}T$, which is a statistical force, into the linear-response theory, we use a gravitational electric field $\bm{E}_{\mathrm{g}}$ originally introduced by Luttinger [23]. This field is a mechanical force corresponding to the temperature gradient as $\bm{E}_{\mathrm{g}}\sim-\bm{\nabla}_{x}T/T$. Similar to an electric field, this field is defined by either the gradient of Luttinger’s gravitational potential or the time derivative of a gravitational vector potential [51, 52]. The former is coupled to the Hamiltonian, while the latter is coupled to the heat current and allows us to evaluate the spin–heat-current correlation function as the spin response to (the gradient of) the temperature gradient.

We calculate the linear response of $\hat{s}_{a}$ to the (gravitational) electric-field gradient,

$$\langle\Delta\hat{s}_{a}\rangle^{(1,1)}(\Omega,\bm{Q})=g_{a}^{\phantom{a}ij}(iQ_{i})E_{j}(\Omega,\bm{Q})+Tg_{qa}^{\phantom{qa}ij}(iQ_{i})E_{\mathrm{g}j}(\Omega,\bm{Q}).$$

(1)

Since the (gravitational) electric field is defined as $\bm{E}_{\mathrm{(g)}}(\Omega,\bm{Q})=-(-i\Omega)\bm{A}_{\mathrm{(g)}}(\Omega,\bm{Q})$ using the (gravitational) vector potential $\bm{A}_{\mathrm{(g)}}$, we focus on the first-order terms with respect to the frequency $\Omega$ and wave number $\bm{Q}$ in the spin–charge- (heat-)current correlation function, which is implied by the superscript $(1,1)$. The temperature $T$ in the second term is introduced to compensate the denominator of $\bm{E}_{\mathrm{g}}\sim-\bm{\nabla}_{x}T/T$. Thus, $g_{qa}^{\phantom{qa}ij}$ characterizes the thermal response $\langle\Delta\hat{s}_{a}\rangle^{(1,1)}=g_{qa}^{\phantom{qa}ij}\partial_{x^{i}}(-\partial_{x^{j}}T)$.

In the presence of time-reversal symmetry, only the dissipative Fermi-surface term contributes to spin accumulation as

	$\displaystyle g_{a}^{\phantom{a}ij}=$	$\displaystyle\int d\epsilon[-f^{\prime}(\epsilon)]g_{0a}^{\phantom{0a}ij}(\epsilon),$		(2a)
	$\displaystyle Tg_{qa}^{\phantom{qa}ij}=$	$\displaystyle\int d\epsilon[-f^{\prime}(\epsilon)]Tg_{q0a}^{\phantom{q0a}ij}(\epsilon),$		(2b)

in which $f(\epsilon)=[e^{(\epsilon-\mu)/T}+1]^{-1}$ is the Fermi distribution function, and

	$\displaystyle g_{0a}^{\phantom{0a}ij}(\epsilon)=$	$\displaystyle\frac{i\hbar^{2}q}{8\pi}\int\frac{d^{d}k}{(2\pi)^{d}}\operatorname{tr}[-2\hat{S}_{a}^{\mathrm{AR}}\hat{G}^{\mathrm{R}}(\hat{V}^{\mathrm{RR}i}\hat{G}^{\mathrm{R}}\hat{V}^{\mathrm{RA}j}-\hat{V}^{\mathrm{RA}j}\hat{G}^{\mathrm{A}}\hat{V}^{\mathrm{AA}i})\hat{G}^{\mathrm{A}}$
		$\displaystyle+\hat{S}_{a}^{\mathrm{AA}}\hat{G}^{\mathrm{A}}(\hat{V}^{\mathrm{AA}i}\hat{G}^{\mathrm{A}}\hat{V}^{\mathrm{AA}j}-\hat{V}^{\mathrm{AA}j}\hat{G}^{\mathrm{A}}\hat{V}^{\mathrm{AA}i})\hat{G}^{\mathrm{A}}+\hat{S}_{a}^{\mathrm{RR}}\hat{G}^{\mathrm{R}}(\hat{V}^{\mathrm{RR}i}\hat{G}^{\mathrm{R}}\hat{V}^{\mathrm{RR}j}-\hat{V}^{\mathrm{RR}j}\hat{G}^{\mathrm{R}}\hat{V}^{\mathrm{RR}i})\hat{G}^{\mathrm{R}}],$		(3a)
	$\displaystyle Tg_{q0a}^{\phantom{q0a}ij}(\epsilon)=$	$\displaystyle\frac{i\hbar^{2}}{8\pi}\int\frac{d^{d}k}{(2\pi)^{d}}\operatorname{tr}[-2\hat{S}_{a}^{\mathrm{AR}}\hat{G}^{\mathrm{R}}(\hat{V}^{\mathrm{RR}i}\hat{G}^{\mathrm{R}}\hat{\mathcal{J}}_{q}^{\mathrm{RA}j}-\hat{\mathcal{J}}_{q}^{\mathrm{RA}j}\hat{G}^{\mathrm{A}}\hat{V}^{\mathrm{AA}i})\hat{G}^{\mathrm{A}}$
		$\displaystyle+\hat{S}_{a}^{\mathrm{AA}}\hat{G}^{\mathrm{A}}(\hat{V}^{\mathrm{AA}i}\hat{G}^{\mathrm{A}}\hat{\mathcal{J}}_{q}^{\mathrm{AA}j}-\hat{\mathcal{J}}_{q}^{\mathrm{AA}j}\hat{G}^{\mathrm{A}}\hat{V}^{\mathrm{AA}i})\hat{G}^{\mathrm{A}}+\hat{S}_{a}^{\mathrm{RR}}\hat{G}^{\mathrm{R}}(\hat{V}^{\mathrm{RR}i}\hat{G}^{\mathrm{R}}\hat{\mathcal{J}}_{q}^{\mathrm{RR}j}-\hat{\mathcal{J}}_{q}^{\mathrm{RR}j}\hat{G}^{\mathrm{R}}\hat{V}^{\mathrm{RR}i})\hat{G}^{\mathrm{R}}].$		(3b)

We have omitted the arguments of $\epsilon$ and $\bm{k}$ for simplicity. Equation (3a) was derived previously [49], while Eq. (3b) is a new result. Here, $q$ is the electron charge, and $d$ is the spatial dimension. $\hat{G}^{\alpha}(\epsilon,\bm{k})=[\epsilon-\hat{\Sigma}^{\alpha}(\epsilon,\bm{k})-\hat{\mathcal{H}}(\bm{k})]^{-1}$ ($\alpha=\mathrm{R},\mathrm{A}$) is the renormalized Green’s function with the self-energy $\hat{\Sigma}^{\alpha}(\epsilon,\bm{k})$. $\hat{S}_{a}^{\alpha\beta}(\epsilon,\bm{k})$, $\hat{V}^{\alpha\beta i}(\epsilon,\bm{k})$, and $\hat{\mathcal{J}}_{q}^{\alpha\beta j}(\epsilon,\bm{k})$ are the renormalized vertices of the spin, velocity, and heat current, respectively. See Eqs. (11) and (12) for a specific case. The bare heat current vertex should be defined by the product of the time derivative and velocity [53, 54], rather than that of the Hamiltonian and velocity. We do not specify either the Hamiltonian or the self-energy, and hence Eq. (3) can be applied to any disordered or interacting system as far as the time-reversal symmetry is preserved.

Hereafter, we restrict ourselves to disordered systems without inelastic scattering. Since elastic scattering does not carry energy, the renormalized heat current vertex $\hat{\mathcal{J}}_{q}^{\alpha\beta j}(\epsilon,\bm{k})$ is reduced to $\hat{\mathcal{J}}_{q}^{\alpha\beta j}(\epsilon,\bm{k})=(\epsilon-\mu)\hat{V}^{\alpha\beta j}(\epsilon,\bm{k})$. Then, we obtain $Tg_{q0a}^{\phantom{q0a}ij}(\epsilon)=(\epsilon-\mu)g_{0a}^{\phantom{0a}ij}(\epsilon)/q$, and Eq. (2) becomes the generalized Mott relation, which is reduced to the Mott relation near zero temperature using the Sommerfeld expansion. This is our main result. The only assumptions are the presence of time-reversal symmetry and the absence of inelastic scattering.

First, we consider the three-dimensional isotropic Luttinger model [50] for $p$-type semiconductors, which was studied in the context of the SH effect [55, 56]. The Hamiltonian is expressed as

$$\hat{\mathcal{H}}(\bm{k})=\frac{\hbar^{2}}{2m}\left[\left(\gamma_{1}+\frac{5}{2}\gamma_{2}\right)k^{2}-2\gamma_{2}(\bm{k}\cdot\bm{\Sigma})^{2}\right],$$

(4)

in which $\bm{\Sigma}$ is the $j=3/2$ matrix. We can construct $\Gamma$ matrices satisfying $\{\Gamma_{a},\Gamma_{b}\}=2\delta_{ab}$ ($a=1,\dots,5$) from products of $\bm{\Sigma}$ [56]. Using these $\Gamma$ matrices, this Hamiltonian can be rewritten as

$$\hat{\mathcal{H}}(\bm{k})=\frac{\hbar^{2}k^{2}}{2m}\left[\gamma_{1}+2\gamma_{2}\sum_{a=1}^{5}h^{a}(\bm{k})\Gamma_{a}\right],$$

(5)

in which

	$\displaystyle h^{1}(\bm{k})=$	$\displaystyle-\sqrt{3}\sin\theta\cos\theta\sin\phi,$		(6a)
	$\displaystyle h^{2}(\bm{k})=$	$\displaystyle-\sqrt{3}\sin\theta\cos\theta\cos\phi,$		(6b)
	$\displaystyle h^{3}(\bm{k})=$	$\displaystyle-\frac{\sqrt{3}}{2}\sin^{2}\theta\sin 2\phi,$		(6c)
	$\displaystyle h^{4}(\bm{k})=$	$\displaystyle-\frac{\sqrt{3}}{2}\sin^{2}\theta\cos 2\phi,$		(6d)
	$\displaystyle h^{5}(\bm{k})=$	$\displaystyle-\frac{1}{2}(3\cos^{2}\theta-1).$		(6e)

The eigenvalues are $\epsilon_{\sigma}(\bm{k})=(\gamma_{1}+2\sigma\gamma_{2})\hbar^{2}k^{2}/2m$ with $\sigma=\pm 1$, each of which is doubly degenerate owing to the time-reversal and inversion symmetries. $\sigma=+1$ is called a light hole, while $\sigma=-1$ is called a heavy hole.

We consider $\delta$-function nonmagnetic disorder within the first Born approximation. The bare retarded Green’s function is expressed as

	$\displaystyle\hat{g}^{\mathrm{R}}(\epsilon,\bm{k})=$	$\displaystyle\frac{1}{\epsilon+i\eta-\hat{\mathcal{H}}(\bm{k})}$
	$\displaystyle=$	$\displaystyle\frac{1}{2}[g_{+}^{\mathrm{R}}(\epsilon,\bm{k})+g_{-}^{\mathrm{R}}(\epsilon,\bm{k})]$
		$\displaystyle+\frac{1}{2}[g_{+}^{\mathrm{R}}(\epsilon,\bm{k})-g_{-}^{\mathrm{R}}(\epsilon,\bm{k})]\sum_{a=1}^{5}h^{a}(\bm{k})\Gamma_{a},$		(7)

with $g_{\sigma}^{\mathrm{R}}(\epsilon,\bm{k})=[\epsilon+i\eta-\epsilon_{\sigma}(\bm{k})]^{-1}$. In the first Born approximation, the imaginary part of the self-energy is

$$\hat{\Gamma}(\epsilon)=-\operatorname{Im}\left[n_{\mathrm{i}}v_{\mathrm{i}}^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\hat{g}^{\mathrm{R}}(\epsilon,\bm{k})\right]=\pi n_{\mathrm{i}}v_{\mathrm{i}}^{2}D(\epsilon),$$

(8)

in which

	$\displaystyle D(\epsilon)=$	$\displaystyle\frac{1}{2}\sum_{\sigma}\int_{0}^{\infty}\frac{k^{2}dk}{2\pi^{2}}\frac{-1}{\pi}\operatorname{Im}[g_{\sigma}^{\mathrm{R}}(\epsilon,\bm{k})]$
	$\displaystyle=$	$\displaystyle\frac{m}{4\pi^{2}\hbar^{2}}\sum_{\sigma}\frac{1}{(\gamma_{1}+2\sigma\gamma_{2})^{3/2}}\sqrt{\frac{2m\epsilon}{\hbar^{2}}},$		(9)

is the density of states. Below, we denote $\hat{\Gamma}(\epsilon)$ as $\Gamma(\epsilon)$. The renormalized retarded Green’s function is given as

	$\displaystyle\hat{G}^{\mathrm{R}}(\epsilon,\bm{k})=$	$\displaystyle\frac{1}{\epsilon+i\Gamma(\epsilon)-\hat{\mathcal{H}}(\bm{k})}$
	$\displaystyle=$	$\displaystyle\frac{1}{2}[G_{+}^{\mathrm{R}}(\epsilon,\bm{k})+G_{-}^{\mathrm{R}}(\epsilon,\bm{k})]$
		$\displaystyle+\frac{1}{2}[G_{+}^{\mathrm{R}}(\epsilon,\bm{k})-G_{-}^{\mathrm{R}}(\epsilon,\bm{k})]\sum_{a=1}^{5}h^{a}(\bm{k})\Gamma_{a},$		(10)

with $G_{\sigma}^{\mathrm{R}}(\epsilon,\bm{k})=[\epsilon+i\Gamma(\epsilon)-\epsilon_{\sigma}(\bm{k})]^{-1}$.

We take into account the vertex corrections. First, $\hat{V}^{\mathrm{RR}i}(\epsilon,\bm{k})=\hat{V}^{\mathrm{AA}i}(\epsilon,\bm{k})=\hat{v}^{i}(\bm{k})=\hbar^{-1}\partial_{k_{i}}\hat{\mathcal{H}}(\bm{k})$, since $\delta$-function disorder does not carry a wave number. Second, $\hat{V}^{\mathrm{RA}y}(\epsilon,\bm{k})$ is obtained by solving

	$\displaystyle\hat{V}^{\mathrm{RA}y}(\epsilon,\bm{k})=$	$\displaystyle\hat{v}^{y}(\bm{k})+n_{\mathrm{i}}v_{\mathrm{i}}^{2}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}$
		$\displaystyle\times\hat{G}^{\mathrm{R}}(\epsilon,\bm{k}^{\prime})\hat{V}^{\mathrm{RA}y}(\epsilon,\bm{k}^{\prime})\hat{G}^{\mathrm{A}}(\epsilon,\bm{k}^{\prime}),$		(11)

but results in $\hat{V}^{\mathrm{RA}y}(\epsilon,\bm{k})=\hat{v}^{y}(\bm{k})$ owing to the inversion symmetry [57]. $\hat{S}_{z}^{\mathrm{AR}}(\epsilon,\bm{k})$ is obtained by solving

	$\displaystyle\hat{S}_{z}^{\mathrm{AR}}(\epsilon,\bm{k})=$	$\displaystyle\hat{s}_{z}+n_{\mathrm{i}}v_{\mathrm{i}}^{2}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}$
		$\displaystyle\times\hat{G}^{\mathrm{A}}(\epsilon,\bm{k}^{\prime})\hat{S}_{z}^{\mathrm{AR}}(\epsilon,\bm{k}^{\prime})\hat{G}^{\mathrm{R}}(\epsilon,\bm{k}^{\prime}).$		(12)

We put $\hat{s}_{z}=(\hbar/3)\Sigma_{z}$, because the contribution of the spin $s=1/2$ to the total angular momentum $j=3/2$ is $1/3$ [55]. The solution takes the form of $\hat{S}_{z}^{\mathrm{AR}}(\epsilon,\bm{k})=(\hbar/3)S_{z}^{z}(\epsilon)\Sigma_{z}$. Thus, $g_{0z}^{\phantom{0z}xy}(\epsilon)$ is obtained as

	$\displaystyle g_{0z}^{\phantom{0z}xy}(\epsilon)=$	$\displaystyle\frac{\gamma_{2}\hbar^{5}q}{12\pi m^{2}}\int_{0}^{\infty}\frac{k^{2}dk}{2\pi^{2}}\left(S_{z}^{z}(\epsilon)\{(\gamma_{1}+2\gamma_{2})[(G_{+}^{\mathrm{A}})^{2}G_{-}^{\mathrm{R}}+G_{-}^{\mathrm{A}}(G_{+}^{\mathrm{R}})^{2}]-(\gamma_{1}-2\gamma_{2})[(G_{-}^{\mathrm{A}})^{2}G_{+}^{\mathrm{R}}+G_{+}^{\mathrm{A}}(G_{-}^{\mathrm{R}})^{2}]\right.$
		$\displaystyle+(\gamma_{1}+\gamma_{2})G_{+}^{\mathrm{A}}G_{+}^{\mathrm{R}}(G_{-}^{\mathrm{A}}+G_{-}^{\mathrm{A}})-(\gamma_{1}-5\gamma_{2})G_{-}^{\mathrm{A}}G_{-}^{\mathrm{R}}(G_{+}^{\mathrm{A}}+G_{+}^{\mathrm{A}})\}$
		$\displaystyle\left.-(2\gamma_{1}+3\gamma_{2})[(G_{+}^{\mathrm{A}})^{2}G_{-}^{\mathrm{A}}+(G_{+}^{\mathrm{R}})^{2}G_{-}^{\mathrm{R}}]+(2\gamma_{1}-7\gamma_{2})[G_{+}^{\mathrm{A}}(G_{-}^{\mathrm{A}})^{2}+G_{+}^{\mathrm{R}}(G_{-}^{\mathrm{R}})^{2}]\right)k^{2}$
	$\displaystyle=$	$\displaystyle\frac{\gamma_{2}\hbar^{3}q}{12m}\frac{S_{z}^{z}(\epsilon)}{\Gamma(\epsilon)}\sum_{\sigma}\int_{0}^{\infty}\frac{k^{2}dk}{2\pi^{2}}\frac{-1}{\pi}\operatorname{Im}\left(\left\{3\sigma-\frac{2\gamma_{2}\hbar^{2}k^{2}/m}{\gamma_{2}\hbar^{2}k^{2}/m-i\sigma\Gamma(\epsilon)}+\frac{\gamma_{1}\gamma_{2}(\hbar^{2}k^{2}/m)^{2}}{[\gamma_{2}\hbar^{2}k^{2}/m-i\sigma\Gamma(\epsilon)]^{2}}\right\}G_{\sigma}^{\mathrm{R}}\right)$
		$\displaystyle+\frac{\hbar^{3}q}{12m}\sum_{\sigma}\int_{0}^{\infty}\frac{k^{2}dk}{2\pi^{2}}\frac{1}{\pi}\operatorname{Re}\left\{\frac{\gamma_{2}S_{z}^{z}(\epsilon)k\partial_{k}G_{\sigma}^{\mathrm{R}}}{\gamma_{2}\hbar^{2}k^{2}/m-i\sigma\Gamma(\epsilon)}-[2\gamma_{1}+(5\sigma-2)\gamma_{2}](G_{\sigma}^{\mathrm{R}})^{2}+\frac{2\sigma(\gamma_{1}-\gamma_{2})G_{\sigma}^{\mathrm{R}}}{\gamma_{2}\hbar^{2}k^{2}/m}\right\},$		(13)

with

	$\displaystyle S_{z}^{z}(\epsilon)=$	$\displaystyle\left[1-\frac{n_{\mathrm{i}}v_{\mathrm{i}}^{2}}{10}\int_{0}^{\infty}\frac{k^{2}dk}{2\pi^{2}}(3G_{+}^{\mathrm{A}}G_{+}^{\mathrm{R}}+3G_{-}^{\mathrm{A}}G_{-}^{\mathrm{R}}+2G_{+}^{\mathrm{A}}G_{-}^{\mathrm{R}}+2G_{-}^{\mathrm{A}}G_{+}^{\mathrm{R}})\right]^{-1}$
	$\displaystyle=$	$\displaystyle\left(1-\pi n_{\mathrm{i}}v_{\mathrm{i}}^{2}\sum_{\sigma}\int_{0}^{\infty}\frac{k^{2}dk}{2\pi^{2}}\left\{\frac{3}{10\Gamma(\epsilon)}\frac{-1}{\pi}\operatorname{Im}(G_{\sigma}^{\mathrm{R}})+\frac{1}{5}\frac{1}{\pi}\operatorname{Re}\left[\frac{\sigma G_{\sigma}^{\mathrm{R}}}{\gamma_{2}\hbar^{2}k^{2}/m-i\sigma\Gamma(\epsilon)}\right]\right\}\right)^{-1}.$		(14)

We have omitted the arguments of $\epsilon$ and $\bm{k}$ for simplicity.

In the limit of $n_{\mathrm{i}}v_{\mathrm{i}}^{2}\rightarrow+0$, we find

$$g_{0z}^{\phantom{0z}xy}(\epsilon)=\frac{q}{12\pi^{2}}\frac{\hbar S_{z}^{z}(\epsilon)}{2\Gamma(\epsilon)}\sum_{\sigma}\frac{\gamma_{1}+(3\sigma-2)\gamma_{2}}{(\gamma_{1}+2\sigma\gamma_{2})^{3/2}}\sqrt{\frac{2m\epsilon}{\hbar^{2}}},$$

(15)

with $S_{z}^{z}(\epsilon)=5/2$. Since the density of states (9) is also proportional to $\sqrt{2m\epsilon/\hbar^{2}}$, $g_{0z}^{\phantom{0z}xy}(\epsilon)$ is constant. Hence, the electrical spin accumulation $g_{z}^{\phantom{z}xy}$ does not depend on $\mu$ or $T$, and the thermal one $g_{qz}^{\phantom{qz}xy}$ vanishes. Note that this vanishing response results from the accidental cancellation of the energy dependence. The thermal response can be nonzero for long-range nonmagnetic or spin-orbit-coupled disorder or with higher-order approximations.

Second, we consider the two-dimensional Rashba model for $n$-type semiconductor heterostructures, which was also studied in the context of the SH effect [58]. The Hamiltonian is expressed as

$$\hat{\mathcal{H}}(\bm{k})=\frac{\hbar^{2}k^{2}}{2m}+\hbar\alpha(k_{y}\sigma_{x}-k_{x}\sigma_{y}),$$

(16)

where $\bm{\sigma}$ is the Pauli matrix corresponding to the spin operator $\hat{\bm{s}}=(\hbar/2)\bm{\sigma}$. In the case of $\delta$-function nonmagnetic disorder within the first Born approximation, $g_{0z}^{\phantom{0z}xy}(\epsilon)$ was already obtained as [49]

$$g_{0z}^{\phantom{0z}xy}(\epsilon)=\frac{q\tau_{0}}{8\pi}\begin{cases}1&(\epsilon>0),\\ 1+2\epsilon/m\alpha^{2}&(\epsilon<0).\end{cases}$$

(17)

Here, $\tau_{0}=\hbar/2\Gamma_{0}$ is the relaxation time with $\Gamma_{0}=mn_{\mathrm{i}}v_{\mathrm{i}}^{2}/2\hbar^{2}$. The chemical potential dependences of $g_{(q)z}^{\phantom{(q)z}xy}$ are shown in Fig. 2. Apart from the thermal excitation effect, $g_{qz}^{\phantom{qz}xy}$ is nonzero only below the Rashba crossing, where $g_{0z}^{\phantom{0z}xy}(\epsilon)$ depends on $\epsilon$. The plateau at the lowest temperature in Fig. 2(b) takes the value of $2\pi^{2}/3$ according to the Mott relation.

It has been believed that the spin current can be measured via the dampinglike (DL) spin torque using the ferromagnetic resonance in a ferromagnet–heavy-metal heterostructure [59, 60]. The unit vector of the magnetization $\bm{n}$ obeys the Landau-Lifshitz-Gilbert equation,

$$\dot{\bm{n}}=\gamma\bm{H}_{\mathrm{eff}}\times\bm{n}+\alpha_{\mathrm{G}}\bm{n}\times\dot{\bm{n}}+\bm{\tau}^{\prime},$$

(18)

in which $\gamma\bm{H}_{\mathrm{eff}}$ is the effective magnetic field including the gyromagnetic ratio, $\alpha_{\mathrm{G}}$ is the Gilbert damping, and $\bm{\tau}^{\prime}$ is an additional spin torque. The spin current is assumed to yield the DL spin transfer torque $\bm{\tau}^{\prime}\propto J_{s}\bm{n}\times(\bm{\sigma}\times\bm{n})$ with the direction of the spin polarization of the spin current $\bm{\sigma}$. On the other hand, spin accumulation yields the fieldlike spin torque $\bm{\tau}^{\prime}=\bm{n}\times J\langle\Delta\hat{\bm{s}}\rangle/\hbar$ if we consider the exchange coupling $\hat{\mathcal{H}}^{\prime}=-J\bm{n}\cdot\hat{\bm{s}}/\hbar$ near the interface. The DL spin torque may arise in combination with the Gilbert damping. In addition, since both the inversion and time-reversal symmetries are broken near the interface, the magnetoelectric effect is allowed and contributes to the DL spin torque, as already pointed out in Ref. [6]. Thus, it is difficult to extract the spin current from the DL spin torque.

Onsager’s reciprocity relations readily follow from Eq. (3). Since spin is conjugate to the Zeeman field, namely, $\hat{\mathcal{H}}^{\prime}=-\bm{B}\cdot\hat{\bm{s}}$, the inverse SH and SN effects are described as

	$\displaystyle\langle\Delta\hat{J}^{j}\rangle^{(1,1)}(\Omega,\bm{Q})=$	$\displaystyle g_{a}^{\phantom{a}ij}(iQ_{i})\dot{B}^{a}(\Omega,\bm{Q}),$		(19a)
	$\displaystyle\langle\Delta\hat{J}_{q}^{j}\rangle^{(1,1)}(\Omega,\bm{Q})=$	$\displaystyle Tg_{qa}^{\phantom{qa}ij}(iQ_{i})\dot{B}^{a}(\Omega,\bm{Q}).$		(19b)

Equation (19a) is similar to the previous formulation of the inverse Edelstein effect [61]. The external force $\dot{\bm{B}}$ can be realized by the spin pumping in a ferromagnet–heavy-metal heterostructure as follows: We can read $\bm{B}=J\bm{n}/\hbar$ from the aforementioned exchange coupling near the interface. Assuming a simple precession along the $a$ axis, there is a dc component in $\dot{B}^{a}=J\alpha_{\mathrm{G}}(\bm{n}\times\dot{\bm{n}})^{a}/\hbar$, which is interpreted as the spin current in the spin pumping theory [62, 63, 64]. A dc component of $\dot{\bm{B}}$ also arises from the second-order perturbation with respect to the exchange coupling. Indeed, the previous microscopic calculations revealed that the charge current is induced by $\bm{n}\times\dot{\bm{n}}$ [65, 66, 67] and $\partial_{x^{i}}(\bm{n}\times\dot{\bm{n}})$ [67, 68], which can be attributed to the inverse Edelstein and SH effects, respectively.

Finally, let us comment on the importance of time-reversal symmetry. In the presence of this symmetry, as assumed above, the only process to exhibit spin accumulation is a dissipative Fermi-surface term. That is why we do not suffer from a magnetization correction, which is a thermodynamic quantity. If a system lacks time-reversal symmetry, a nondissipative Fermi-sea term is allowed. Note that the SH effect in such a system is called the magnetic SH effect [69]. As in the Nernst and thermal Hall effects [18, 19, 20, 21, 22] and the gravitomagnetoelectric effect [24, 25], the Kubo formula of thermal spin accumulation diverges towards zero temperature. Such unphysical divergence would be corrected by the spin magnetic octupole moment, which is a future problem.

To summarize, we have studied the spin response to the gradient of a temperature gradient and proved the generalized Mott relation to the electric counterpart assuming the presence of time-reversal symmetry and the absence of inelastic scattering only. Since such a gradient appears at the boundaries, this response describes the thermal spin accumulation. In sharp contrast to previous theories on the SN effect, our theory is free from the ambiguity regarding the definition of the spin current or the magnetization correction. We have found that thermal spin accumulation is absent in the isotropic Luttinger model for $p$-type semiconductors but present in the Rashba model for $n$-type semiconductor heterostructures. Onsager’s reciprocity relations are evident in our theory.

This work was supported by the Japan Society for the Promotion of Science KAKENHI (Grants No. JP21H01816 and No. JP22K03498).