Intrinsic spin Nernst effect of magnons in a noncollinear antiferromagnet

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.

Many spintronics concepts also apply to antiferromagnets [31]. In particular, collinear antiferromagnets can exhibit the spin Seebeck effect [32], spin pumping [33], spin-orbit torque [34], the spin Nernst effect [25][26][27][28]30], etc. Noncollinear antiferromagnets (NAFMs) have attracted considerable attention in recent years, as such systems support nontrivial band structure topology. The anomalous Hall effect [35] and spin Hall effect [36,37] have been realized in Mn 3 X (X = Ge, Sn, Ga, Ir, Rh, or Pt) systems, where electrons act as charge or spin carriers. Furthermore, the thermal Hall effect, mediated by magnons, is also identified in NAFM insulators [30,[38][39][40][41]. Nevertheless, the magnon-mediated spin transport in NAFMs [42][43][44] has not yet been well explored, especially in the context of the topology of magnon bands. 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 better control of the spin current. NAFMs typically possess different ground states [45][46][47] 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 [48]. Studies of spin currents in such systems can provide a new venue for probing these materials [49]. Motivated by these interesting possibilities, we initiate a discussion on the magnon-mediated spin current physics in noncollinear antiferromagnets and hope to stimulate subsequent research on, e.g., spin transport in topological magnon insulators [50], optical generation of magnon-mediated spin currents [51,52], 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 [53]. Magnons driven by temperature gradients require accounting for the effects associated with the orbital magnetization [25,54]. We start by discussing the definition of spin current in particlehole space by following Refs. [25,53], where spin nonconservation 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, where J γ λ is the spin current with polarization γ , ( j S n,k ) γ λ β is the spin Berry curvature of magnons, and c 1 ( is an auxiliary function stemming from the Bose-Einstein statistics of magnons. We apply our theory to the kagome antiferromagnet KFe 3 (OH) 6 (SO 4 ) 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.

II. 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 where 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.
Because of Bose-Einstein statistics, the eigenvalue problem has to be solved for the matrix σ 3 H k [55], where here and in what follows we use Pauli matrices σ i acting in the particle-hole space. Here H k is the Hamiltonian matrix in momentum space, which can be diagonalized by a paraunitary matrix T k , i.e., is the matrix describing eigenvalues and T k satisfies T † k σ 3 T k = σ 3 . We now write a general theory applicable to bosonic systems where the Bloch wave function corresponding to the band dispersion ε n,k is given by |ψ n,k = e ik·r |u n,k . We can then introduce a notation [56] σ 3 H k u R n,k =ε n,k u R n,k , u L n,k σ 3 H k =ε n,k u L n,k , where in terms of the magnonic Hamiltonian |u R n,k = T n,k and u L n,k | = T † n,k σ 3 are the right and left eigenstates of the pseudo-Hermitian Hamiltonian, andε n,k = (σ 3 E k ) nn . Hereafter, we will only refer to the right eigenstates |u R n,k = |u n,k . The normalization relation reads u n,k |σ 3 |u m,k = (σ 3 ) nm . Moreover, the Hamiltonian (2) possesses particle-hole symmetry (PHS) so that the Hamiltonian obeys σ 1 H k σ 1 = H * −k , which leads to relationsε n+N,k = −ε n,−k and |u n,k = e iφ n σ 1 |u n+N,−k * , where φ n is a redundant phase factor.
Because the temperature gradient is a statistical force and does not directly enter the Hamiltonian, we introduce a perturbation corresponding to a pseudogravitational potential, χ (r), to account for the temperature gradient [25,54,57], With the perturbation, the total Hamiltonian is amended to H= 1 2 dr˜ † (r)Ĥ˜ (r), where˜ (r)=(1 + r · ∇χ/2) (r). To linear order, the system will respond to a temperature gradient in the same way as to a perturbation with χ (r) = −T (r)/T . We now introduce an arbitrary matrixÔ and a local observable O(r) = 1 2 † (r)Ô (r). In what follows, we will mostly considerÔ =ˆ α , which corresponds to the magnon spin density operator given byˆ α = −σ 0 ⊗ Diag(n α 1 , . . . , n α N ), where α = x, y, z, σ 0 describes the particlehole space, and n α i (i = 1, . . . , N ) is the projection ofn i along α-axis with unit vectorn i being the ground-state direction of spin at position i in each unit cell. The time evolution of this operator can be obtained from the Heisenberg equation applied to the total Hamiltonian (details in Appendix A) [25] Here j O =˜ † (r)ĵ O˜ (r) and S O =˜ † (r)Ŝ O˜ (r) correspond to the local current and source densities, respectively, To linear order in the temperature gradient, the above densities are explicitly decomposed as ρ θ = where for θ one needs to substitute either j O or S O . We will use a four-vector convention in which θ 0 = S O and θ = j O . The nonvanishing source term indicates the nonconservation of the observable; for instance, when O(r) corresponds to spin density, the source term represents torque density. We note in passing that the source term dipole P O can be defined as Thus, a conserved current can be defined as J O = j O + P O to restore the continuity equation [53]. The current term j O coincides with the conventional definition in the literature of the spin Hall effect [24]. In general, based on Eq. (5) we can interpret 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 α The thermal response to linear order in the temperature gradient reads 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 α = (S θ α β + M θ α β )∇ β χ , where S θ α β and M θ α β correspond to the first and second terms in Eq. (6).
In the spirit of the Kubo response calculation [4,25], the nonequilibrium part can be described by Here where ω m is the bosonic Matsubara frequency. J q is the averaged heat current operator defined as . This heat current expression can be inferred from the energy conservation equationρ E + ∇ · j q = 0, where ρ E is the energy density of the system. After performing the linear response calculation, the intrinsic nonequilibrium coefficient reads Here where (. . . ) nm = u n,k | . . . |u m,k and g(ε nk ) is the Bose-Einstein distribution. Here ( θ n,k ) α β is the generalized Berry curvature calculated for operatorθ α . This Berry curvature respects the sum rule 2N n=1 ( θ n,k ) α β = 0, and PHS results in The contribution corresponding to ρ [1] θ α is expressed as To calculate this term, we can identify a thermodynamic expression for M θ α β by following Refs. [58][59][60][61][62]. We introduce a perturbation coupled with a four-component fictitious field If the field varies very slowly on the scale of the lattice constant, we can identify a thermodynamic expression where 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 being an auxiliary quantity and K = + T S. We assume that the fictitious field takes the form h α (r) = (h 0 α /q) sin(q · r), with q = qê β (β = x, y, z in three dimensions and β = x, y in two dimensions). The auxiliary quantity is calculated by picking up the appropriate Fourier component where δK (r) is the variation due to the fictitious field, which can be obtained from perturbation theory [59]. Combining Eqs. (13) and (12), we obtain (see details in Appendix B) By combining the nonequilibrium part in Eq. (8) with Eq. (14) and canceling the orbital part (corresponding to a bound current), we obtain the thermal response formula which constitutes the main result of this paper: Note that we express our result using particle bands (n 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 k = H −k , the Berry curvature of the torque term satisfies Together with the relation ε n,k = ε n,−k , this results in the vanishing of all torque components in Eq. (15).

III. 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ĵ γ ,λ = 1 4 (v λ σ 3ˆ γ +ˆ γ σ 3vλ ) 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 where the first and third terms represent nearest and secondnearest 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 D i j = D pni j + D zẑ , where D p and D z correspond to the in-plane and out-of-plane DMI strength, andn i j is 013079-3 an in-plane unit vector corresponding to the direction of the in-plane DMI. The in-plane DMI can only arise when M z symmetry is broken [35]; 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 outof-plane canting angle η to spin order with magnitude η = ) [30]. Here we consider the q = 0 phase with spin order as shown in Fig. 1. The magnetic moments orient according ton i = (cos η cos φ i , cos η sin φ i , sin η), where φ i is the angle formed by the in-plane projection of moment with the x axis. Specifically, φ A = π/2, φ B = 7π/6, and φ C = −π/6. For the spin Nernst response, we iden-tifyÔ discussed above as the spin operator in the magnon basis (r) . The spin conductivity tensor of a spin-polarized current in a noncollinear antiferromagnet [37,63] is restricted to a certain form by the magnetic space group of the system. Given that the intrinsic spin Nernst tensor in relation J γ λ = α γ λβ ∇ β T is even under the time reversal transformation, the symmetry constraints become where the matrix R represents a symmetry element R (in Cartesian coordinates) entering the antiunitary symmetry RT or unitary symmetry R of the system (see Appendix C). As an example, we focus on a system with two symmetries: mirror reflection with respect to the y-z plane combined with time reversal, M x T , and threefold rotation about the z axis, C 3z . The shape of the spin Nernst tensor corresponding to the constraints in Eq. (17) becomes Here, the M x T symmetry can be replaced by C 2x 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 Mn 3 X (X = Rh, Ir, or Pt) [36].
We apply our theory to a single layer of potassium iron jarosite, KFe 3 (OH) 6 (SO 4 ) 2 , for which the material parameters are J 1 = 3.18 meV, J 2 = 0.11 meV, |D p |/J 1 = 0.062, and D z /J 1 = −0.062 [30,64]. 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 [65]. The numerically obtained form of the spin Nernst conductivities agrees with Eq. (18). 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, η = 1.9 • [30]. By applying a 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 threedimensional structures as a temperature gradient of 20 K/mm should result in a spin current of the order of 10 −11 J/m 2 according to Fig. 3, where α 3D = α/c, with c being the interlayer distance. Finally, we note that the spin Nernst effect reported in Ref. [42] differs from the intrinsic effect reported here as the former has the symmetry of the extrinsic effect.

ACKNOWLEDGMENTS
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.

APPENDIX A: TIME EVOLUTION OF A LOCAL OBSERVABLE
To derive the time evolution equation for a local observable O(r) = 1 2 † (r)Ô (r), we first prepare a basic knowledge of the Hamiltonian operator and commutators in particle-hole space by following Refs. [25,54]. The total Hamiltonian can be generally expressed as H = 1 2 dr˜ † (r)Ĥ˜ (r) withĤ = δ H δ e ip·δ , in which e ip·δ is the translation operator that satisfies e ip·δ f (r) = f (r + δ). Here δ is the vector shift between unit cells,˜ (r) = (1 + r · ∇χ/2) (r). Based on the basic commutators between bosons [a i (r), a † j (r )] = δ i j δ r,r , [a i (r), a j (r )] = 0, we can construct commutators in the particle-hole basis where σ 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 follows, Here we used the simplified notation ξ (r) = 1 + r · ∇χ/2. We also took advantage of particle-hole symmetry, i.e., n (r) = (σ 1 ) nl † l (r) and σ 1Ô σ 1 =Ô, 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: Here we used the notationv δ = iδH δ e ip·δ andv = i δ δH δ e ip·δ = i[Ĥ, r]. In the last line, we take the limit δ → 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.

APPENDIX B: 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 where S θ α β and M θ α β 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.

Dipole moment contribution
In the main text, we have shown the relation . The thermodynamic definition of grand potential reads = E − T S − μN, where S, μ, and N are the entropy, chemical potential, and particle number of the system and E is the energy, which reads Here we use the relation † n,k m,k = (σ 3 ) nm g(ε n,k ) with m,k = l (T k ) ml 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 from which we can identify the Maxwell relation To get rid of calculations involving entropy S, we first introduce an auxiliary quantitỹ with K = + T S = E (μ = 0). By utilizing Eq. (B4), we obtain and hence the dipole moment contribution can be calculated asM If we regard the fictitious field term as a perturbation, the variation of K to linear order reads where |ψ nk = e ik·r √ V |u nk is the Bloch wave function of the system. If we assume a special form of the fictitious field with q = qê β , where α, β = x, y, z in three dimensions or α, β = x, y in two dimensions, the auxiliary quantity can be identified by picking up the appropriate Fourier component As an example, we calculateM θ y x by taking q 1 = qê x and h α (r) = h q sin(q 1 · r)δ α,y . Applying perturbation theory to linear order under the Bloch representation, we find Here it is implied that we will use the operator under Bloch representation henceforth, i.e.,Ĥ → H k = e −ik·rĤ e ik·r , θ α → θ α,k = e −ik·rθ α e ik·r . This step is guaranteed by the requirement that the operatorθ α is well defined in a periodic system. By using the results above, we obtaiñ Taking the limit, we get for m = n, [g(ε mk )ε mk − g(ε n,k )ε n,k ](σ 3 ) nn (σ 3 ) mm i u n,k |σ 3 |∂ k x u m,k u m,k |θ y |u n,k ε n,k −ε m,k + c.c.
The calculation of all other components is fully analogous to what we have done. The general result will bẽ Note the Berry curvature defined in Eq. (B18) exists in both particle and hole space. Finally, by using Eq. (B7) we obtain Here we used the relation 1

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

Total response
When we add the Kubo formula and dipole moment contributions together, the total response reads Using the relation −g(−η) = 1 + g(η), we havec(x) =c(−x). Therefore, the response function can be reduced to Here we used the properties of Berry curvature shown in Eqs. (B24) and (B27) and the relation − ∞ ε n η dg(η) dη = 1 β c 1 [g(ε n )].

Property of Berry curvature
Here we provide two useful properties of Berry curvature defined in (B18).
(1) Summation rule: 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 k satisfies At the same time, the particle-hole symmetry of the Hamiltonian enforces the relation which is clearly satisfied when we consider the current and source term response for a given operatorÔ. Using the particle-hole symmetry property of the eigenstates and eigenvalues, we are able to show

Calculation of spin Berry curvature
In Fig. 4, we calculate the spin Berry curvature for kagome antiferromagnet KFe 3 (OH) 6 (SO 4 ) 2 using

APPENDIX C: 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)K for antiunitary operation (containing time reversal) with 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) † = U (g) T . On the other hand, the point group symmetries do not mix particle and hole symmetry, such that [σ 3 , U (g)] = 0. For the unitary case, we assume where M(g) is the matrix acting on momentum variables. We can deduce that |ψ n,M(g)k = U (g)|ψ n,k , ε M(g)k = ε k .
As a consequence, by inserting the symmetry operation in the matrix elements of an observable, we find If the operatorÂ satisfies and we combine this with the element's symmetry relation, we can obtain a transformation relation for the spin Nernst response coefficient where the plus and minus signs correspond to unitary and antiunitary symmetry and R s/v (g) stands for the transformation matrix for the spin and velocity operator, respectively.
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ĝ 1 = M yz T , and the threefold rotation about the z axisĝ 2 = C 3z . It is straightforward to obtain the matrix representation in Cartesian coordinates of these two symmetry operations: (C12) By applying these symmetries to Eq. (C11), the spin Nernst response tensor (only considering in-plane driven response) can be fixed to where α y yx and α z yx correspond α 1 , α 2 in the main text individually. Here we comment that the g 1 symmetry can also be replaced by C 2x T , the twofold rotation about the x axis plus time reversal.