Nonlinear magnon spin Nernst effect in antiferromagnets and strain-tunable pure spin current

In this Letter, we study the spin Nernst effect (SNE) of magnons in the nonlinear response regime. We derive the formula for the nonlinear magnon spin Nernst current by solving the Boltzmann equation and find out that it is described by an extended Berry curvature dipole of magnons. The nonlinear magnon SNE is expected to occur in various N\'eel antiferromagnets without Dzyaloshinskii-Moriya interaction. In particular, the nonlinear spin Nernst current in the honeycomb and diamond lattice antiferromagnets can be controlled by strain/pressure.


I. INTRODUCTION
The Berry phase and curvature play an essential role in modern condensed matter physics; e.g., they are responsible for polarization, orbital magnetism, and various types of Hall effects [1]. In particular, it is well-known that in the linear response regime, the Hall effect [2,3] and the spin Hall effect [4][5][6][7] are described by the integral of Berry curvature (BC). Recently, transport phenomena have been further explored by taking into account the nonlinear response contributions. A remarkable study is the nonlinear Hall effect [8]. Even in materials with timereversal symmetry, transverse electric current can emerge as the second-order response to an electric field [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. It originates from a dipole moment of the BC in the crystal momentum space, named the Berry curvature dipole (BCD) [8], which appears in the systems breaking inversion and rotational symmetries.
Despite a number of studies on magnon systems, representative transport phenomena, e.g., thermal Hall effect and SNE of magnons, have been considered mostly in the magnets with the Dzyaloshinskii-Moriya interaction (DMI) or noncollinear spin configurations, which give rise to complex hopping matrix elements in the magnon * kondo-hiroki290@g.ecc.u-tokyo.ac.jp Hamiltonian. This is because their manifestation as a linear response requires the integration of the magnon BC over the whole Brillouin zone to be nonzero. In order to relax the restrictions on such magnon transport phenomena, we should explore the region beyond the linear response regime as one option. As in the case of electronic systems, by investigating the magnon nonlinear response, nontrivial transport phenomena are expected to be discovered in magnets even without the DMI or noncollinear spin configuration. So far, the nonlinear Hall effect [63], spin Seebeck effect [64], and optical response [65,66] of magnons have been proposed. However, research on the nonlinear response of magnons is just beginning to emerge.
In this paper, we study the magnon second-order response to the temperature gradient by solving the Boltzmann equation and show that the nonlinear magnon spin Nernst current can be described by an extended BCD. As shown later, the extended BCD is easily found in collinear Néel AFMs without the DMI, whereas the manifestation of BC, which is responsible for the linear magnon SNE and thermal Hall effect, requires the DMI or noncollinear spin configuration. Therefore, it is expected that various Néel AFMs even without the DMI can exhibit the nonlinear SNE of magnons, while the magnon thermal Hall effect, as well as the linear magnon SNE, is absent. As a demonstration, we apply the obtained formula to the strained honeycomb lattice AFM without the DMI. A remarkable result is that the direction of the spin current can be controlled by tuning the strain. We also investigate several other models for Néel AFMs and find the presence of the extended BCD, which results in the nonlinear magnon SNE.

II. EXPRESSION OF NONLINEAR SPIN NERNST CURRENT
First, we derive the formula for the magnon spin Nernst current up to the second-order response to the temperature gradient. We begin with the expression of the transverse magnon current in Ref. [26]: where E n (k), Ω n (k), ρ(E, T (x)) are the energy eigenvalue, BC of the nth band with the wave vector k, and a distribution function for magnons with the energy E under the temperature T (x), respectively. Here, we assume a steady state where both ends of the system are in contact with heat baths at different temperatures. In such a case, it is known that the distribution of temperature can be written as a linear function of a position [67].
Bearing this in mind, we henceforth assume that the temperature gradient is applied in the x-direction as T (x) = T 0 − x∇T [68]. The constants T 0 and ∇T denote the average temperature and the temperature gradient, respectively.
Next, we consider the Boltzmann equation to derive the formula for the nonlinear response from Eq. (1). The Boltzmann equation in the relaxation time approximation [69][70][71] is written as follows: where τ and ρ 0 (E, T (x)) are the relaxation time of magnons and the equilibrium distribution function defined as ρ 0 (E, T (x)) = [e E/T (x) −1] −1 , respectively. Here, x andk are the time-derivatives of the position and the wave vector, respectively. On the left-hand side of Eq. (2), the first and third terms drop out because we assume the steady state and the system without external field hereafter. Thus, Eq. (2) can be deformed as Here, the velocityẋ is written as (1/ )∂ kx E n (k). Replacing ∂/∂x with −∇T ∂/∂T 0 due to T (x) = T 0 − x∇T , we obtain the x-derivative of Eq. (3) up to the second order in ∇T as follows: The second term on the right-hand side of Eq. (4) vanishes in the whole space since it is an odd function of x.
Hereafter, we take the Boltzmann constant as k B = 1.
Substituting Eq. (4) to Eq. (1), we obtain the expression of the transverse magnon current as up to a second-order response to the temperature gradient ∇T : Here, by using the function c 1 (ρ 0 ) := (1 + ρ 0 )ln(1 + ρ 0 ) − ρ 0 lnρ 0 , we can write down the following equation: By using Eq. (6) and replacing the sum over k with the integral over the first Brillouin zone, the second term in the right-hand side of Eq. (5) can be rewritten as follows: In the second equality of Eq. (7), we replaced the T 0derivative acting on the function By substituting Eqs. (6) and (7) to Eq. (5), we obtain the transverse magnon current up to the second-order of the temperature gradient ∇T as follows: We note that the first term on the right-hand side has been obtained in Ref. [26]. The second term describes the nonlinear response in magnon systems.
In the following, we focus on the magnons in Néel AFMs, assuming that the spins are aligned in the zdirection. In terms of magnons, the systems have PTsymmetry due to the perfect staggered magnetization; i.e., the magnon Hamiltonian H(k) satisfies the following equation where P and T are parity and time-reversal operators, respectively. The matrix Σ z is defined as Σ z = σ z ⊗ 1 N , where σ i (i = x, y, z), N , and 1 N are the i-component of the Pauli matrix, the number of the sublattices in the unit cell, and N -dimensional identity matrix, respectively. We note that in bosonic Bogoliubovde Gennes (BdG) systems which contain pairing terms, we have to diagonalize not Hermitian matrix H(k) but non-Hermitian matrix Σ z H(k) appearing in Eq. (9) to leave the commutation relations of boson operators unchanged [72]. As long as the systems conserve the total spin in the z-direction S z , we always have two degenerate magnon states related by PT -operator.
We here write the eigenvectors and BC of magnons with the up (down) spin dipole moment as ψ n,↑(↓) (k) and Ω n,↑(↓) (k), respectively. Since each magnon excitation carries the spin angular momentum ± in systems with conservation of the total spin along the z-direction, the magnon spin Nernst current J S y = (J y,↑ − J y,↓ ) can be written as follows: where the magnon BC is defined as Ω n, , slightly different from BC in fermionic systems due to the non-Hermiticity of the matrix Σ z H(k) in the magnon BdG systems [72]. We remark that the spin Nernst current as a second-order response is described in terms of not BCD D , which we call extended BCD. One can see the structure similar to this term in the formula for the nonlinear Nernst effect of electrons [19]. We also note that the magnon thermal Hall current defined in Ref. [26] is absent because Ω n,↑ (k) = −Ω n,↓ (k) holds due to PT -symmetry.

III. MODEL
As the first model exhibiting the nonlinear magnon SNE, we consider a honeycomb lattice AFM with strain, which is illustrated in Fig. 1. The Hamiltonian of the AFM is written as follows: where S i = (S x i , S y i , S z i ) denotes the spin at site i. Here, ij 2 and ij 1 are the nearest-neighbor vertical bonds and the other ones shown in Fig. 1, respectively. The third term is an easy-axis anisotropy in the z-direction. By applying Holstein-Primakoff and Fourier transformations, we can obtain the magnon Hamiltonian as follows: Here, d and γ(k) are defined as d = 2J 1 S + J 2 S + 2κS and γ(k) = 2J 1 Se iky/2 √ 3 cos(k x /2) + J 2 Se −iky/ √ 3 , respectively, where S is the spin magnitude. The operator b ↑(↓) (k) is the annihilation operator of magnons with the spin dipole moments upward (downward), i.e., magnons from spins pointing downward (upward). The parity and time-reversal operators of this model are defined as P = 1 2 ⊗ σ x and T = K, respectively. Here, K is the complex conjugation operator. One can easily confirm that the Hamiltonian H(k) satisfies the PT -symmetry in Eq. (9).
Owing to the simple and typical model (11), we can find the candidate materials of the honeycomb AFMs. For example, the honeycomb AFMs 2-Cl-3,6-F 2 -V [73] and Mn[C 10 H 6 (OH)(COO)] 2 ×2H 2 O [74] would be candidates modeled by Eq. (11) with 0.7 < J 2 /J 1 < 1.0 and J 2 = 2J 1 , respectively. We note that the materials possess bond-dependences inherently. Thus, the nonlinear magnon SNE is expected to exhibit even without the strain. Figure 2 shows the band structure, BC, and the extended BCD of magnons with up spin dipole moment in the strained honeycomb AFM. Those of magnons with down spin dipole moment are determined by which appears due to breaking inversion and rotational symmetries. Figure 3 shows the coefficient of the nonlinear magnon SNE given by Eq. (10) in the honeycomb AFM as a function of the coupling constant J 1 . As expected, the coefficient becomes zero for J 1 = 0 and J 1 = J 2 (corresponding to the system with the three-fold rotational symmetry restored). It is noteworthy that spin current flows in the +y-direction and the −y-direction in the cases of J 1 < J 2 and J 1 > J 2 , respectively. It implies that the direction of the spin current can be controlled by tuning the strain. Here, in which directions the transverse current is driven can be understood intuitively in terms of the balance of the coupling constants of the nearest-neighbor bonds; i.e., the transverse magnon current tends to flow in the direction of the stronger nearest-neighbor bonds corresponding eventually to the +y/−y-direction in total. We also note that magnon SNE in the linear response regime is absent due to the system without the DMI [75]. In addition, the system does not exhibit the magnon thermal Hall effect due to Ω ↓ (k) = −Ω ↑ (k) and E ↓ (k) = E ↑ (k).  Fig. 1(a). In (b), (d), and (f), we take the parameters to be J 1 S = 1.5, J 2 S = 1.0, and κS = 0.01, corresponding to Fig. 1(b). The coupling constant J 2 , the easy-axis anisotropy, and the average temperature are taken to be J 2 S = 1.0, and κS = 0.01, and T 0 = 0.1, respectively. Here, we take the factor τ /( V T 0 ) to be unity.

V. ORDER ESTIMATION OF NONLINEAR SPIN NERNST CURRENT
Let us discuss the order of the nonlinear magnon SNE, by comparing it to the linear one observed in the honeycomb antiferromagnet MnPS 3 [60] with the DMI. We consider the following Hamiltonian for the antiferromagnet MnPS 3 : where the third term is the DMI between the second nearest neighbor spins. The sign convention of the DMI ξ ij is shown in Fig. 4. Other terms are the same as those in model (11). By applying Holstein-Primakoff and Fourier transformations, we obtain the magnon Hamiltonian for model (13), which is written as in the same form as in Eq. (12) with the additional term ∆(k) derived from the DMI, i.e., FIG. 4: Honeycomb lattice corresponding to the antiferromagnet MnPS 3 expressed by Eq. (13). The sign convention ξ ij = +1 (= −ξ ji ) for i → j is indicated by the orange arrows. The spin configuration is the same as in Fig. 1.
The experiment [60] shows a good agreement with the theoretical study [59], in which the parameters of MnPS 3 are considered as J 1 = J 2 = 1.54 meV, D = 0.36 meV, κS = 0.0086 meV, and S = 5/2 [76]. To estimate the order of linear spin Nernst current, we ignore the second and the third nearest-neighbor Heisenberg interactions taken into account in Ref. [59], whose coupling constants are less than one-fourth of the nearest-neighbor ones. Since the antiferromagnet MnPS 3 possesses the nonnegligible DMI, the spin current observed in MnPS 3 [60] is mainly attributed to the linear SNE, which is described by the first term in Eq. (10). Here, we write the part of the integral in this term as follows: Figure 5(a) shows the numerical result of I 1 as a function of temperature T 0 in model (13). From the figure, the order of the spin current at T 0 = 20 K can be written as follows: Next, we estimate the order of nonlinear SNE in model (13) where the DMI is set to zero, which is equivalent to model (11). We here assume that the coupling constant J 1 in Eq. (11) is changed as J 1 = 1.54 → 2.0 meV. In such a case, i.e., for D = 0 and J 1 = J 2 , the system exhibits not the linear magnon SNE but the nonlinear one, which is described by the second term in Eq. (10). Here, we define the part of the integral of this term with the prefactor 1/T 0 as I 2 , i.e., The numerical result of I 2 as a function of T 0 is shown in Fig. 5(b). From this figure, we can evaluate the nonlinear spin current at T 0 = 20 K as Since the linear spin Nernst current in Eq. (17) was observed with the electric voltage V L ∼ 1 µV through the inverse spin Hall effect at T 0 = 20 K (see Fig. 3(c) in Ref. [60]), we can estimate the voltage V N L by the nonlinear spin Nernst current in Eq. (19) in the following by taking their ratio: Here, we assume that the lifetime of magnons and the applied temperature gradient are τ ∼ 10 n s and ∇T ∼ 10 −6 K · nm −1 which is estimated by the experiment in Ref. [60], respectively. In experiments, the minimum voltage detectable through the inverse Hall effect is roughly 10 −3 µV [77], and thus we can detect the nonlinear SNE if the magnon lifetime is τ 1 ns. In model (11), the velocity of magnons at the Γ point is estimated as v = (∂/ ∂k i )E ↑ (k) ∼ 10 12 nm · s −1 . Then, the corresponding mean free path for V N L ∼ 10 −3 µV (τ ∼ 1 ns) is l ∼ 1 µm, which is achievable in magnets. , the other parameters J 2 , κ, and S are taken as J 2 = 1.54 meV, κS = 0.0086 meV, and S = 5/2, respectively. We note that I 1 is dimensionless.

VI. NONLINEAR MAGNON SPIN NERNST EFFECT IN VARIOUS ANTIFEROMAGNETS
In this part, we showcase the magnon extended BCD and nonlinear spin Nernst current in several AFMs; square lattice AFMs with bond dependences and diamond lattice AFM under pressure [see Fig. 6]. The forms of the Hamiltonians in these three cases are the same as Eq. (11). Figure 6 shows which nearest-neighbor bonds correspond to ij 1 and ij 2 in the Hamiltonians. By applying Holstein-Primakoff and Fourier transformations to Eq. (11), we obtain magnon Hamiltonians for these AFMs which are the same form as Eq. (12) with different d and γ(k). In the case of the square lattice AFMs with staggered-bond dependence [see Fig. 6(a)], d and γ(k) in Eq. (12) is defined as follows: Those of the square lattice AFMs with zigzag-bond dependence [see Fig. 6(b)] are given by In the case of diamond lattice AFM [see Fig. 6(c)], d and γ(k) are written as follows: where a 0 = 0, 1 , and a 3 = 0, 0, 3 4 2 3 . One can easily make sure that these AFMs also have PTsymmetry thanks to the same form of the magnon Hamiltonian as that expressed in Eq. (12).
in the k x = 0 plane, respectively. They are relevant for nonlinear magnon SNE in three dimensions, as seen later [Eq. (27)]. As in the case of honeycomb lattice AFM, the energy eigenvalue, BC, and the extended BCD of magnons with down spin dipole moment in the square (diamond) lattice AFM are deter- , respectively. As shown in Fig. 7, the extended BCD mainly appears around the Γ point, which contributes to the nonlinear SNE. Figure 8 shows the coefficients of nonlinear magnon SNE in these systems as a function of the coupling constant J 2 . For the case of the diamond lattice AFM, we generalize the formula Eq. (10) to that in three dimensions; i.e., As shown in Fig. 8, the nonlinear SNE of magnons occurs in the square and diamond lattice AFMs while the coefficients become zero for the symmetric case, i.e., J 1 = J 2 . We here emphasize that Fig. 8(c) implies the pressuretunable spin current can generate in the diamond lattice AFM as well as the honeycomb lattice AFM. The quasi-two-dimensional antiferromagnet Cu(en)(H 2 O) 2 SO 4 is a candidate material of the square lattice antiferromagnet with zigzag bond dependence, in which the coupling constants are estimated as J 1 = 10J 2 ∼ 0.30 meV [78]. We can also find candidate materials of the diamond lattice antiferromagnets exhibiting Néel order, such as CoRh 2 O 4 [79], MnAl 2 O 4 [80], and (ET)Ag 4 (CN) 5 [81]. In particular, (ET)Ag 4 (CN) 5 is a molecular compound, and can be distorted by applying pressure without difficulty. Thus, the pressure-tunable spin current is expected in the AFM.

VII. SUMMARY
In this paper, we have derived the formula for the magnon spin Nernst current as a second-order response and found that it is characterized by the extended BCD. We then have applied the obtained formula to the strained honeycomb lattice AFMs and found out that the direction of spin current can be controlled by tuning the strain. We have also calculated the extended BCD of magnons and confirmed that the nonlinear magnon SNE appears in the square lattice AFMs with bond dependences and the diamond lattice AFM under pressure. Even without the DMI, the nonlinear magnon SNE is expected to be brought about in various Néel AFMs when the inversion and rotational symmetries are broken by such as strain/pressure. Here, we note that distortion can induce the DMI, and thus the systems have the potential for exhibiting the linear SNE. However, if materials consist of light elements, they can show mainly not the linear but nonlinear SNE because the strain-induced DMI would be negligible.
Our study reveals that the pure spin current can be generated in various Néel AFMs. Owing to the simple setup, we can find a number of candidate antiferromagnetic materials exhibiting nonlinear magnon SNE; e.g., AFMs on a honeycomb lattice [73,74,[82][83][84][85][86], square lattice with zigzag bond dependence [78], and diamond lattice [79][80][81]. We here emphasize that our proposal for the nonlinear magnon SNE provides one of a few possible ways to generate the spin current in materials composed of light elements such as organic materials [87,88] where the DMI is negligible. Since such organic materials [73,81] are easily deformed mechanically, we expect to have a controllable pure spin current by strain or pressure, which expands the possibilities of applications in spintronics.