Intrinsic spin Nernst effect of magnons in a noncollinear antiferromagnet

TL;DR

This paper develops a theoretical framework to analyze the intrinsic magnon spin current in noncollinear antiferromagnetic insulators, predicting a measurable spin Nernst effect influenced by symmetry considerations.

Contribution

It introduces a new definition of magnon spin current in noncollinear antiferromagnets and formulates a linear response theory accounting for symmetry effects.

Findings

Magnon spin current can be non-conserved but averages to zero under certain symmetries.

The theory predicts a measurable spin Nernst response in specific materials.

Symmetry considerations critically influence the magnon spin response.

Abstract

We investigate the intrinsic magnon spin current in a noncollinear antiferromagnetic insulator. We introduce a definition of the magnon spin current in a noncollinear antiferromagnet and find that it is in general non-conserved, but for certain symmetries and spin polarizations the averaged effect of non-conserving terms can vanish. We formulate a general linear response theory for magnons in noncollinear antiferromagnets subject to a temperature gradient and analyze the effect of symmetries on the response tensor. We apply this theory to single-layer potassium iron jarosite KFe$_3$(OH)$_6$(SO$_4$)$_2$ and predict a measurable spin current response.

Equations

$$J_{\gamma \lambda }=\frac{2k_B}{V}\sum _{n=1}^N\sum _{\mathbf{k}}({\Omega }_{n,\mathbf{k}}^{j_S}{)}_{\beta }^{\gamma \lambda }{c}_1[g({\epsilon }_{n,\mathbf{k}})]{\nabla }_{\beta }T,$$

$$H_0=\frac{1}{2}\int d\mathbf{r}{\Psi }^{†}(\mathbf{r})\hat{H}\Psi (\mathbf{r}),$$

$${\sigma }_3{H}_{\mathbf{k}}|u_{n,\mathbf{k}}^R⟩={\bar{\epsilon }}_{n,\mathbf{k}}|u_{n,\mathbf{k}}^R⟩,$$

$$⟨u_{n,\mathbf{k}}^L|{\sigma }_3{H}_{\mathbf{k}}={\bar{\epsilon }}_{n,\mathbf{k}}⟨u_{n,\mathbf{k}}^L|,$$

$$H^{\prime }=\frac{1}{4}\int d\mathbf{r}{\Psi }^{†}(\mathbf{r})(\chi \hat{H}+\hat{H}\chi )\Psi (\mathbf{r}).$$

$$\frac{\partial \mathcal{O}(\mathbf{r})}{\partial t}=i[H,\mathcal{O}(\mathbf{r})]=-\nabla \cdot {\mathbf{j}}_O+S_O.$$

$${\Theta }_{\alpha }=⟨{\Theta }_{\alpha }^{[0]}{⟩}_{neq}+⟨{\Theta }_{\alpha }^{[1]}{⟩}_{eq},$$

$$⟨{\Theta }_{\alpha }^{[0]}{⟩}_{\text{neq}}=\underset{\omega \to 0}{\lim }\frac{1}{i\omega }[{\Pi }_{\alpha \beta }(\omega )-{\Pi }_{\alpha \beta }(0)]{\nabla }_{\beta }\chi .$$

$$S_{\beta }^{{\theta }_{\alpha }}=\frac{1}{V}\sum _{n\mathbf{k}}-\left[({\Omega }_{n,\mathbf{k}}^{\theta }{)}_{\beta }^{\alpha }{\bar{\epsilon }}_{n,\mathbf{k}}+(m_{n,\mathbf{k}}^{\theta }{)}_{\beta }^{\alpha }\right]g({\bar{\epsilon }}_{n,\mathbf{k}}).$$

$$({\Omega }_{n,\mathbf{k}}^{\theta }{)}_{\beta }^{\alpha }=\sum _{m(\ne n)}({\sigma }_3{)}_{nn}({\sigma }_3{)}_{mm}\frac{2\text{Im}[({\theta }_{\alpha ,\mathbf{k}}{)}_{nm}(v_{\beta ,\mathbf{k}}{)}_{mn}]}{({\bar{\epsilon }}_{n,\mathbf{k}}-{\bar{\epsilon }}_{m,\mathbf{k}}{)}^2},$$

$$(m_{n,\mathbf{k}}^{\theta }{)}_{\beta }^{\alpha }=\sum _{m(\ne n)}({\sigma }_3{)}_{nn}({\sigma }_3{)}_{mm}\frac{-\text{Im}[({\theta }_{\alpha ,\mathbf{k}}{)}_{nm}(v_{\beta ,\mathbf{k}}{)}_{mn}]}{{\bar{\epsilon }}_{n,\mathbf{k}}-{\bar{\epsilon }}_{m,\mathbf{k}}}$$

$$M_{\beta }^{{\theta }_{\alpha }}=\frac{1}{2V}⟨\int d\mathbf{r}{\Psi }^{†}(\mathbf{r})({\hat{\theta }}_{\alpha }{r}_{\beta }+r_{\beta }{\hat{\theta }}_{\alpha })\Psi (\mathbf{r}){⟩}_{eq}.$$

$$M_{\beta }^{{\theta }_{\alpha }}=-\underset{h_{\alpha }\to 0}{\lim }\frac{1}{V}\frac{\partial \Omega }{\partial ({\partial }_{r_{\beta }}{h}_{\alpha })}$$

$${\tilde{M}}_{\beta }^{{\theta }_{\alpha }}=\frac{\partial (\beta M_{\beta }^{{\theta }_{\alpha }})}{\partial \beta }$$

$${\tilde{M}}_{\beta }^{{\theta }_{\alpha }}=\underset{h_{\alpha }\to 0}{\lim }\underset{q\to 0}{\lim }\frac{-2}{h_{\alpha }^{0}V}\int d\mathbf{r}\delta K(\mathbf{r})\cos (\mathbf{q}\cdot \mathbf{r})$$

$$M_{\beta }^{{\theta }_{\alpha }}=\frac{1}{V}\sum _{n\mathbf{k}}[({\Omega }_{n,\mathbf{k}}^{\theta }{)}_{\beta }^{\alpha }{\int }_0^{{\bar{\epsilon }}_{n\mathbf{k}}}d\eta g(\eta )+(m_{n,\mathbf{k}}^{\theta }{)}_{\beta }^{\alpha }g({\bar{\epsilon }}_{n,\mathbf{k}})].$$

$${\Theta }_{\alpha }=\frac{2k_B}{V}\sum _{n=1}^N\sum _{\mathbf{k}}({\Omega }_{n,\mathbf{k}}^{\theta }{)}_{\beta }^{\alpha }{c}_1[g({\epsilon }_{n,\mathbf{k}})]{\nabla }_{\beta }T.$$

$$H=\sum _{⟨ij⟩}{J}_1{\mathbf{S}}_i\cdot {\mathbf{S}}_j+{\mathbf{D}}_{ij}\cdot ({\mathbf{S}}_i\times {\mathbf{S}}_j)+\sum _{⟨⟨ij⟩⟩}{J}_2{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$${\alpha }_{\lambda \beta }^{\gamma }=\det (R)R_{\gamma {\gamma }^{\prime }}{R}_{\lambda {\lambda }^{\prime }}{R}_{\beta {\beta }^{\prime }}{\alpha }_{{\lambda }^{\prime }{\beta }^{\prime }}^{{\gamma }^{\prime }},$$

$$\displaystyle[\hat{\alpha}^{x},\hat{\alpha}^{y},\hat{\alpha}^{z}]=\left[\left(\begin{array}[]{cc}-\alpha_{1}&0\\ 0&\alpha_{1}\end{array}\right),\right.\left(\begin{array}[]{cc}0&\alpha_{1}\\ \alpha_{1}&0\end{array}\right),\left.\left(\begin{array}[]{cc}0&-\alpha_{2}\\ \alpha_{2}&0\end{array}\right)\right].$$

$$[{\Psi }_i(\mathbf{r}),{\Psi }_j^{†}({\mathbf{r}}^{\prime })]=({\sigma }_3{)}_{ij}{\delta }_{\mathbf{r},{\mathbf{r}}^{\prime }},[{\Psi }_i(\mathbf{r}),{\Psi }_j({\mathbf{r}}^{\prime })]=i({\sigma }_2{)}_{ij}{\delta }_{\mathbf{r},{\mathbf{r}}^{\prime }},[{\Psi }_i^{†}(\mathbf{r}),{\Psi }_j^{†}({\mathbf{r}}^{\prime })]=-i({\sigma }_2{)}_{ij}{\delta }_{\mathbf{r},{\mathbf{r}}^{\prime }}$$

$$\frac{\partial \mathcal{O}(\mathbf{r})}{\partial t}=$$

$$=$$

$$\frac{\partial \mathcal{O}(\mathbf{r})}{\partial t}=$$

$$=$$

$${\Theta }_{\alpha }=L_{\beta }^{{\theta }_{\alpha }}{\nabla }_{\beta }\chi =(S_{\beta }^{{\theta }_{\alpha }}+M_{\beta }^{{\theta }_{\alpha }}){\nabla }_{\beta }\chi ,$$

$$E=⟨H{⟩}_{eq}=\frac{1}{2}\sum _{\mathbf{k},n=1}^{2N}({\sigma }_3{)}_{nn}g({\bar{\epsilon }}_{n,\mathbf{k}})⟨{\psi }_{n,\mathbf{k}}|\hat{H}|{\psi }_{n,\mathbf{k}}⟩.$$

$$d\Omega =-SdT-⟨{\Theta }_{\alpha }^{[0]}⟩d{h}_{\alpha }-M_{\beta }^{{\theta }_{\alpha }}d({\partial }_{r_{\beta }}{h}_{\alpha }),$$

$${\left(\frac{\partial M_{\beta }^{{\theta }_{\alpha }}}{\partial T}\right)}_{h_{\alpha },{\partial }_{r_{\beta }}{h}_{\alpha }}={\left[\frac{\partial S}{\partial ({\partial }_{r_{\beta }}{h}_{\alpha })}\right]}_{T,h_{\alpha }}.$$

$${\tilde{M}}_{\beta }^{{\theta }_{\alpha }}=-\frac{1}{V}\frac{\partial K}{\partial ({\partial }_{r_{\beta }}{h}_{\alpha })}$$

$$M_{\beta }^{{\theta }_{\alpha }}={\tilde{M}}_{\beta }^{{\theta }_{\alpha }}+T\frac{\partial M_{\beta }^{{\theta }_{\alpha }}}{\partial T}$$

Full text

In recent years, advances in research on topological properties of electron systems [1] have encouraged explorations of manifestations of topology in many other systems, e.g., magnonic [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], acoustic [14, 15], photonic [16], etc.. The possibility of coupling between various degrees of freedom has led to new visions for spintronics [17, 18], and resulted in new subfields such as spin caloritronics [19], in which spin carriers are manipulated by exciting heat flows. The study of spin currents is fundamental for the field of spintronics, and the spins carried by magnons possess certain advantages over electrons, e.g., low dissipation. At the same time, magnons exhibit rich and fascinating physics associated with the topology of magnonic bands, e.g., the thermal Hall effect has been observed in collinear ferromagent Lu2V2O7 [20]. The spin Nernst effect [21, 22], akin to the spin Hall effect [23], can also be realized in magnon systems [4, 24, 25, 26, 27, 28, 29].

Many spintronics concepts also apply to antiferromagnets [30]. In particular, collinear antiferromagnets can exhibit the spin Seebeck effect [31], spin pumping [32], spin-orbit torque [33], spin Nernst effect [24, 25, 26, 27, 29], etc. Noncollinear antiferromagnets (NAFMs) have attracted considerable attention in recent years, as such systems support nontrivial band structure topology. The anomalous Hall effect [34] and spin Hall effect [35, 36] have been realized in Mn3X (X=Ge,Sn,Ga,Ir,Rh, and Pt) systems, where electrons act as charge or spin carriers. Furthermore, the thermal Hall effect, mediated by magnons, is also identified in NAFM insulators [29, 37, 38, 39, 40]. Nevertheless, the magnon-mediated spin transport in NAFMs [41, 42, 43] has not yet been well explored, especially in the context of the topology of magnon bands. As a new class of materials, NAFMs feature rich magnetic point group symmetries, chirality, and easily tunable properties (e.g., by magnetic field). As a result, studies of spin currents in NAFMs can bring new vitality to spintronics, especially in the context of spin caloritronics. In contrast to the unique spin polarization of a magnon current in the collinear system, the spin current in NAFMs can be arbitrarily polarized, which allows a better control of the spin current. NAFMs typically possess different ground states [44, 45, 46] depending on the ambient environment, e.g., external field, temperature, substrates, and one can envisage using the spin current as a probe of the ground state. Meanwhile, many NAFM materials can also hold exotic quantum effects [47]. Studies of spin currents in such systems can provide a new venue for probing these materials [48]. Motivated by these interesting possibilities, we initiate a discussion on the magnon-mediated spin current physics in noncollinear antiferromagnets therein and hope to stimulate subsequent research on, e.g., spin transport in topological magnon insulators [49], optical generation of magnon-mediated spin currents [50, 51], and many others, as has been discussed above.

In this paper, we formulate a linear response theory of magnon-mediated spin transport induced by temperature gradients in a noncollinear antiferromagnet, concentrating on the intrinsic contribution not reliant on magnon lifetime. The difficulty in considering a NAFM is similar to a typical spin Hall system in which spin is not conserved [52]. Magnons driven by temperature gradients require accounting for the effects associated with the orbital magnetization [53, 24]. We start by discussing the definition of spin current in particle-hole space by following Refs. [52, 24], where spin non-conservation is signaled by a source term. Next, we develop a linear response theory to temperature gradients for a general observable, i.e., the source term (torque) or spin current, and discuss the symmetry constraints. One of our main results is the expression for the intrinsic spin Nernst response in noncollinear antiferromagnetic insulators,

[TABLE]

where $J_{\gamma\lambda}$ is the spin current with polarization $\gamma$, $(\Omega_{n,\mathbf{k}}^{j_{S}})^{\gamma\lambda}_{\beta}$ is the spin Berry curvature of magnons, and $c_{1}(x)=(1+x)\ln(1+x)-x\ln(x)$ is an auxiliary function stemming from the Bose-Einstein statistics of magnons. We apply our theory to the kagome antiferromagnet KFe3(OH)6(SO4)2 (see Fig. 1) and show that the in-plane Dzyaloshinskii-Moriya interaction (DMI) leads to a measurable spin Nernst response. Our study opens a way for future studies of fascinating physics related to spin flows in noncollinear antiferromagnets, e.g., in the context of different magnetic orders and material realizations.

Spin Nernst response. We consider a general antiferromagnet with noncollinear ordering. To capture its magnonic excitations at low temperatures, we perform the Holstein-Primakoff transformation in the limit of large $S$, which leads us to a general Hamiltonian

[TABLE]

where $\Psi(\mathbf{r})=(b_{1}(\mathbf{r})\cdots b_{N}(\mathbf{r}),b^{\dagger}(\mathbf{r})\cdots b_{N}^{\dagger}(\mathbf{r}))^{T}$ is the bosonic field, and $N$ the number of atoms in each unit cell. The particle-hole space representation is necessary to describe the anomalous coupling between magnons in an antiferromagnet.

Due to Bose-Einstein statistics, the eigenvalue problem has to be solved for the matrix $\sigma_{3}H_{\mathbf{k}}$ [54], where here and in what follows we use Pauli matrices $\sigma_{i}$ acting in the particle-hole space. Here $H_{\mathbf{k}}$ is the Hamiltonian matrix in momentum space, which can be diagonalized by a paraunitary matrix $T_{\mathbf{k}}$, i.e., $T_{\mathbf{k}}^{\dagger}H_{\mathbf{k}}T_{\mathbf{k}}=\mathcal{E}_{\mathbf{k}}$, where $\mathcal{E}_{\mathbf{k}}=\text{Diag}(\varepsilon_{1,\mathbf{k}}\cdots\varepsilon_{N,\mathbf{k}},\varepsilon_{1,-\mathbf{k}}\cdots\varepsilon_{N,-\mathbf{k}})$ is the matrix describing eigenvalues, and $T_{\mathbf{k}}$ satisfies $T^{\dagger}_{\mathbf{k}}\sigma_{3}T_{\mathbf{k}}=\sigma_{3}$. We now write a general theory applicable to bosonic systems where the Bloch wave function corresponding to the band dispersion $\varepsilon_{n,\mathbf{k}}$ is given by $|\psi_{n,\mathbf{k}}\rangle=e^{i\mathbf{k}\cdot\mathbf{r}}|u_{n,\mathbf{k}}\rangle$. We can then introduce a notation [55]

[TABLE]

where in terms of the magnonic Hamiltonian $|u_{n,\mathbf{k}}^{R}\rangle=T_{n,\mathbf{k}}$ and $\langle u^{L}_{n,\mathbf{k}}|=T_{n,\mathbf{k}}^{\dagger}\sigma_{3}$ are the right and left eigenstates of the pseudo-Hermitian Hamiltonian, and $\bar{\varepsilon}_{n,\mathbf{k}}=(\sigma_{3}\mathcal{E}_{\mathbf{k}})_{nn}$. Hereafter, we will only refer to the right eigenstates $|u^{R}_{n,\mathbf{k}}\rangle=|u_{n,\mathbf{k}}\rangle$. The normalization relation reads $\langle u_{n,\mathbf{k}}|\sigma_{3}|u_{m,\mathbf{k}}\rangle=(\sigma_{3})_{nm}$. Moreover, the Hamiltonian (2) possesses particle-hole symmetry (PHS) so that the Hamiltonian obeys $\sigma_{1}H_{\mathbf{k}}\sigma_{1}=H_{-\mathbf{k}}^{\ast}$, which leads to relations $\bar{\varepsilon}_{n+N,\mathbf{k}}=-\bar{\varepsilon}_{n,-\mathbf{k}}$ and $|u_{n,\mathbf{k}}\rangle=e^{i\phi_{n}}\sigma_{1}|u_{n+N,-\mathbf{k}}\rangle^{\ast}$, where $\phi_{n}$ is a redundant phase factor.

Because the temperature gradient is a statistical force and doesn’t directly enter the Hamiltonian, we introduce a perturbation corresponding to a pseudo-gravitational potential, $\chi(\mathbf{r})$, to account for the temperature gradient [56, 53, 24],

[TABLE]

With the perturbation, the total Hamiltonian is amended to $H=\frac{1}{2}\int d\mathbf{r}\tilde{\Psi}^{\dagger}(\mathbf{r})\hat{H}\tilde{\Psi}(\mathbf{r})$, where $\tilde{\Psi}(\mathbf{r})=(1+\mathbf{r}\cdot\bm{\nabla}\chi/2)\Psi(\mathbf{r})$. To linear order, the system will respond to a temperature gradient in the same way as to a perturbation with $\chi(\mathbf{r})=-T(\mathbf{r})/T$. We now introduce an arbitrary matrix $\hat{O}$ and a local observable $\mathcal{O}(\mathbf{r})=\frac{1}{2}\Psi^{\dagger}(\mathbf{r})\hat{O}\Psi(\mathbf{r})$. In what follows, we will mostly consider $\hat{O}=\hat{S}^{\alpha}$, which corresponds to the magnon spin density operator given by $\hat{S}^{\alpha}=-\sigma_{0}\otimes\text{Diag}(\left<S^{\alpha}_{1}\right>/S_{1},\cdots,\left<S^{\alpha}_{N}\right>/S_{N})$, where $\alpha=x,y,z$, $\sigma_{0}$ describes the particle-hole space, and averages of spins within a unit cell have been taken. The time evolution of this operator can be obtained from the Heisenberg equation applied to the total Hamiltonian (details in Supplemental Material) [24]

[TABLE]

Here $\mathbf{j}_{O}=\tilde{\Psi}^{\dagger}(\mathbf{r})\hat{\mathbf{j}}_{O}\tilde{\Psi}(\mathbf{r})$ and $S_{O}=\tilde{\Psi}^{\dagger}(\mathbf{r})\hat{S}_{O}\tilde{\Psi}(\mathbf{r})$ correspond to the local current and source densities, respectively, with $\hat{\mathbf{j}}_{O}=\frac{1}{4}(\hat{\mathbf{v}}\sigma_{3}\hat{O}+\hat{O}\sigma_{3}\hat{\mathbf{v}})$, $\hat{S}_{O}=-\frac{i}{2}(\hat{O}\sigma_{3}\hat{H}-\hat{H}\sigma_{3}\hat{O})$, and $\hat{\mathbf{v}}=i[\hat{H},\mathbf{r}]$. To linear order in the temperature gradient, the above densities are explicitly decomposed as $\rho_{\theta}=\rho_{\theta}^{[0]}+\rho_{\theta}^{[1]}$, with $\rho_{\theta}^{[0]}=\Psi^{\dagger}(\mathbf{r})\hat{\theta}\Psi(\mathbf{r})$, $\rho_{\theta}^{[1]}=\frac{1}{2}\Psi^{\dagger}(\mathbf{r})(\hat{\theta}r_{\beta}+r_{\beta}\hat{\theta})\Psi(\mathbf{r})\nabla_{\beta}\chi$, where for $\theta$ one needs to substitute either $\mathbf{j}_{O}$ or $S_{O}$. We will use a four-vector convention in which $\theta_{0}=S_{O}$ and $\bm{\theta}=\mathbf{j}_{O}$. The non-vanishing source term indicates the non-conservation of the observable, for instance, when $\mathcal{O}(\mathbf{r})$ corresponds to spin density, the source term represents torque density. We note in passing that the source term dipole $\mathbf{P}_{O}$ can be defined as $S_{O}=-\bm{\nabla}\cdot\mathbf{P}_{O}$ for vanishing total source $\frac{1}{V}\int d\mathbf{r}S_{O}=0$, where $V$ is the volume and $\mathbf{P}_{O}=\mathbf{r}S_{O}$. Thus, a conserved current can be defined as $\bm{\mathcal{J}}_{O}=\mathbf{j}_{O}+\mathbf{P}_{O}$ to restore the continuity equation [52]. The current term $\mathbf{j}_{O}$ coincides with the conventional definition in the literature of the spin Hall effect [23]. In general, based on Eq. (5) we can interpret $\mathbf{j}_{O}$ as a spin current and $S_{O}$ as the torque. In our discussion below, we concentrate on the spin current term.

We consider spatially averaged quantities $\Theta_{\alpha}=\Theta_{\alpha}^{[0]}+\Theta_{\alpha}^{[1]}$ with $\Theta_{\alpha}^{[0,1]}=\frac{1}{V}\int d\mathbf{r}\rho^{[0,1]}_{\theta_{\alpha}}(\mathbf{r})$. The thermal response to linear order in the temperature gradient reads

[TABLE]

where on the right hand side the first term is evaluated with respect to nonequilibrium states from the Kubo linear response calculation, while the second term corresponds to orbital magnetization in the system and is evaluated with respect to the equilibrium state. In total, we can express the linear response as $\Theta_{\alpha}=(S^{\theta_{\alpha}}_{\beta}+M^{\theta_{\alpha}}_{\beta})\nabla_{\beta}\chi$, where $S^{\theta_{\alpha}}_{\beta}$ and $M^{\theta_{\alpha}}_{\beta}$ correspond to the first and second terms in Eq. (6).

In the spirit of the Kubo response calculation [4, 24], the nonequilibrium part can be described by

[TABLE]

Here $\Pi_{\alpha\gamma}(i\omega_{m})=-\int_{0}^{\beta}d\tau e^{i\omega_{m}\tau}\langle T_{\tau}\Theta^{[0]}_{\alpha}(\tau)J_{\gamma}^{q}(0)\rangle$, with $\beta=1/(k_{B}T)$, where $\omega_{m}$ is the bosonic Matsubara frequency. $\mathbf{J}^{q}$ is the averaged heat current operator defined as $\mathbf{J}^{q}=\frac{1}{V}\int d\mathbf{r}\mathbf{j}^{q}(\mathbf{r})$, where the heat current density $\mathbf{j}^{q}=\frac{1}{4}\Psi^{\dagger}(\mathbf{r})(\hat{H}\sigma_{3}\hat{\mathbf{v}}+\hat{\mathbf{v}}\sigma_{3}\hat{H})\Psi(\mathbf{r})$. This heat current expression can be inferred from the energy conservation equation $\dot{\rho}_{E}+\bm{\nabla}\cdot\mathbf{j}^{q}=0$, where $\rho_{E}$ is the energy density of the system. After performing the linear response calculation, the intrinsic nonequilibrium coefficient reads

[TABLE]

Here

[TABLE]

where $(\dots)_{nm}=\langle u_{n,\mathbf{k}}|\dots|u_{m,\mathbf{k}}\rangle$ and $g(\bar{\varepsilon}_{n\mathbf{k}})$ is the Bose-Einstein distribution. Here $(\Omega_{n,\mathbf{k}}^{\theta})^{\alpha}_{\beta}$ is the generalized Berry curvature calculated for operator $\hat{\theta}_{\alpha}$. This Berry curvature respects the sum rule $\sum_{n=1}^{2N}(\Omega_{n,\mathbf{k}}^{\theta})^{\alpha}_{\beta}=0$, and PHS results in the relation $(\Omega_{n,\mathbf{k}}^{\theta})^{\alpha}_{\beta}=(\Omega_{n+N,-\mathbf{k}}^{\theta})^{\alpha}_{\beta}$ ($1\leq n\leq N$).

The contribution corresponding to $\rho_{\theta_{\alpha}}^{[1]}$ is expressed as

[TABLE]

To calculate this term, we can identify a thermodynamic expression for $M^{\theta_{\alpha}}_{\beta}$ by following Refs. [57, 58, 59, 60, 61]. We introduce a perturbation coupled with a four-component fictitious field $h_{\alpha}(\mathbf{r})$, i.e., $\hat{H}_{1}=\hat{H}_{0}-[\hat{\theta}_{\alpha}h_{\alpha}(\mathbf{r})+h_{\alpha}(\mathbf{r})\hat{\theta}_{\alpha}]$. If the field varies very slowly on the scale of the lattice constant, we can identify a thermodynamic expression

[TABLE]

where $\Omega$ is the thermodynamic grand potential of the system. If we regard the local fictitious field and its gradient as independent variables, we can assert a Maxwell relation $\left(\partial M^{\theta_{\alpha}}_{\beta}/\partial T\right)_{h_{\alpha},\partial_{r_{\beta}}h_{\alpha}}=\left[\partial S/\partial(\partial_{r_{\beta}}h_{\alpha})\right]_{T,h_{\alpha}}$, where $S$ is the entropy. Taking both Eq. (11) and the Maxwell relation into account, we are led to

[TABLE]

with $\tilde{M}^{\theta_{\alpha}}_{\beta}=-\frac{1}{V}\frac{\partial K}{\partial(\partial_{r_{\beta}}h_{\alpha})}$ being an auxiliary quantity and $K=\Omega+TS$. We assume that the fictitious field takes the form $h_{\alpha}(\bm{r})=(h^{0}_{\alpha}/q)\sin(\mathbf{q}\cdot\mathbf{r})$, with $\mathbf{q}=q\hat{\mathbf{e}}_{\beta}$ ($\beta=x,y,z$ in three dimensions and $\beta=x,y$ in two dimensions). The auxiliary quantity is calculated by picking up the appropriate Fourier component

[TABLE]

where $\delta K(\mathbf{r})$ is the variation due to the fictitious field, which can be obtained from perturbation theory [58]. Combining Eq. (13) and (12), we obtain (see details in Supplemental Material)

[TABLE]

By combining the nonequilibrium part in Eq. (8) with Eq. (Intrinsic spin Nernst effect of magnons in a noncollinear antiferromagnet) and cancelling the orbital part (corresponding to a bound current), we obtain the thermal response formula which constitutes the main result of this paper:

[TABLE]

Note that we express our result using particle bands ($n\leq N$) by utilizing PHS.

It is useful to identify the symmetry constraints leading to a vanishing source term response. In general, for the averaged torque density this can happen for only some of the torque components. However, for an inversion symmetric system, i.e., $H_{\mathbf{k}}=H_{-\mathbf{k}}$, the Berry curvature of the torque term satisfies $(\Omega_{n,\mathbf{k}}^{S_{O}})_{\beta}=-(\Omega_{n,-\mathbf{k}}^{S_{O}})_{\beta}$. Together with the relation $\varepsilon_{n,\mathbf{k}}=\varepsilon_{n,-\mathbf{k}}$, this results in the vanishing of all torque components in Eq. (15).

Spin Nernst effect in kagome antiferromagnet. We use the result in Eq. (15) to calculate the spin Nernst response tensor in a noncollinear kagome antiferromagnet in Eq. (1) where the spin Berry curvature is calculated with respect to operator $\hat{j}_{\gamma,\lambda}=\frac{1}{4}(\hat{v}_{\lambda}\sigma_{3}\hat{S}^{\gamma}+\hat{S}^{\gamma}\sigma_{3}\hat{v}_{\lambda})$ corresponding to the spin current. We can immediately identify that the spin Berry curvature in Eq. (1) is even under the time reversal transformation. As a result, the spin Nernst conductivity is also even under the time reversal transformation, and this result will be used in the symmetry analysis below. Furthermore, in a kagome antiferromagnet, due to the presence of inversion symmetry, the averaged torque density (source term) vanishes. We consider the Hamiltonian

[TABLE]

where the first and third terms represent nearest and second-nearest neighbor Heisenberg exchange, and the second term represents nearest neighbor Dzyaloshinskii-Moriya interaction (DMI) with both in-plane and out-of-plane DMI vectors, as shown in Fig. 1. The DMI vector can be expressed as $\mathbf{D}_{ij}=D_{p}\hat{n}_{ij}+D_{z}\hat{z}$, where $D_{p}$ and $D_{z}$ correspond to the in-plane and out-of-plane DMI strength, and $\hat{n}_{ij}$ is an in-plane unit vector corresponding to the direction of the in-plane DMI. The in-plane DMI can only arise when $\mathcal{M}_{z}$ symmetry is broken [34], i.e., time-reversal followed by mirror symmetry with respect to the kagome plane is not a symmetry in such a case. This introduces a small out-of-plane canting angle $\eta$ to spin order with magnitude $\eta=\frac{1}{2}\tan^{-1}(\frac{-2D_{p}}{\sqrt{3}(J_{1}+J_{2})-D_{z}})$ [29]. Here we consider the $\mathbf{q}=0$ phase with spin order as shown in Fig. 1. The magnetic moments orient according to $\left<\mathbf{S}_{i}\right>=S(\cos\eta\cos\phi_{i},\cos\eta\sin\phi_{i},\sin\eta)$, where $\phi_{i}$ is the angle formed by the in-plane projection of moment with the $x$ axis. Specifically, $\phi_{A}=\pi/2$, $\phi_{B}=7\pi/6$, and $\phi_{C}=-\pi/6$. For the spin Nernst response, we identify $\hat{O}$ discussed above as the spin operator in the magnon basis $\Psi(\mathbf{r})=[b_{A}(\mathbf{r}),b_{B}(\mathbf{r}),b_{C}(\mathbf{r}),b^{\dagger}_{A}(\mathbf{r}),b^{\dagger}_{B}(\mathbf{r}),b^{\dagger}_{C}(\mathbf{r})]^{T}$, i.e., $\hat{S}^{\alpha}=-\sigma_{0}\otimes\text{Diag}(\left<S^{\alpha}_{A}\right>,\left<S^{\alpha}_{B}\right>,\left<S^{\alpha}_{C}\right>)/S$. The spin conductivity tensor of a spin-polarized current in a noncollinear antiferromagnet [36, 62] is restricted to a certain form by the magnetic space group of the system. Given that the intrinsic spin Nernst tensor in relation $J_{\gamma\lambda}=\alpha^{\gamma}_{\lambda\beta}\nabla_{\beta}T$ is even under the time reversal transformation, the symmetry constraints become

[TABLE]

where the matrix $R$ represents a symmetry element $\mathcal{R}$ (in Cartesian coordinates) entering the antiunitary symmetry $\mathcal{R}\mathcal{T}$ or unitary symmetry $\mathcal{R}$ of the system (see Supplemental Material). As an example, we focus on a system with two symmetries: mirror reflection with respect to the $y-z$ plane combined with time reversal, $\mathcal{M}_{x}\mathcal{T}$, and threefold rotation about the $z$ axis, $\mathcal{C}_{3z}$. The shape of the spin Nernst tensor corresponding to the constraints in Eq. (17) becomes

[TABLE]

Here, the $\mathcal{M}_{x}\mathcal{T}$ symmetry can be replaced by $\mathcal{C}_{2x}\mathcal{T}$, twofold rotation about the $x$ axis and time-reversal, which will lead us to the same result. We note that our results are consistent with the spin Hall response tensors in Mn3X (X= Rh, Ir and Pt) [35].

We apply our theory to a single layer of potassium iron jarosite, KFe3(OH)6(SO4)2, for which the material parameters are $J_{1}=3.18\text{meV}$, $J_{2}=0.11\text{meV}$, $|D_{p}|/J_{1}=0.062$, $D_{z}/J_{1}=-0.062$ [29, 63]. We note, however, that the magnon dispersion in this material can also be explained by $J_{2}=0$, in which case the flat band is broadened by fluctuations [64]. The numerically obtained form of the spin Nernst conductivities agrees with Eq. (24). In Fig. 2, we plot the magnon bands and the spin Berry curvature for the $y$ polarization of the spin. The spin Berry curvature is peaked at avoided crossings, which give the largest contribution to the spin Nernst effect. The integral of the ordinary Berry curvature gives the Chern numbers $-3$, $1$, and $2$, from the bottom to the top bands in Fig. 2. In Fig. 3, we show the spin Nernst response coefficients as a function of temperature for the $y$ and $z$ spin polarizations. The spin Nernst response sharply increases at temperatures sufficient to excite magnons in the Brillouin zone where the spin Berry curvature is large. The $z$ direction polarized spin current is two orders of magnitude smaller than the current with in-plane spin polarization, which is due to the fact that the canting angle is fairly small, $\eta=1.9^{\circ}$ [29]. By applying magnetic field, the canting angle and the spin Nernst response with the $z$ polarization direction can be substantially increased. The predicted spin currents should be easily detectable in three dimensional structures as a temperature gradient of 20 K/mm should result in a spin current of the order of $10^{-11}$ J/m2 according to Fig. 3, where $\alpha^{3D}=\alpha/c$, with $c$ being the interlayer distance. Finally, we note that the spin Nernst effect reported in Ref. [41] differs from the intrinsic effect reported here as the former has the symmetry of the extrinsic effect.

Conclusions. We have developed a theory of magnon-mediated intrinsic spin currents in insulating noncollinear antiferromagnets and applied this theory to potassium iron jarosite KFe3(OH)6(SO4)2. Our results are applicable to two- and three-dimensional systems, promising to reveal fascinating physics in other layered quasi-2D antiferromagnets, e.g., silver iron jarosite AgFe3(OH)6(SO$4$)2 [65], chromium jarosite KCr3(OH)6(SO$4$)2 [66], vesignieite BaCu3V2O8(OH)2 [67], and 3D pyrochlore antiferromagnets LiGaGr4O8 and LiInGr4O8 [68, 69]. Besides exploring material candidates, one can also study the effect of magnetic order on the spin Nernst effect, e.g., in kagome antiferromagnets other possible spin chiralities exist [46, 40]. Recently proposed antiferromagnetic skyrmions with noncollinear magnetic order [70] can also be explored using our theory.

We gratefully acknowledge stimulating discussions with Kirill Belashchenko, Ding-Fu Shao and Liang Dong. This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0014189.

I Time evolution of a local observable

To derive the time evolution equation for a local observable $\mathcal{O}(\mathbf{r})=\frac{1}{2}\Psi^{\dagger}(\mathbf{r})\hat{O}\Psi(\mathbf{r})$, we first prepare a basic knowledge of the Hamiltonian operator and commutators in particle-hole space by following Ref. [24, 53]. The total Hamiltonian can be generally expressed as $H=\frac{1}{2}\int d\mathbf{r}\tilde{\Psi}^{\dagger}(\mathbf{r})\hat{H}\tilde{\Psi}(\mathbf{r})$ with $\hat{H}=\sum_{\bm{\delta}}H_{\bm{\delta}}e^{i\hat{\mathbf{p}}\cdot\bm{\delta}}$, in which $e^{i\hat{\mathbf{p}}\cdot\bm{\delta}}$ is the translation operator that satisfies $e^{i\hat{\mathbf{p}}\cdot\bm{\delta}}f(\mathbf{r})=f(\mathbf{r}+\bm{\delta})$. Here $\bm{\delta}$ is the vector shift between unit cells, $\tilde{\Psi}(\mathbf{r})=(1+\mathbf{r}\cdot\bm{\nabla}\chi/2)\Psi(\mathbf{r})$. Based on the basic commutators between bosons $[a_{i}(\mathbf{r}),a^{\dagger}_{j}(\mathbf{r^{\prime}})]=\delta_{ij}\delta_{\mathbf{r},\mathbf{r}^{\prime}}$, $[a_{i}(\mathbf{r}),a_{j}(\mathbf{r^{\prime}})]=0$, we can construct commutators in the particle-hole basis

[TABLE]

where $\sigma_{i}$ ($i=1,2,3$) are Pauli matrices acting in particle-hole space. Now we use the above Hamiltonian and commutators to perform a local observable time evolution calculation in two steps. First, we work out the Heisenberg equation commutation as follow,

[TABLE]

Here we used the simplified notation $\xi(\mathbf{r})=1+\mathbf{r}\cdot\bm{\nabla}\chi/2$. We also took advantage of particle-hole symmetry, i.e., $\Psi_{n}(\mathbf{r})=(\sigma_{1})_{nl}\Psi^{\dagger}_{l}(\mathbf{r})$ and $\sigma_{1}\hat{O}\sigma_{1}=\hat{O}$, where the second relation results from the first one. Next, we reduce the above result to a continuous expression by properly sending the shift vector to an infinitely small value.

[TABLE]

Here we used the notation $\hat{v}_{\bm{\delta}}=i\bm{\delta}H_{\bm{\delta}}e^{i\hat{\mathbf{p}}\cdot\bm{\delta}}$ and $\hat{\mathbf{v}}=i\sum_{\bm{\delta}}\bm{\delta}H_{\bm{\delta}}e^{i\hat{\mathbf{p}}\cdot\bm{\delta}}=i[\hat{H},\mathbf{r}]$. In the last line, we take the limit $\bm{\delta}\rightarrow 0$ to obtain the continuous expression. We can easily read out the current and source term as discussed in the main text from the final result.

II Linear response theory

We provide a fully quantum mechanical derivation in this section. As shown in the main text, the linear thermal response for a given observable can be expressed as

[TABLE]

where $S^{\theta_{\alpha}}_{\beta}$ and $M^{\theta_{\alpha}}_{\beta}$ are the Kubo response and dipole moment contribution, respectively. In the following, we will first calculate the dipole moment part from a thermodynamic point of view and then combine it with the Kubo part to arrive at the final response formula. All calculations will be performed in the full particle-hole space, and we will express the final result in terms of particle space in the end.

II.1 Dipole moment contribution

In the main text, we have shown the relation $M^{\theta_{\alpha}}_{\beta}=-\lim_{h_{\alpha}\rightarrow 0}\frac{1}{V}\frac{\partial\Omega}{\partial(\partial_{r_{\beta}}h_{\alpha})}$. The thermodynamic definition of grand potential reads $\Omega=E-TS-\mu N$, where $S$, $\mu$, $N$ are the entropy, chemical potential, particle number of the system, and $E$ is the energy which reads

[TABLE]

Here we use the relation $\langle\Gamma^{\dagger}_{n,\mathbf{k}}\Gamma_{m,\mathbf{k}}\rangle=(\sigma_{3})_{nn}g(\bar{\varepsilon}_{n,\mathbf{k}})$ with $\Gamma_{m,\mathbf{k}}=\sum_{l}(T_{\mathbf{k}})_{ml}\Psi_{\mathbf{k},l}$. Below we will assume the chemical potential to be zero. If we regard the local fictitious field and its gradient as independent variables, the variation of the grand potential can be identified as

[TABLE]

from which we can identify the Maxwell relation

[TABLE]

To get rid of calculations involving entropy $S$, we first introduce an auxiliary quantity

[TABLE]

with $K=\Omega+TS=E$ ($\mu=0$). By utilizing Eq. (31), we obtain

[TABLE]

and hence the dipole moment contribution can be calculated as

[TABLE]

If we regard the fictitious field term as a perturbation, the variation of $K$ to linear order reads

[TABLE]

where $|\psi_{n\mathbf{k}}\rangle=\frac{e^{i\mathbf{k}\cdot\mathbf{r}}}{\sqrt{V}}|u_{n\mathbf{k}}\rangle$ is the Bloch wave function of the system. If we assume a special form of the fictitious field

[TABLE]

with $\mathbf{q}=q\hat{\mathbf{e}}_{\beta}$, where $\alpha,\beta=x,y,z$ in three dimensions or $\alpha,\beta=x,y$ in two dimensions, the auxiliary quantity can be identified by picking up the appropriate Fourier component

[TABLE]

As an example, we calculate $\tilde{M}^{\theta_{y}}_{x}$ by taking $\mathbf{q}_{1}=q\hat{\mathbf{e}}_{x}$ and $h_{\alpha}(\mathbf{r})=\frac{h}{q}\sin(\mathbf{q}_{1}\cdot\mathbf{r})\delta_{\alpha,y}$. Applying perturbation theory to linear order under the Bloch representation, we find

[TABLE]

and

[TABLE]

Here it is implied that we will use the operator under Bloch representation henceforth, i.e., $\hat{H}\rightarrow H_{\mathbf{k}}=e^{-i\mathbf{k}\cdot\mathbf{r}}\hat{H}e^{i\mathbf{k}\cdot\mathbf{r}}$, $\hat{\theta}_{\alpha}\rightarrow\theta_{\alpha,\mathbf{k}}=e^{-i\mathbf{k}\cdot\mathbf{r}}\hat{\theta}_{\alpha}e^{i\mathbf{k}\cdot\mathbf{r}}$. This step is guaranteed by the requirement that the operator $\hat{\theta}_{\alpha}$ is well defined in a periodic system. By using the results above we obtain

[TABLE]

Taking the limit, we get for $m\neq n$,

[TABLE]

For $m=n$, we have

[TABLE]

Above $v_{x}=\partial_{k_{x}}H$. In total, we have

[TABLE]

The calculation of all other components is fully analogous to what we have done. The general result will be

[TABLE]

with

[TABLE]

Note the Berry curvature defined in Eq. (45) exists in both particle and hole space. Finally, by using Eq. (34) we obtain

[TABLE]

Here we used the relation $\frac{1}{\beta}\int_{0}^{\beta}d\bar{\beta}g(\bar{\varepsilon}_{n,\mathbf{k}})\bar{\varepsilon}_{n,\mathbf{k}}=\int_{0}^{\bar{\varepsilon}_{n\mathbf{k}}}d\eta g(\eta)$ with $g(\eta)=\frac{1}{e^{\bar{\beta}\eta}-1}$ and $\frac{d}{d\bar{\beta}}[\bar{\beta}g(\bar{\varepsilon}_{n,\mathbf{k}})]=g(\bar{\varepsilon}_{n,\mathbf{k}})+g^{\prime}(\bar{\varepsilon}_{n,\mathbf{k}})\bar{\varepsilon}_{n,\mathbf{k}}$.

II.2 Kubo type response

The intrinsic part of the Kubo linear response in particle-hole space is

[TABLE]

II.3 Total response

Adding the Kubo formula and dipole moment contributions together, the total response reads

[TABLE]

where $\tilde{c}_{1}(x)=\int_{0}^{x}d\eta\eta\frac{dg(\eta)}{d\eta}$ with $g(\eta)=\frac{1}{e^{\beta\eta}-1}$. Using the relation $-g(-\eta)=1+g(\eta)$, we have $\tilde{c}(x)=\tilde{c}(-x)$. Therefore, the response function can be reduced to

[TABLE]

Here we used the properties of Berry curvature shown in Eq. (51) and (54), and the relation $-\int_{\varepsilon_{n}}^{\infty}\eta\frac{dg(\eta)}{d\eta}=\frac{1}{\beta}c_{1}[g(\varepsilon_{n})]$.

II.4 Property of Berry curvature

Here we provide two useful properties of Berry curvature defined in (45).

(1) Summation rule:

[TABLE]

In the middle step, we utilized the property that the band indices $m,n$ can be interchanged.

(2) Mapping between particle and hole space.

We note that the velocity operator $v_{\mathbf{k}}$ satisfies

[TABLE]

At the same time, the particle-hole symmetry of the Hamiltonian enforces the relation

[TABLE]

which is clearly satisfied when we consider the current and source term response for a given operator $\hat{O}$. Using the particle-hole symmetry property of the eigenstates and eigenvalues, we are able to show

[TABLE]

II.5 Calculation of Berry curvature

In Fig. 4, we calculate the Berry curvature for kagome antiferromagnet KFe3(OH)6(SO4)2 using

[TABLE]

III Symmetry analysis on the spin Nernst tensor in kagome antiferromagnet

We perform a detailed analysis on the effect of the (magnetic) point group on the Nernst response tensor. Suppose the Hamiltonian respects symmetry $g$ with matrix representation $U(g)$ for unitary operation and $U(g)\mathcal{K}$ for anti-unitary operation (containing time-reversal) with $\mathcal{K}$ being the complex conjugate operator. Here $U(g)$ corresponds to the point group operation on spin mode orbitals, which is a unitary matrix that satisfies $U(g)^{\dagger}=U(g)^{T}$. On the other hand, the point group symmetries don’t mix particle and hole symmetry, such that $[\sigma_{3},U(g)]=0$. For the unitary case, we assume

[TABLE]

where $M(g)$ is the matrix acting on momentum variables. We can deduce that

[TABLE]

As a consequence, by inserting the symmetry operation in the matrix elements of an observable, we find

[TABLE]

Similarly, for the anti-unitary case,

[TABLE]

such that

[TABLE]

These relations will lead to

[TABLE]

If the operator $\hat{A}$ satisfies

[TABLE]

and we combine this with the element’s symmetry relation, we can obtain a transformation relation for the spin Nernst response coefficient

[TABLE]

where the plus and minus sign correspond to unitary and anti-unitary symmetry, and $R^{s/v}(g)$ stands for the transformation matrix for the spin and velocity operator, respectively. Moreover, suppose the involved non-magnetic point group symmetry $U(g)$ corresponds to a spatial operation with matrix form $R$ in Cartesian coordinates. If $\hat{g}$ is a unitary operation,

[TABLE]

For anti-unitary operation,

[TABLE]

Plugging Eq. (63), (64) into Eq. (62), we find

[TABLE]

In the kagome AF, we focus on two symmetries of the system: the mirror reflection with respect to the $y-z$ plane plus time-reversal $\hat{g}_{1}=\mathcal{M}_{yz}\mathcal{T}$, and the threefold rotation about the $z$ axis $\hat{g}_{2}=\mathcal{C}_{3z}$. It is straightforward to obtain the matrix representation in Cartesian coordinates of these two symmetry operations

[TABLE]

By applying these symmetries to Eq. (65), the spin Nernst response tensor (only considering in-plane driven response) can be fixed to

[TABLE]

where $\alpha^{y}_{yx}$ and $\alpha^{z}_{yx}$ correspond $\alpha_{1}$, $\alpha_{2}$ in the main text individually. Here we comment that the $g_{1}$ symmetry can also be replaced by $\mathcal{C}_{2x}\mathcal{T}$, the twofold rotation about the $x$ axis plus time reversal.