Effects of magnetic anisotropy on spin and thermal transports in classical antiferromagnets on the square lattice

TL;DR

This study investigates how magnetic anisotropy influences spin and thermal transport in classical antiferromagnets on a square lattice, revealing distinct behaviors at phase transitions and highlighting the topological nature of the Kosterlitz-Thouless transition.

Contribution

It provides the first detailed analysis of spin and thermal transport across different anisotropy-induced phase transitions in the classical XXZ model using hybrid Monte-Carlo and spin-dynamics simulations.

Findings

Spin current conductivity diverges at the KT transition with exponential behavior.

Thermal current does not show significant anomalies at phase transitions.

Enhanced spin transport near the KT transition is linked to vortex dynamics.

Abstract

Transport properties of the classical antiferromagnetic XXZ model on the square lattice have been theoretically investigated, putting emphasis on how the occurrence of a phase transition is reflected in spin and thermal transports. As is well known, the anisotropy of the exchange interaction $\Delta\equiv J_z/J_x$ plays a role to control the universality class of the transition of the model, i.e., either a second-order transition at $T_N$ into a magnetically ordered state or the Kosterlitz-Thouless (KT) transition at $T_{KT}$, which respectively occur for the Ising-type ($\Delta >1$) and $XY$-type ($\Delta <1$) anisotropies, while for the isotropic Heisenberg case of $\Delta=1$, a phase transition does not occur at any finite temperature. It is found by means of the hybrid Monte-Carlo and spin-dynamics simulations that the spin current probes the difference in the ordering properties, while the thermal current does not. For the $XY$-type anisotropy, the longitudinal spin-current conductivity $\sigma^s_{xx}$ ($=\sigma^s_{yy}$) exhibits a divergence at $T_{KT}$ of the exponential form, $\sigma^s_{xx} \propto \exp\big[ B/\sqrt{T/T_{KT}-1 }\, \big]$ with $B={\cal O}(1)$, while for the Ising-type anisotropy, the temperature dependence of $\sigma^s_{xx}$ is almost monotonic without showing a clear anomaly at $T_{N}$ and such a monotonic behavior is also the case in the Heisenberg-type spin system. The significant enhancement of $\sigma^s_{xx}$ at $T_{KT}$ is found to be due to the exponential rapid growth of the spin-current-relaxation time toward $T_{KT}$, which can be understood as a manifestation of the topological nature of a vortex whose lifetime is expected to get longer toward $T_{KT}$. Possible experimental platforms for the spin-transport phenomena associated with the KT topological transition are discussed.

Equations

$${\cal H}=-J\sum_{\langle i,j\rangle}\Big{(}S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j}+\Delta S^{z}_{i}S^{z}_{j}\Big{)},$$

$$\frac{d{\mathbf{S}}_i}{dt}$$

$${\mathbf{H}}_i^{\mathrm{eff}}$$

$$\int d\mathbf{r}\mathbf{r}\frac{\partial }{\partial t}\mathcal{O}(\mathbf{r},t)=-\int d\mathbf{r}\mathbf{r}\nabla \cdot {\mathbf{j}}_{\mathcal{O}}(\mathbf{r},t)=\int d\mathbf{r}{\mathbf{j}}_{\mathcal{O}}(\mathbf{r},t).$$

$${\mathbf{J}}_{\mathcal{O}}(t)=\int d\mathbf{r}{\mathbf{j}}_{\mathcal{O}}(\mathbf{r},t)=\int d\mathbf{r}\mathbf{r}\frac{\partial }{\partial t}\mathcal{O}(\mathbf{r},t).$$

$${\mathbf{J}}_s^z(t)$$

$${\mathbf{J}}_{th}(t)$$

$$\left(\begin{array}[]{c}{\bf j}_{s}\\ {\bf j}_{th}\end{array}\right)=\left(\begin{array}[]{cc}L^{s,s}&L^{s,th}\\ L^{th,s}&L^{th,th}\end{array}\right)\left(\begin{array}[]{c}\nabla H\\ \nabla T/T\end{array}\right)$$

$${\sigma }^s=L^{s,s},\kappa =T^{-1}{L}^{th,th}.$$

$$L^{a,b}_{\mu\nu}(\omega)=\int_{0}^{\infty}dt\,e^{-i\omega t-\eta t}\int_{0}^{1/T}d\lambda\,\big{\langle}j_{a,\nu}(-i\hbar\lambda)j_{b,\mu}(t)\big{\rangle},$$

$$L^{a,b}_{\mu\nu}(0)=\frac{1}{T}\int_{0}^{\infty}dt\,\big{\langle}j_{a,\nu}(0)\,j_{b,\mu}(t)\big{\rangle}.$$

$${\sigma }_{\mu \nu }^s$$

$${\kappa }_{\mu \nu }$$

$$\left\{\begin{array}[]{l}S^{x}_{i}=\tilde{S}^{z}_{i}\sin(\theta_{i})+\tilde{S}^{x}_{i}\cos(\theta_{i})\\ S^{z}_{i}=\tilde{S}^{z}_{i}\cos(\theta_{i})-\tilde{S}^{x}_{i}\sin(\theta_{i})\\ S^{y}_{i}=\tilde{S}^{y}_{i}\end{array}\right.,$$

$${\cal H}=\frac{J}{2}\sum_{i}\sum_{j\in N(i)}\Big{[}\tilde{S}^{x}_{i}\tilde{S}^{x}_{j}-\tilde{S}^{y}_{i}\tilde{S}^{y}_{j}+\Delta\tilde{S}^{z}_{i}\tilde{S}^{z}_{j}\Big{]}.$$

$$\left\{\begin{array}[]{l}\tilde{S}^{z}_{i}=S-\hat{a}^{\dagger}_{i}\hat{a}_{i}\\ \tilde{S}^{x}_{i}+i\tilde{S}^{y}_{i}=\sqrt{2S}\Big{(}1-\frac{\hat{a}^{\dagger}_{i}\hat{a}_{i}}{2S}\Big{)}^{\frac{1}{2}}\hat{a}_{i}=\sqrt{2S}\,\hat{a}_{i}+{\cal O}(S^{-\frac{1}{2}})\\ \tilde{S}^{x}_{i}-i\tilde{S}^{y}_{i}=\sqrt{2S}\hat{a}^{\dagger}_{i}\Big{(}1-\frac{\hat{a}^{\dagger}_{i}\hat{a}_{i}}{2S}\Big{)}^{\frac{1}{2}}=\sqrt{2S}\,\hat{a}^{\dagger}_{i}+{\cal O}(S^{-\frac{1}{2}})\\ \end{array}\right.$$

$${\hat{a}}_i^{†}=\frac{1}{\sqrt{N}}\sum _{\mathbf{q}}{\hat{a}}_{\mathbf{q}}^{†}{e}^{-i\mathbf{q}\cdot {\mathbf{r}}_i},{\hat{a}}_i=\frac{1}{\sqrt{N}}\sum _{\mathbf{q}}{\hat{a}}_{\mathbf{q}}{e}^{i\mathbf{q}\cdot {\mathbf{r}}_i},$$

$$\mathcal{H}$$

$$\left\{\begin{array}[]{l}\hat{a}_{\bf q}=u_{\bf q}\,\hat{b}_{\bf q}+v_{\bf q}\,\hat{b}^{\dagger}_{-{\bf q}},\\ u_{\bf q}=u_{-{\bf q}}=\frac{1}{2}\Big{[}\Big{(}\frac{A_{\bf q}+B_{\bf q}}{A_{\bf q}-B_{\bf q}}\Big{)}^{1/4}+\Big{(}\frac{A_{\bf q}-B_{\bf q}}{A_{\bf q}+B_{\bf q}}\Big{)}^{1/4}\Big{]},\\ v_{\bf q}=v_{-{\bf q}}=\frac{1}{2}\Big{[}\Big{(}\frac{A_{\bf q}+B_{\bf q}}{A_{\bf q}-B_{\bf q}}\Big{)}^{1/4}-\Big{(}\frac{A_{\bf q}-B_{\bf q}}{A_{\bf q}+B_{\bf q}}\Big{)}^{1/4}\Big{]},\end{array}\right.$$

$$\mathcal{H}$$

$$A_{\mathbf{q}}$$

$$B_{\mathbf{q}}$$

$${\gamma }_{\mathbf{q}}$$

$${\bf J}_{th}=\big{(}4|J|S\big{)}^{2}\sum_{\bf q}\tilde{\varepsilon}_{\bf q}\,\tilde{{\bf v}}_{\bf q}\,\hat{b}_{\bf q}^{\dagger}\hat{b}_{\bf q}+{\cal O}\big{(}S^{1}\big{)},$$

$${\bf J}^{z}_{s}=\left\{\begin{array}[]{l}2|J|S\sum_{\bf q}\tilde{{\bf v}}_{\bf q}\Big{[}\frac{A_{\bf q}}{B_{\bf q}}(\hat{b}_{\bf q}^{\dagger}\hat{b}_{-{\bf q}+{\bf Q}}^{\dagger}+\hat{b}_{\bf q}\hat{b}_{-{\bf q}+{\bf Q}})\\ \qquad-(\hat{b}_{\bf q}^{\dagger}\hat{b}_{{\bf q}+{\bf Q}}+\hat{b}_{\bf q}\hat{b}_{{\bf q}+{\bf Q}}^{\dagger})\Big{]}+{\cal O}\big{(}S^{0}\big{)}\quad(\Delta\geq 1)\\ {\cal O}\big{(}S^{1/2}\big{)}\qquad\qquad\qquad\qquad\qquad\qquad\quad(\Delta<1),\\ \end{array}\right.$$

$${\tilde{\epsilon }}_{\mathbf{q}}$$

$${\tilde{\mathbf{v}}}_{\mathbf{q}}$$

$$\displaystyle\big{\langle}J_{th,\nu}(0)\,J_{th,\mu}(0)\big{\rangle}=\sum_{{\bf q},{\bf q}^{\prime}}\varepsilon_{\bf q}\varepsilon_{{\bf q}^{\prime}}v_{{\bf q},\mu}v_{{\bf q}^{\prime},\nu}\big{\langle}\hat{b}_{\bf q}^{\dagger}\hat{b}_{\bf q}\hat{b}_{{\bf q}^{\prime}}^{\dagger}\hat{b}_{{\bf q}^{\prime}}\big{\rangle}$$

$$\displaystyle\qquad=\delta_{\mu,\nu}\sum_{\bf q}\big{[}\varepsilon_{{\bf q}}\,v_{{\bf q},\mu}\big{]}^{2}f_{\rm B}(\varepsilon_{\bf q})\big{[}1+2f_{\rm B}(\varepsilon_{\bf q})\big{]},$$

$$\displaystyle\big{\langle}\hat{b}_{\bf q}^{\dagger}\hat{b}_{\bf q}\hat{b}_{{\bf q}^{\prime}}^{\dagger}\hat{b}_{{\bf q}^{\prime}}\big{\rangle}=\frac{T^{2}}{Z}\frac{\partial^{2}\,Z}{\partial\varepsilon_{\bf q}\partial\varepsilon_{{\bf q}^{\prime}}},$$

$$\displaystyle Z={\rm Tr}\Big{[}\exp\big{(}-\frac{1}{T}\sum_{\bf q}\varepsilon_{\bf q}\hat{b}_{\bf q}^{\dagger}\hat{b}_{\bf q}\big{)}\Big{]}=\prod_{\bf q}\big{[}-f_{\rm B}(-\varepsilon_{\bf q})\big{]}$$

$$f_{\mathrm{B}}(x)\to \frac{T}{x},$$

$$\big{\langle}J_{th,\nu}(0)\,J_{th,\mu}(0)\big{\rangle}_{\rm cl}=\delta_{\mu,\nu}2\,T^{2}\,\sum_{\bf q}\big{[}v_{{\bf q},\mu}\big{]}^{2}.$$

$$\varepsilon_{\bf q}\simeq\left\{\begin{array}[]{l}\sqrt{\Delta_{gp}^{2}+\frac{(4|J|S)^{2}}{2}|{\bf q}|^{2}}\qquad(\Delta>1)\\ 2\sqrt{2}|J|S\,|{\bf q}|\qquad\qquad\quad(\Delta=1)\\ 2|J|S\sqrt{1+\Delta}\,|{\bf q}|\qquad\quad(\Delta<1)\\ \end{array}\right..$$

Full text

I Introduction

Transport phenomena in magnetic systems reflect dynamical properties of interacting spins, such as magnetic excitations and fluctuations. Of recent particular interest is spin transport which is becoming available as a probe to study magnetic properties thanks to the development of experimental methods in the context of spintronics Spincurrent-mag_Frangou_16 ; Spincurrent-mag_Qiu_16 ; Spincurrent-mag_Wang_17 ; Spincurrent-mag_Frangou_17 ; Spincurrent-mag_Gladii_18 ; Spincurrent-mag_Ou_18 . This demands to explore the fundamental physics underlying the association between the spin transport and magnetic phase transitions. In this paper, we theoretically investigate transport properties of two-dimensional antiferromagnetic insulators, putting emphasis on the effects of magnetic anisotropy which plays a role of controlling the universality class of the system.

A minimal model of two-dimensional antiferromagnets with magnetic anisotropy would be the classical nearest-neighbor (NN) antiferromagnetic XXZ model on the square lattice. The spin Hamiltonian of the system is given by

[TABLE]

where $S^{\alpha}_{i}$ is $\alpha$-component of a classical spin ${\bf S}_{i}$ at a lattice site $i$, $\langle i,j\rangle$ denotes the summation over all the NN pairs, $J<0$ is an antiferromagnetic exchange interaction, and $\Delta>0$ is a dimensionless parameter characterizing the magnetic anisotropy. The ground state of this system is the conventional two-sublattice antiferromagnetic order, whereas the finite-temperature properties depend on the magnetic anisotropy $\Delta$. In the isotropic case of $\Delta=1$, Eq. (1) is nothing but the isotropic Heisenberg model, so that a phase transition does not occur at any finite temperature. In the anisotropic case of $\Delta>1$ ($\Delta<1$), the system belongs to the Ising ($XY$) universality class and exhibits a magnetic (Kosterlitz-Thouless topological KT_KT_73 ) phase transition at a finite temperature $T_{N}$ ($T_{KT}$). The purpose of this work is to clarify how the difference in the ordering properties among the three cases, $\Delta>1$, $\Delta<1$, and $\Delta=1$, is reflected in the transport properties. Our main focus is on whether a signature of a phase transition shows up in the spin and thermal transports or not.

In the ordered phase at lower temperatures, spin and thermal currents should be carried by spin waves or magnons. With increasing temperature, thermally-activated nontrivial excitations and fluctuations would come into play. In particular, in the case of the $XY$-type anisotropy ($\Delta<1$), free vortices dissociated at higher temperature above $T_{KT}$ may strongly affect the current relaxation, because the topological object of the vortex is generally robust against weak perturbations, resulting in a relatively long lifetime compared with the damping of the spin-wave mode SpnDyn-XY_Huber_82 ; SpnDyn-XY_Mertens_87 ; SpnDyn-XY_Gouvea_89 ; SpnDyn-XY_Evertz_96 . As we will demonstrate below, this is actually the case for the spin-current relaxation. In this paper, we will investigate temperature dependences of the conductivities of the spin and thermal currents in the Ising-type ($\Delta>1$), $XY$-type ($\Delta<1$), and Heisenberg-type ($\Delta=1$) spin systems by means of the hybrid Monte-Carlo (MC) and spin-dynamics simulations.

Our result is summarized in Fig. 1. The longitudinal thermal conductivity $\kappa_{xx}$, which is the response to the temperature gradient $\nabla T$, is insensitive to the difference in the ordering properties. $\kappa_{xx}$ increases toward $T=0$ as a power function of temperature $T$ in all the three cases of $\Delta>1$, $\Delta<1$, and $\Delta=1$, without showing a clear anomaly at $T_{N}$ and $T_{KT}$. In contrast, the longitudinal spin-current conductivity $\sigma^{s}_{xx}$, which is the response to the magnetic-field gradient $\nabla H$, exhibits temperature dependences characteristic of the three different universality classes. For the $XY$-type anisotropy, $\sigma^{s}_{xx}$ exhibits a divergent sharp peak at $T_{KT}$, while for the Ising-type anisotropy, the temperature dependence of $\sigma^{s}_{xx}$ is monotonic without showing a clear anomaly at $T_{N}$. In the Heisenberg-type isotropic case, $\sigma^{s}_{xx}$ shows an exponential increase toward $T=0$. The significant enhancement of $\sigma^{s}_{xx}$ at $T_{KT}$ is due to the spin-current relaxation getting slower toward $T_{KT}$, which can be understood as a manifestation of the topological nature of the vortex whose lifetime is expected to get longer toward $T_{KT}$.

This paper is organized as follows: In Sec. II, the theoretical framework for transport phenomena in magnetic insulators will be given. We derive the expressions for the conductivities of the spin and thermal currents, and explain the details of our simulations. In Sec. III, low-temperature transport properties will be discussed based on analytical calculations within the linear spin-wave theory. Numerical results on the thermal and spin transports will be shown in Secs. IV and V, respectively. We end this paper with summary and discussions in Sec. VI.

II Theoretical framework for transport phenomena in magnets

In this section, starting from the introduction to the equation of motion for the spin dynamics, we first derive the spin and thermal currents by using this spin-dynamics equation, and then, formulate the spin-current conductivity $\sigma^{s}_{\mu\nu}$ and the thermal conductivity $\kappa_{\mu\nu}$ within the linear response theory. Subsequently, we will explain numerical methods to integrate the equation of motion, taking account of temperature effects.

II.1 Spin dynamics

The spin dynamics, i.e., the time evolution of the spins for the Hamiltonian (1), is determined by the following semiclassical equation of motion:

[TABLE]

where $N(i)$ denotes all the NN sites of $i$. Since Eq. (II.1) is a classical analogue of the Bloch equation, namely, the Heisenberg equation for the spin operator, all the static and dynamical magnetic properties purely intrinsic to the Hamiltonian (1) should be described by the combined use of Eqs. (1) and (II.1). A familiar alternative way to examine the spin dynamics is solving the Landau-Lifshitz-Gilbert (LLG) equation LLG_Landau_35 which includes a damping term originally introduced phenomenologically. In this work, we use Eq. (II.1) instead of the LLG equation for the following two reasons: (i) In the LLG equation, the damping, which is characterized by a dimensionless parameter $\alpha$, may be either intrinsic or extrinsic to the spin Hamiltonian. Equation (II.1), on the other hand, corresponds to the LLG equation without the phenomenological damping term, so that any relaxation described by Eq. (II.1) has its origin in the Hamiltonian (1). As our focus in the present paper is on fundamental aspects intrinsic to the spin Hamiltonian (1), we use Eq. (II.1); (ii) As we will see in the following subsection, the conventional forms of the spin and thermal currents SpinDyn_Huber_74 ; SpinDyn_Jencic_prb_15 ; MHall_Mook_prb_16 ; MHall_Mook_prb_17 ; Thermal_Huber_ptp_68 ; SpinDyn_Zotos_prb_05 ; SpinDyn_Sentef_07 ; SpinDyn_Pires_09 ; SpinDyn_Chen_13 are derived from the Heisenberg equation or its classical analogue without the damping term, so that it is self-consistent to use Eq. (II.1) rather than the LLG equation with the additional damping term.

II.2 Conductivities of spin and thermal currents

In this subsection, we will derive the spin current ${\bf J}_{s}^{z}$ and the thermal current ${\bf J}_{th}$, and then, introduce the spin-current conductivity $\sigma^{s}_{\mu\nu}$ and the thermal conductivity $\kappa_{\mu\nu}$. We shall start from the general discussion on a current in the continuum limit. Suppose that a conserved physical quantity ${\cal O}=\int d{\bf r}\,{\cal O}({\bf r},t)$ should satisfy the continuity equation $\frac{\partial}{\partial t}{\cal O}({\bf r},t)+\nabla\cdot{\bf j}_{\cal O}({\bf r},t)=0$ with associated local current density ${\bf j}_{\cal O}({\bf r},t)$. By multiplying the both side of the equation by ${\bf r}$ and integrating over the whole ${\bf r}$ region, we obtain

[TABLE]

Thus, the net current ${\bf J}_{\cal O}(t)$ is given by book_Mahan

[TABLE]

In the present XXZ model given by Eq. (1), the conserved quantities are the $z$ component of the magnetization $M^{z}=\sum_{i}S^{z}_{i}$ and the total energy ${\cal H}=\sum_{i}{\cal H}_{i}$ with {\cal H}_{i}=\frac{-J}{2}\sum_{j\in N(i)}\big{(}S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j}+\Delta S^{z}_{i}S^{z}_{j}\big{)}, so that the associated currents, namely, the spin and thermal currents (${\bf J}^{z}_{s}$ and ${\bf J}_{th}$) are given by

[TABLE]

[TABLE]

where Eq. (II.1) has been used in going from the first line to the second line for each current. The obtained result is essentially the same as the previously obtained expressions SpinDyn_Huber_74 ; SpinDyn_Jencic_prb_15 ; MHall_Mook_prb_16 ; MHall_Mook_prb_17 ; Thermal_Huber_ptp_68 ; SpinDyn_Zotos_prb_05 ; SpinDyn_Sentef_07 ; SpinDyn_Pires_09 ; SpinDyn_Chen_13 . We note that Eqs. (5) and (6) for the classical spin systems can also be applied for quantum spin systems by merely replacing ${S}_{i}^{\alpha}$ with the associated spin operator $\hat{S}_{i}^{\alpha}$. Indeed, one can verify that with the use of the Heisenberg equation instead of Eq. (II.1), the same expressions as Eqs. (5) and (6) are obtained.

Next, we turn to the conductivities of the spin and thermal currents. We first introduce the theoretical framework for the quantum mechanical systems, and then, take the classical limit. In general, driving forces for the spin and thermal currents are magnetic-field and temperature gradients, $\nabla H$ and $\nabla T$, respectively [see Figs. 1 (a) and (b)], so that the linear response equations are given by

[TABLE]

with the spin and thermal current-densities ${\bf j}_{s}$ and ${\bf j}_{th}$ MHall_Mook_prb_16 ; MHall_Mook_prb_17 ; book_Mahan . Then, the spin-current conductivity $\sigma^{s}$ and the thermal conductivity $\kappa$ are expressed as

[TABLE]

Note that in the present model without a magnetic field, $L^{s,th}=L^{th,s}=0$ is satisfied because these quantities are odd with respect to spins. In the linear response theory KuboFormular_Kubo_57 , the coefficients $L^{a,b}$ can be calculated from the formula

[TABLE]

where $\langle...\rangle$ denotes the thermal average in the equilibrium state. Now, we will take the classical limit of Eq. (9). In the classical system, by making $\hbar\rightarrow 0$ KuboFormular_Kubo_57 , we have

[TABLE]

Thus, in the present classical XXZ model, we obtain the following expressions for the spin-current and thermal conductivities:

[TABLE]

where we have used the relation between the total current and its current density, ${\bf j}_{s}={\bf J}^{z}_{s}/L$ and ${\bf j}_{th}={\bf J}_{th}/L$, with $L$ being a linear system size SpinDyn_Zotos_prb_05 ; SpinDyn_Jencic_prb_15 ; MHall_Mook_prb_16 ; MHall_Mook_prb_17 . Now, the problem is reduced to calculate the time correlations of the spin and thermal currents at various temperatures. In the present square lattice, the total number of spin $N_{\rm spin}$ and the system size $L$ are related by $L^{2}=N_{\rm spin}\,a^{2}$, where $a$ is a lattice constant. Noting that the time $t$ is measured in units of $|J|^{-1}$, it turns out that $\sigma^{s}_{\mu\nu}$ is a dimensionless quantity and $\kappa_{\mu\nu}$ has the dimension of $|J|$. Although in Eqs. (5) and (6), the currents themselves involve the dimension of length, the conductivities in the present two-dimensional system do not, so that the length scale of the lattice constant $a$ is not relevant and thus, we take $a=1$ throughout this paper except for the case where $a$ is explicitly written.

II.3 Numerical method

The time evolutions of ${\bf J}^{z}_{s}$ and ${\bf J}_{th}$ are determined microscopically by the spin-dynamics equation (II.1), so that we numerically integrate Eq. (II.1) and calculate the time correlations $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(t)\rangle$ and $\langle J_{th,\nu}(0)\,J_{th,\mu}(t)\rangle$ at each time step. In the numerical integration of Eq. (II.1), we use the second order symplectic method which guarantees the exact energy conservation Symplectic_Krech_98 ; Sqomega_Okubo_jpsj_10 ; Symplectic_Furuya_11 . We have confirmed that numerical results shown below are not altered if the 4th order Runge-Kutta method is used instead of the symplectic method. To properly evaluate the integral over time in Eq. (II.2), we perform long-time integrations typically up to $t=100\,|J|^{-1}\,-\,800\,|J|^{-1}$ with the time step $\delta t=0.01\,|J|^{-1}$ until the time correlations \big{\langle}J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(t)\big{\rangle} and \big{\langle}J_{th,\nu}(0)\,J_{th,\mu}(t)\big{\rangle} are completely lost.

Since Eq. (II.1) does not have a phenomenological dissipation term, the thermal fluctuations are only one possible cause for the current relaxation. Although Eq. (II.1) itself is deterministic, such a temperature effect can be incorporated by using temperature-dependent equilibrium spin configurations as the initial states for the equation of motion (II.1). In order to thermalize the system to given temperature $T$, we perform MC simulations for the spin Hamiltonian (1). The thermal average is taken as the average over initial equilibrium spin configurations generated in the MC simulations. In this work, at each temperature, we prepared 2000-4000 equilibrium spin configurations by picking up a spin snapshot in every 1000 MC sweeps after 105 MC sweeps for thermalization, where one MC sweep consists of the 1 heat-bath sweep and successive 10-30 over-relaxation sweeps.

By carefully analyzing the system-size dependences of the spin-current conductivity $\sigma^{s}_{\mu\nu}$ and the thermal conductivity $\kappa_{\mu\nu}$ at given temperatures, we will discuss the temperature dependences of $\sigma^{s}_{\mu\nu}$ and $\kappa_{\mu\nu}$ in the thermodynamic limit ($L\rightarrow\infty$) of our interest.

III Analytical results in the low-temperature limit: calculations based on the linear spin-wave theory

Before discussing numerical results, we should know how $\kappa_{\mu\nu}$ and $\sigma^{s}_{\mu\nu}$ should behave in the low-temperature limit. In this section, we will analytically investigate the temperature dependences of $\kappa_{\mu\nu}$ and $\sigma^{s}_{\mu\nu}$ based on the linear spin-wave theory (LSWT). As a low-temperature ordered state is a starting point in LSWT, one might be afraid that LSWT cannot be applied to the Heisenberg case because of the absence of the long-range order at any finite temperature. As long as there is a long-range order at $T=0$, however, the spin-wave expansions could still be done locally within the regions smaller than the spin-correlation length $\xi_{s}$ MagnonDamping_Tyc_89 . Thus, in the Heisenberg case, we introduce a lower cutoff in the momentum space which corresponds to the inverse spin-correlation length $\xi_{s}^{-1}$, and take the temperature dependence of $\xi_{s}\sim a\exp[b_{H}|J|/T]$ into account, where $b_{H}\simeq 2\pi$ is a universal constant Heisenberg_Polyakov_75 .

In this section, we will start from the theory of the corresponding quantum spin system, and then, take the classical limit of relevant physical quantities. By performing the spin-wave expansion, one can obtain the magnon representation of the Hamiltoninan (1) and the spin and thermal currents in Eqs. (5) and (6). Since in Eq. (II.2), the time correlation functions $\langle J_{th,\nu}(0)\,J_{th,\mu}(t)\rangle$ and $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(t)\rangle$ are essential for $\kappa_{\mu\nu}$ and $\sigma^{s}_{\mu\nu}$, we will first examine the associated static thermodynamic quantities, the equal-time correlation function $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle$ and $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(0)\rangle$. Then, the dynamical quantities, i.e., $\kappa_{\mu\nu}$ and $\sigma^{s}_{\mu\nu}$ due to the magnon propagation, will be calculated, putting emphasis on their temperature dependences in the classical limit. As we will see below, the temperature dependence of the thermal conductivity $\kappa_{\mu\nu}$ is almost independent of the magnetic anisotropy $\Delta$, while the spin-current conductivity $\sigma^{s}_{\mu\nu}$ is sensitive to the ordering properties controlled by $\Delta$.

III.1 Magnon representation

Although our target system in the present paper is the classical XXZ model, we consider, for convenience, the corresponding quantum spin system throughout this subsection. The magnon representation of the Hamiltonian (1) and the spin and thermal currents in Eqs. (5) and (6) can be derived by using the spin-wave expansions. In the Ising case of $\Delta>1$, the quantization axis of spin is in the $z$ direction, so that we introduce the transformation from the laboratory frame to the rotated frame with $y$ being the rotation axis,

[TABLE]

where $\theta_{i}={\bf Q}\cdot{\bf r}_{i}$ and ${\bf Q}=(\pi,\pi)$ is the ordering vector of the two-sublattice antiferromagnetic order. Then, the Hamiltonian reads

[TABLE]

By using the Holstein-Primakoff transformation

[TABLE]

with $\hat{a}^{\dagger}_{i}$ and $\hat{a}_{i}$ being respectively the bosonic creation and annihilation operators and the Fourier transformation of these operators

[TABLE]

we obtain

[TABLE]

where $A_{\bf q}=-4JS\Delta$, $B_{\bf q}=-4JS\gamma_{\bf q}$, and \gamma_{\bf q}=\frac{1}{2}\big{[}\cos(q_{x})+\cos(q_{y})\big{]}. The above Hamiltonian for the $\hat{a}_{\bf q}$ magnons can be diagonalized with the help of the Bogoliubov transformation

[TABLE]

where $\hat{b}^{\dagger}_{\bf q}$ and $\hat{b}_{\bf q}$ are the creation and annihilation operators for magnons. In the $XY$ ($\Delta<1$) and the Heisenberg ($\Delta=1$) cases, we take the quantization axis in the $x$ and $z$ directions, respectively. The diagonalized magnon Hamiltonian in the three cases, $\Delta>1$, $\Delta=1$, and $\Delta<1$, is summarized as follows:

[TABLE]

where we have dropped constant and higher-order terms. Note that in the $XY$ and Heisenberg cases of $\Delta\leq 1$, the magnon is a gapless excitation, while in the Ising case of $\Delta>1$, the magnon excitation has the gap $\Delta_{gp}=4|J|S\sqrt{\Delta^{2}-1}$.

In the same manner, the thermal and spin currents in Eqs. (6) and (5) can be expressed by the $\hat{b}_{\bf q}$ magnons as follows:

[TABLE]

[TABLE]

where

[TABLE]

Since ${\bf v}_{\bf q}=\nabla_{\bf q}\varepsilon_{\bf q}$ represents the magnon velocity, the thermal current ${\bf J}_{th}$ can be regarded as the energy flow carried by the magnons. In contrast to the thermal current ${\bf J}_{th}$ having the common magnon-representation independent of the magnetic anisotropy $\Delta$, the spin current in the Ising and Heisenberg cases ($\Delta\geq 1$) is expressed in the form fundamentally different from the one in the $XY$ case ($\Delta<1$). The former has the leading order contribution of the order of {\cal O}\big{(}S^{1}\big{)}, while the latter does not. In the $XY$ case, the spin current due to the magnon propagation is of the order of {\cal O}\big{(}S^{1/2}\big{)}. As the spin-wave expansion is the expansion with respect to $1/S$, such a higher order term is dropped in spirits of LSWT, so that ${\bf J}^{z}_{s}$ vanishes in the low-temperature ordered phase of the $XY$-type spin systems. The difference between $\Delta\geq 1$ and $\Delta<1$ cases stems from the difference in the direction of the quantization axis of spin: for $\Delta\geq 1$, the quantization axis is in the $z$ direction, whereas for $\Delta<1$, it is in the $xy$ plane which is perpendicular to the spin polarization of the spin current ${\bf J}^{z}_{s}$. Remember that although the spin current has its foundation on the conservation of the magnetization, only the $z$ component of the magnetization is conserved in the present anisotropic XXZ model in Eq. (1).

III.2 Static physical quantities

As the magnon Hamiltonian (III.1) is already diagonalized, one can easily calculate the thermal average of the current-related static quantities, $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle$ and $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(0)\rangle$.

First, we consider the thermal average of the equal-time correlation function for the thermal current, $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle$. With the use of the magnon representation in Eq. (22), we have

[TABLE]

where we have used the formula

[TABLE]

with the Bose-Einstein distribution function $f_{\rm B}(x)=(e^{x/T}-1)^{-1}$. Note that in Eq. (III.2), the off-diagonal term of $\mu\neq\nu$ vanishes after the summation over ${\bf q}$ because $v_{{\bf q},\mu}\propto\sin(q_{\mu})$ is an odd function of ${\bf q}$.

Now, we shall move on to the classical spin system. In the classical limit of

[TABLE]

the equal-time correlation for the classical spins $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle_{\rm cl}$ is obtained as

[TABLE]

At this point, the $T^{2}$ dependence of $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle_{\rm cl}$ is clear at least in the Ising and $XY$ cases. In the Heisenberg case, however, the additional temperature dependence due to the spin-correlation length $\xi_{s}$ comes in through the summation over ${\bf q}$. As we mentioned in the beginning of this section, $\xi_{s}$ enters in the form of the lower cutoff in the ${\bf q}$ space, i.e., $\xi_{s}^{-1}\leq|{\bf q}|$. For completeness, we shall evaluate the summation over ${\bf q}$ in Eq. (28) in all the three cases. Since the dominant contribution comes from the low-energy excitation near $|{\bf q}|\simeq 0$, we have

[TABLE]

Then, the ${\bf q}$-summation can be replaced with the following integral over $\varepsilon_{\bf q}$,

[TABLE]

where the density of states $D(\varepsilon_{\bf q})$ and the higher energy cutoff $\varepsilon_{\rm max}$ are given by $D(\varepsilon_{\bf q})=2\varepsilon_{\bf q}/(4|J|S)^{2}$ and $\varepsilon_{\rm max}\sim 4|J|S\Delta$ for $\Delta\geq 1$, and $D(\varepsilon_{\bf q})=[4/(1+\Delta)]\varepsilon_{\bf q}/(4|J|S)^{2}$ and $\varepsilon_{\rm max}\sim 4|J|S$ for $\Delta<1$. Note that in the Heisenberg case of $\Delta=1$, the low-energy cutoff $\varepsilon_{\rm min}$ possesses the temperature dependence via the spin-correlation length $\xi_{s}/a\sim\exp[b_{H}|J|/T]$. As we will see below, this additional temperature dependence coming from $\xi_{s}$ is negligibly small for the thermal transport, but not for the spin transport. By using Eq. (III.2) and performing the integral over $\varepsilon_{\bf q}$, we can evaluate the ${\bf q}$-summation in Eq. (28) to yield

[TABLE]

As the correction $(a/\xi_{s})^{2}$ in the $\Delta=1$ case is negligibly small, $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle$ in the classical limit exhibits the $T^{2}$ behavior at low temperatures, being independent of the magnetic anisotropy $\Delta$.

Next, we calculate the equal-time correlation function for the spin current in the classical limit, $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(0)\rangle_{\rm cl}$. Since in the $XY$ case of $\Delta<1$, the spin current is absent within the leading-order magnon contribution [see Eq. (23)], we only consider the $\Delta\geq 1$ case in which after some manipulations, we have

[TABLE]

Note that Eq. (III.2) is obtained for the quantum spin system. Now, we take the classical limit of Eq. (III.2). As the relations, $A_{\pm{\bf q}+{\bf Q}}=A_{\bf q}$, $B_{\pm{\bf q}+{\bf Q}}=-B_{\bf q}$, and ${\bf v}_{\pm{\bf q}+{\bf Q}}=\pm{\bf v}_{\bf q}$, are satisfied for $\Delta\geq 1$, the classical limit Eq. (27) yields

[TABLE]

By using Eq. (III.2), one can evaluate the summation over ${\bf q}$ in Eq. (38). The final result is summarized as follows:

[TABLE]

Note that in the $XY$ case of $\Delta<1$, \big{\langle}J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(0)\big{\rangle}_{\rm cl} is zero because the spin current is absent within the leading-order magnon contribution [see Eq. (23)]. In the Ising case of $\Delta>1$, the equal-time correlation of the spin current $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(0)\rangle_{\rm cl}$ has the same $T^{2}$ dependence as $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle_{\rm cl}$. In the Heisenberg case of $\Delta=1$, on the other hand, $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(0)\rangle_{\rm cl}$ includes a non-negligible correction term coming from the temperature-dependent $\xi_{s}$, i.e., $T^{2}\ln(\xi_{s}/a)\sim b_{H}|J|\,T$, and thus, takes the form of $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(0)\rangle_{\rm cl}/L^{2}\sim\delta_{\mu\nu}(const\,T^{2}+T)$. The correction term ($\propto T$) becomes the leading order contribution at lower temperatures, which is in sharp contrast to $\langle J_{th,\nu}(0)\,J_{th,\mu}(0)\rangle_{\rm cl}$ with the irrelevant correction terms [see Eq. (III.2)]. As we will see below, such a situation is also the case for the dynamical quantities.

III.3 Dynamical physical quantities

In the classical spin systems, the conductivities $\kappa_{\mu\nu}$ and $\sigma^{s}_{\mu\nu}$ are obtained from the time-correlation of the associated currents [see Eq. (II.2)]. Here, we consider the current dynamics brought by the magnon propagation in the presence of the magnon-magnon scatterings. In order to calculate the thermal average of the time correlation, it is convenient to start from the quantum mechanical system and take the classical limit of Eq. (27) afterwards. In the quantum mechanical system, the dynamical correlation function $L^{a,b}_{\mu\nu}(\omega)$ in Eq. (9) can be expressed in the following form book_AGD :

[TABLE]

Here, $Q^{a}_{\mu\nu}(i\omega_{l})$ is a response function and $\omega_{n}=2\pi nT$ is the bosonic Matsubara frequency. Then, the thermal conductivity $\kappa_{\mu\nu}$ and the spin-current conductivity $\sigma^{s}_{\mu\nu}$ are given by

[TABLE]

We first calculate the thermal conductivity $\kappa_{\mu\nu}$. For the thermal current carried by the magnons in Eq. (22), the response function $Q^{th}_{\mu\nu}(i\omega_{n})$ is given by book_AGD

[TABLE]

where ${\cal D}^{R}_{\bf q}(x)$ ({\cal D}^{A}_{\bf q}(x)=\big{[}{\cal D}^{R}_{\bf q}(x)\big{]}^{*}) is the retarded (advanced) magnon Green’s function obtained by analytic continuation $i\omega_{m}\rightarrow\omega+i0$ in the temperature Green’s function ${\cal D}_{\bf q}(i\omega_{m})$ defined by

[TABLE]

With the use of Eq. (III.3), the thermal conductivity in the quantum system is formally expressed as

[TABLE]

Here, the magnon Green’s function ${\cal D}_{\bf q}^{R}(x)$ is given by

[TABLE]

where the dimensionless coefficient $\alpha$ represents the magnon damping which corresponds to the Gilbert damping in the LLG equation MagnonGreen_Yamaguchi_17 ; MagnonTrans_Tatara_15 . In general, the damping $\alpha$ originates from the interactions associated with spins in solids, so that it may be brought not only by the magnon-magnon scatterings but also, for example, by magnon-phonon scatterings. In the present work, however, the starting point is the spin Hamiltonian (1) and no further assumption is made. Thus, $\alpha$ is of purely magnetic origin and brought by the magnon-magnon scatterings. Since the temperature dependence of $\alpha$ has already been calculated in the typical case of $\Delta=1$ MagnonDamping_Tyc_89 ; MagnonDamping_Harris_71 , we will skip the microscopic derivation of $\alpha$ in this paper.

In the classical spin system, the concrete expression of Eq. (47) can straightforwardly be derived, as shown below. Substituting Eq. (48) into Eq. (47) and taking the classical limit of $f_{\rm B}^{\prime}(x)=-T/x^{2}$, we obtain the following expression for the thermal conductivity in the classical spin systems $\kappa_{\mu\nu}^{\rm cl}$ as

[TABLE]

where the equation

[TABLE]

has been used. The summation over ${\bf q}$ can be evaluated in the same manner as that for the static physical quantities. With the use of Eq. (III.2), we obtain

[TABLE]

Only the longitudinal components of the thermal conductivity $\kappa_{\mu\mu}^{\rm cl}$ are non-vanishing. When the magnon damping is sufficiently small such that $\alpha\ll 1$, it follows that $\kappa_{\mu\nu}^{\rm cl}\propto 1/\alpha$, which agrees with the results obtained in other theoretical approaches MagnonTrans_Tatara_15 ; MagnonTrans_Jiang_13 .

One can see from Eq. (III.3) that in the Heisenberg case of $\Delta=1$, although the spin-correlation length $\xi_{s}$ rapidly increases toward $T=0$, such a temperature effect is irrelevant at lower temperatures because $\xi_{s}$ enters in $\kappa_{\mu\mu}^{\rm cl}$ in the form of $1/\xi_{s}$. Thus, in all the three ($\Delta>1$, $\Delta=1$, and $\Delta<1$) cases, the temperature dependence of $\kappa_{\mu\mu}^{\rm cl}\propto 1/\alpha$ is governed by the magnon damping factor $\alpha$.

The damping of the antiferromagnetic magnon due to multi-magnon scatterings has already been calculated by using Feynman diagram techniques in Refs. MagnonDamping_Tyc_89 ; MagnonDamping_Harris_71 . The temperature dependence of $\alpha$ in the classical Heisenberg antiferromagnet essentially follows the $T^{2}$ form, i.e., $\alpha\propto T^{2}$, which results from the leading-order scattering process involving four magnons. In the $XY$-type and Ising-type classical spin systems, although the concrete expression of $\alpha$ is not available, the same temperature dependence $\alpha\propto T^{2}$ is expected because the same types of the Feynman diagrams (the same leading-order scattering processes) contribute to the magnon damping. Of course, there must be quantitative differences among the three cases. In particular, for the Ising-type anisotropy of $\Delta>1$, the magnon excitation is gapped, so that the phase space satisfying the energy conservation in the calculation of the relevant Feynman diagrams would be shrunk with increasing $\Delta$, resulting in a smaller value of $\alpha$. Apart from such a quantitative difference which may become serious for strong Ising-type anisotropies, the longitudinal thermal conductivity $\kappa_{\mu\mu}^{\rm cl}$ in the classical limit should behave as $\kappa_{\mu\mu}^{\rm cl}\propto 1/\alpha\propto 1/T^{2}$ in all the three ($\Delta>1$, $\Delta=1$, and $\Delta<1$) cases.

Now, we will move on to the calculation of the spin-current conductivity $\sigma^{s}_{\mu\nu}$ based on Eq. (III.3). As in the case of the thermal current, starting from the magnon representation of the spin current in Eq. (23), we can write down the response function $Q^{s}_{\mu\nu}(i\omega_{n})$ as

[TABLE]

Then, the spin-current conductivity $\sigma^{s}_{\mu\nu}$ is formally written as

[TABLE]

In the same manner as that for $\kappa_{\mu\nu}$, we take the classical limit of Eq. (III.3). By substituting Eq. (48) into Eq. (III.3), taking the classical limit of $f_{\rm B}^{\prime}(x)=-T/x^{2}$, and using Eq. (50) and the formula

[TABLE]

we have the spin-current conductivity in the classical spin systems $\sigma^{s,{\rm cl}}_{\mu\nu}$ as follows:

[TABLE]

By further using the approximation Eq. (III.2), we finally obtain

[TABLE]

In contrast to the thermal conductivity $\kappa_{\mu\nu}^{\rm cl}$, the spin-current conductivity $\sigma^{s,{\rm cl}}_{\mu\nu}$ reflects the difference in the ordering properties. First of all, in the $XY$ case of $\Delta<1$, $\sigma^{s,{\rm cl}}_{\mu\nu}$ is zero because the spin current is absent within the leading-order magnon contribution [see Eq. (23)]. In the Ising case of $\Delta>1$, as one can see from Eq. (58), the temperature dependence of $\sigma^{s,{\rm cl}}_{\mu\mu}$ is determined by that of $T/\alpha$. Since for relatively weak anisotropies, $\alpha\propto T^{2}$ is expected to be satisfied, the longitudinal spin-current conductivity should exhibit the following temperature dependence: $\sigma^{s,{\rm cl}}_{\mu\mu}\propto T/\alpha\propto T^{-1}$. In the Heisenberg case of $\Delta=1$, one can see from Eq. (58) that the spin-correlation length $\xi_{s}$ enters in the form of $\xi_{s}/\alpha$, so that the longitudinal spin-current conductivity should diverge toward $T=0$ in the exponential form of $\sigma^{s,{\rm cl}}_{\mu\mu}\propto\xi_{s}\,T/\alpha\sim\exp[b_{H}|J|/T]$.

In the following sections, we will show numerical results on $\kappa_{\mu\nu}$ and $\sigma^{s}_{\mu\nu}$, the low-temperature properties of which are qualitatively consistent with the above analytical results. It should be noted that the transport properties near the phase transition, which is our main focus of the present work, is out of the applicability range of LSWT.

IV Numerical results on the thermal conductivity

In this section, we will discuss the association between the phase transition and the thermal transport based on numerical results obtained in the Ising-type ($\Delta>1$), $XY$-type ($\Delta<1$), and Heisenberg-type ($\Delta=1$) spin systems. In this paper, the parameter values of $\Delta=1.05$ and $\Delta=0.95$ are basically used in the Ising and $XY$ cases, respectively, as typical values slightly deviating from $\Delta=1$ of the isotropic Heisenberg case. From the MC simulations (see Appendix), the transition temperature in each case is estimated to be $T_{N}/|J|\simeq 0.75$ for $\Delta=1.05$ and $T_{KT}/|J|\simeq 0.6$ for $\Delta=0.95$ KT_XXZ_Cuccoli_95 ; KT_XXZ_Lee_05 ; KT_XXZ_Pires_96 .

In Eq. (II.2), the temperature dependence of $\kappa_{\mu\nu}$ is determined by the integrated value of the time correlation of the thermal current $\langle J_{th,\nu}(0)\,J_{th,\mu}(t)\rangle$ except the trivial $T^{-2}$ factor, so that we will start from the temperature dependence of $\langle J_{th,\nu}(0)\,J_{th,\mu}(t)\rangle$. Figure 2 shows the time correlation function normalized by the system size $\langle j_{th,x}(0)\,j_{th,x}(t)\rangle\equiv\langle J_{th,x}(0)\,J_{th,x}(t)\rangle/L^{2}$ at different temperatures in the Ising-type ($\Delta=1.05$), $XY$-type ($\Delta=0.95$), and Heisenberg-type ($\Delta=1$) spin systems. System-size dependence can hardly be seen, suggesting that the thermal transport is a spatially local phenomenon. As for the effect of the magnetic anisotropy, there is no qualitative difference among the three cases. With decreasing temperature, the time correlation decays more slowly in time. In other words, the relaxation time of the thermal current, which we denote as $\tau_{th}$, becomes longer. Thus, the associated thermal conductivity $\kappa_{\mu\nu}$ is expected to follow a common monotonic temperature-dependence.

Figure 3 shows the longitudinal and transverse thermal conductivities as a function of temperature $T$ in the Ising-type ($\Delta=1.05$), $XY$-type ($\Delta=0.95$), and Heisenberg-type ($\Delta=1$) spin systems. Because the $yy$ ($xy$) component of $\kappa_{\mu\nu}$ is equivalent to the $xx$ ($yx$) component in the present square-lattice NN model, only the the $xx$ and $yx$ components, $\kappa_{xx}$ and $\kappa_{yx}$, are shown. One can see from Fig. 3 that in all the three cases, the transverse Hall response $\kappa_{yx}$ is absent at $2\sigma$ precision (see lower panels) and the longitudinal thermal conductivity $\kappa_{xx}$ gradually increases toward $T=0$ (see the upper main panels). Although the phase transition occurs in the anisotropic spin systems, no clear anomaly can be seen in the thermal conductivity at the magnetic transition temperature $T_{N}$ or the KT topological transition temperature $T_{KT}$. Thus, in view of the main focus of this work, our conclusion is that the strong association between the thermal conductivity and the phase transition cannot be observed in the present NN XXZ model in two dimensions. Below in this section, to shed light on the basic properties of the thermal transport in the classical spin systems, we will devote ourselves to the low-temperature behavior of the longitudinal thermal conductivity $\kappa_{xx}$.

For the $XY$-type anisotropy $\Delta<1$, the temperature dependence of $\kappa_{xx}$ in Fig. 3 (b) is not altered qualitatively by the change in $\Delta$. For the Ising-type anisotropy $\Delta>1$, on the other hand, the magnon excitation has the gap $\Delta_{gp}=4|J|S\sqrt{\Delta^{2}-1}$, so that the thermal current, which is the energy flow carried by the manons, and the associated conductivity $\kappa_{xx}$ are expected to be suppressed with increasing $\Delta$. Figure 4 shows the longitudinal thermal conductivity $\kappa_{xx}$ as a function of $T/T_{N}$ for various values of $\Delta>1$. Not only the absolute value of $\kappa_{xx}$ but also the divergent behavior toward $T=0$ is suppressed by the increase of $\Delta$. At least for not so strong Ising-type anisotropy, however, $\kappa_{xx}$ tends to diverge toward $T=0$, roughly showing a power-law behavior. Hereafter, we will discuss the origin of such a power-law-type temperature dependence, focusing on the almost isotropic spin systems.

As one can see from Eq. (II.2), $\kappa_{xx}$ involves the trivial $T^{-2}$ dependence. In order to extract the nontrivial temperature dependence other than the $T^{-2}$ factor, $T^{2}\kappa_{xx}=\int dt\,\langle j_{th,x}(0)\,j_{th,x}(t)\rangle$ is plotted in the insets of the upper panels of Fig. 3 as a function of temperature. In the anisotropic cases of $\Delta\neq 1$, $T^{2}\kappa_{xx}$ tends to saturate to a constant value at the lowest temperature, whereas in the isotropic case of $\Delta=1$, it remains increasing toward $T=0$. Except this difference at the lowest temperature, $T^{2}\kappa_{xx}$ shows a weak monotonic increase below $T/|J|\leq 0.8$ in both the anisotropic and isotropic cases. Thus, the divergent behavior toward $T=0$ in $\kappa_{xx}$ is mainly due to the $T^{-2}$ factor, but in the low-temperature range of our simulations, $\kappa_{xx}$ increases slightly faster than $T^{-2}$ due to the non-trivial contribution originating from the thermal fluctuation, $T^{2}\kappa_{xx}=\int dt\,\langle j_{th,x}(0)\,j_{th,x}(t)\rangle$. The analytical result in Eq. (III.3), on the other hand, shows that the thermal conductivity due to the magnon propagation should behave as $\kappa_{xx}\propto 1/\alpha\propto T^{-2}$. As mentioned above, at least in the temperature range of our simulations, the numerically obtained $\kappa_{xx}$ increases faster than $T^{-2}$. In order to examine the origin of the deviation between the numerical and analytical results on the temperature dependence of $\kappa_{xx}$, we will look into the details of the temperature dependences of the physical quantities related to $\langle j_{th,x}(0)\,j_{th,x}(t)\rangle$.

In Fig. 2, the time correlation $\langle j_{th,x}(0)\,j_{th,x}(t)\rangle$ decays exponentially in the form of $e^{-t/\tau_{th}}$ with the relaxation time of the thermal current $\tau_{th}$, so that we could assume $\langle j_{th,x}(0)\,j_{th,x}(t)\rangle\simeq\langle|j_{th,x}(0)|^{2}\rangle e^{-t/\tau_{th}}$. Then, by carrying out the integral over time in Eq. (II.2), one can estimate the longitudinal thermal conductivity as $\kappa_{xx}\simeq T^{-2}\,\langle|j_{th,x}(0)|^{2}\rangle\,\tau_{th}$. As the data on the static quantity $\langle|j_{th,x}(0)|^{2}\rangle$ can be compared directly with the analytical result given in Eq. (III.2), one can relate $\tau_{th}$ to the magnon damping $\alpha$ via Eq. (III.3). If the equal-time correlation $\langle|j_{th,x}(0)|^{2}\rangle$ follows the $T^{2}$ dependence expected in LSWT, the relaxation time of the thermal current $\tau_{th}$ corresponds to the inverse magnon-damping $1/\alpha$ which is roughly proportional to $T^{-2}$ in the lowest-order approximation MagnonDamping_Tyc_89 ; MagnonDamping_Harris_71 .

Figure 5 shows the temperature dependences of $\langle|j_{th,x}(0)|^{2}\rangle$ and $\tau_{th}$ in the three cases of $\Delta=1.05$, $\Delta=0.95$, and $\Delta=1$, where $\tau_{th}$ is extracted by fitting the $\langle j_{th,x}(0)\,j_{th,x}(t)\rangle$ curve with the exponential form of $e^{-t/\tau_{th}}$. Since $\langle|j_{th,x}(0)|^{2}\rangle$ exhibits a power-law behavior, we fit the low-temperature data with the functional form of $T^{x}$ and find $x=1.8\sim 1.9$. The resultant fitting function $T^{x}$ in each case is represented by a dashed curve together with the obtained value of $x$ in Fig. 5. The exponent $x\simeq 2$ for $\langle|j_{th,x}(0)|^{2}\rangle$ is in good agreement with the analytical result given in Eq. (III.2), so that the origin of the discrepancy in the temperature dependence of $\kappa_{xx}$ between the numerical and analytical results consists in the relaxation time $\tau_{th}$ which should satisfy the relation $\tau_{th}\propto 1/\alpha\propto T^{-2}$. As one can see from the right panels in Fig. 5, however, $\tau_{th}$ diverges toward $T=0$ slightly faster than $T^{-2}$. A rough estimation, which is done by fitting all the low-temperature data for $T/|J|\leq 0.6$ with the functional form of $T^{x}$, yields $\tau_{th}\propto T^{-2.5}$ in all the three cases. The deviation from the expected behavior $1/\alpha\propto T^{-2}$ may be attributed to the temperature range considered. The temperature range available for fitting might be higher than that assumed in the analytical calculation where higher-order multi-magnon-scattering processes are neglected. With further decreasing temperature below the lowest temperature of our simulation, $\tau_{th}$ and resultant $\kappa_{xx}$ should tend to obey the expected power-law form $T^{-2}$. Actually, in the Ising and XY cases, a precursor of such a tendency has already been observed as the saturated behavior in $T^{2}\kappa_{xx}$ (see the insets of Fig. 3).

V Numerical results on the spin-current conductivity

In Sec. III, based on the analytical calculations in LSWT, we find that the effect of the magnetic anisotropy $\Delta$, i.e., the difference in the ordering properties, is reflected in the low-temperature spin-transport. In this section, we will discuss the association between the phase transition and the spin-current conductivity $\sigma^{s}_{\mu\nu}$, based on numerical results.

We shall start from the time correlation function of the spin current $\langle J^{z}_{s,\nu}(0)\,J^{z}_{s,\mu}(t)\rangle$ which yields the nontrivial temperature dependence of $\sigma^{s}_{\mu\nu}$ [see Eq. (II.2)]. Figure 6 shows the time correlation function normalized by the system size $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle=\langle J^{z}_{s,x}(0)\,J^{z}_{s,x}(t)\rangle/L^{2}$ at various temperatures in the typical three cases, Ising-type ($\Delta=1.05$), $XY$-type ($\Delta=0.95$), and Heisenberg-type ($\Delta=1$) spin systems. These $\Delta$ values are the same as those in Figs. 2 and 3. At the high temperature $T/|J|=0.74$, one cannot see a clear difference among the three cases. With decreasing temperature, $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle$ exhibits characteristic behaviors depending on the ordering properties. In the Ising case of $\Delta=1.05$, $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle$ shows an oscillating behavior in the very-short time scale [see the insets of Fig. 6(a)], but its long-time relaxation whose characteristic time scale is denoted by $\tau_{s}$ becomes slower at lower temperatures without showing the system size dependence. In the $XY$ case of $\Delta=0.95$, the time correlation persists for a long time at $T/|J|=0.66$ slightly above $T_{KT}$, showing a large system size dependence, whereas the time correlation is lost within a short time scale at $T/|J|=0.3$ much lower than $T_{KT}$. In the Heisenberg case of $\Delta=1$, $\tau_{s}$ becomes longer with decreasing temperature like in the Ising case, but the system size dependence is quite large. The above difference is reflected in the spin-current conductivity $\sigma^{s}_{\mu\nu}$ through the integration of $\langle j^{z}_{s,\nu}(0)\,j^{z}_{s,\mu}(t)\rangle$ over the whole time range.

Figure 7 shows the temperature dependences of the longitudinal (upper panels) and transverse (lower panels) spin-current conductivities, $\sigma^{s}_{xx}$ and $\sigma^{s}_{yx}$, for $\Delta=1.05$ (a), $\Delta=1.05$ (b), and $\Delta=1$ (c). As one can see from Fig. 7, in all the three cases, the transverse Hall response $\sigma^{s}_{yx}\,(=\sigma^{s}_{xy})$ is absent also for the spin transport as well as the thermal transport. The longitudinal spin-current conductivity $\sigma^{s}_{xx}\,(=\sigma^{s}_{yy})$, on the other hand, exhibits temperature dependences characteristic of the three different universality classes. Here, we briefly summarize the temperature dependence of $\sigma^{s}_{xx}$, and a detailed analysis in each case will be given in the following subsections. In the Ising case of $\Delta=1.05$, $\sigma^{s}_{xx}$ gradually increases with decreasing temperature without showing a clear anomaly at the magnetic transition temperature $T_{N}$. Also, the system size dependence cannot be seen, as is already suggested from the size-independent time-correlation-functions in Fig. 6 (a). In the $XY$ case of $\Delta=0.95$, $\sigma^{s}_{xx}$ exhibits a divergent sharp peak toward the KT transition temperature $T_{KT}$, and becomes vanishingly small at lower temperatures below $T_{KT}$. In the Heisenberg case of $\Delta=1$, $\sigma^{s}_{xx}$ increases exponentially with decreasing temperature, showing a large system-size-dependence at lower temperatures. Below in this section, we will give a detailed description of the association between the longitudinal spin-current conductivity $\sigma^{s}_{xx}$ and the ordering properties of the system.

V.1 Ising-type spin system

In Fig. 7, for the Ising-type anisotropy of $\Delta=1.05$, a clear signature of the magnetic transition at $T_{N}$ cannot be seen in $\sigma^{s}_{xx}$. We will first check that this result is not altered qualitatively by the value of $\Delta$, and subsequently discuss the temperature dependence of $\sigma^{s}_{xx}$ in the long-range-ordered phase below $T_{N}$, making a comparison between the numerical result and the analytical one in Sec. III.

The gap-opening in the magnon excitation due to $\Delta$ is expected to suppress $\sigma^{s}_{xx}$, as is actually the case for the thermal conductivity $\kappa_{xx}$. Figure 8 shows $\sigma^{s}_{xx}$ as a function of $T/T_{N}$ for various values of $\Delta>1$. No clear signature of the magnetic transition can commonly be seen near $T_{N}$, and as is expected, $\sigma^{s}_{xx}$ is suppressed by the increase of $\Delta$. For relatively weak magnetic anisotropies, $\sigma^{s}_{xx}$ increases toward $T=0$ and its temperature dependence is almost compatible with the analytical expectation, $\sigma^{s}_{xx}\propto T/\alpha\propto T^{-1}$, given in Eq. (58). To look into the details of the temperature effect on $\sigma^{s}_{xx}$, we will examine the temperature dependence of the current-related quantities for $\Delta=1.05$.

In Fig. 6 (a), except for the short-time oscillating behavior, the time correlation $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle$ decays exponentially in the form of $e^{-t/\tau_{s}}$, so that we could roughly write $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle\sim\langle|j^{z}_{s,x}(0)|^{2}\rangle\,e^{-t/\tau_{s}}$. Then, from Eq. (II.2), the longitudinal spin-current conductivity $\sigma^{s}_{xx}$ can be evaluated as $\sigma^{s}_{xx}\sim T^{-1}\,\langle|j^{z}_{s,x}(0)|^{2}\rangle\,\tau_{s}$. If the static quantity $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ follows the $T^{2}$ behavior expected in LSWT [see Eq. (III.2)] as is actually the case for the thermal transport, it follows that $\sigma^{s}_{xx}\sim T\tau_{s}$. By comparing this expression to Eq. (58), one notice that $\tau_{s}$ is associated with the magnon damping $\alpha$ via $\tau_{s}\sim 1/\alpha$.

Figure 9 shows the temperature dependences of $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ and $\tau_{s}$, where $\tau_{s}$ is extracted by fitting the tail of $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle$ with $e^{-t/\tau_{s}}$. As one can see from the left panel of Fig. 9, $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ shows a power-law behavior of the form $T^{x}$ in the ordered phase, and the exponent $x$ is obtained by fitting the low-temperature data as $x=2$. The resultant fitting function is represented by a dashed curve in the left panel of Fig. 9. The obtained $T^{2}$ behavior for $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ is in good agreement with the analytical result in Eq. (III.2), so that $\tau_{s}\sim 1/\alpha\propto T^{-2}$ should be satisfied. The numerically obtained $\tau_{s}$ shown in the right panel of Fig. 9 tends to obey the expected power-law form $T^{-2}$, but in the wide low-temperature range of our simulation, it increases toward $T=0$ slightly faster than $T^{-2}$. When we fit all the low-temperature data below $T/|J|=0.6$ with the functional form $T^{x}$, the same temperature dependence as that of the thermal-current-relaxation time $\tau_{th}$ is obtained for the spin-current-relaxation time, namely, $\tau_{s}\propto T^{-2.5}$, indicating that in the Ising-type spin systems, the long-time relaxations of the spin and thermal transports are of the same origin, namely, the magnon damping due to the multi-magnon scatterings.

In the short-time scale, on the other hand, one can see the oscillating behavior in $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle$ [see the insets in Fig. 6 (a)], which is not observed in the thermal-current relaxation. Although the origin of the oscillation is not clear, this suggests that the spin-current relaxation may involve not only the ordinary magnon damping but also other effects of the magnetic excitations. As we will see below, in the $XY$-type spin systems, the vortex excitations come into play in the spin-current relaxation, leading to the divergence of $\sigma^{s}_{xx}$ at the KT transition temperature.

V.2 $XY$-type spin system

In the $XY$ antiferromagnet with the weak anisotropy $\Delta=0.95$, as shown in Fig. 7 (b), the longitudinal spin-current conductivity $\sigma^{s}_{xx}$ (=$\sigma^{s}_{yy}$) is significantly enhanced near $T_{KT}$, but once entering in the low-temperature phase below $T_{KT}$, $\sigma^{s}_{xx}$ becomes vanishingly small. These features are universal in the $XY$-type spin systems, being independent of the values of $\Delta$. Furthermore, even if the antiferromagnetic exchange interaction $J<0$ is replaced with a ferromagnetic one $J>0$, the universality class remains unchanged and the above features in $\sigma^{s}_{xx}$ can be observed. Figure 10 shows the temperature dependence of $\sigma^{s}_{xx}$ in the antiferromagnet ($J<0$) with $\Delta=0.2$ (a) and in the ferromagnet ($J>0$) with $\Delta=0$ (b). In both cases, a divergent sharp peak can clearly be seen near $T_{KT}$. With increasing the system size $L$, the peak height increases and the peak temperature approaches $T_{KT}$ from above, suggesting that in the thermodynamic limit of $L\rightarrow\infty$, $\sigma^{s}_{xx}$ diverges at $T_{KT}$. On crossing $T_{KT}$ from above, $\sigma^{s}_{xx}$ drops to a vanishingly small value. Hereafter, we will discuss the origin of this temperature dependence.

As discussed in Sec. III, in the ordered phase of the $XY$-type spin system, the leading-order magnon-spin-current is absent [see Eq. (23)] because of the orthogonal relation between the quantization axis lying in the $xy$-plane of the spin space and the polarization direction of the spin current which is in the $z$ direction in the present XXZ model. The associated spin-current conductivity $\sigma^{s}_{xx}$, therefore, should be vanishingly small, although higher-order magnon contributions may have a little effect on the spin transport. The low-temperature feature observed below $T_{KT}$ in Figs. 7 (b) and 10 is understood as a manifestation of this nature inherent to the $XY$-type anisotropy. Thus, the non-trivial issue is the significant enhancement of $\sigma^{s}_{xx}$ near $T_{KT}$ observed in the numerical simulations.

Since the system-size-dependent divergent peak near $T_{KT}$ is commonly observed for the $XY$-type anisotropy, we focus on the case of $\Delta=0.95$ as a representative example and discuss the thermodynamic limit ($L\rightarrow\infty$) of $\sigma^{s}_{xx}$. Figure 11 (a) shows the system-size dependence of $\sigma^{s}_{xx}$ at various temperatures. One can see that at temperatures away from $T_{KT}/|J|\simeq 0.6$, $\sigma^{s}_{xx}$ as a function of the system size $L$ saturates to a constant value, which corresponds to $\sigma^{s}_{xx}$ in the thermodynamic limit. The extrapolated $L\rightarrow\infty$ value of $\sigma^{s}_{xx}$ and the corresponding original finite-size data in Fig. 7 (b) are plotted in Fig. 11 (b) on the semi-logarithmic scale. The divergent behavior toward $T_{KT}/|J|\simeq 0.6$ and the sudden drop across $T_{KT}$ can clearly be seen. Noting that the spin correlation length $\xi_{s}$ in the KT transition is known to diverge in the form of \xi_{s}/a\sim\exp\big{[}b_{KT}/\sqrt{T/T_{KT}-1}\big{]} with $b_{KT}\simeq\pi/2$ KT_Kosterlitz_74 , we fit the $L\rightarrow\infty$ data of $\sigma^{s}_{xx}$ at $T\gtrsim T_{KT}$ with the functional form of A\exp\big{[}B/\sqrt{T/T_{KT}-1}\big{]}. The fitting parameters $A$ and $B$ are obtained as $B=2.26\pm 0.10$ and $A=0.008\pm 0.002$. The $\sigma^{s}_{xx}(T)$ curve extrapolated in this way is represented by a dashed curve in Fig. 11 (b). One can see that the obtained exponential form well characterizes the numerically-obtained divergent behavior of $\sigma^{s}_{xx}$, which, together with the obtained $B$-value comparable to $b_{KT}$, suggests that this pronounced spin-transport phenomenon is closely related to the KT transition, or equivalently, the vortex binding-unbinding process.

In the KT transition, the spin correlation length $\xi_{s}$ corresponds to the inter-free-vortex distance. With decreasing temperature above $T_{KT}$, the inter-free-vortex distance increases, so that it becomes difficult for a single vortex to find out a partner free anti-vortex to form a vortex pair. This means that in terms of the time evolution, the single free vortex wanders for a longer time until it collides with the partner free anti-vortex. Thus, the lifetime of the single free vortex should get longer on approaching $T_{KT}$ from above. Once across $T_{KT}$, all the vortices are paired up and a single vortex cannot be found any more. Bearing this fundamental physics of the KT transition in our mind, we examine the temperature dependences of $\tau_{s}$ and $\langle|j^{z}_{s,x}(0)|^{2}\rangle$.

Figure 12 shows the temperature dependences of $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ (a) and $\tau_{s}$ (b). The spin-current relaxation time $\tau_{s}$ is determined by fitting the long-time tail of $\langle j^{z}_{s,x}(0)j^{z}_{s,x}(t)\rangle$ in Fig. 8 (b) with the exponential form of $e^{-t/\tau_{s}}$. One can see from Fig. 12 that on approaching $T_{KT}$ from above, $\tau_{s}$ is significantly enhanced, while $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ only shows a slight increase. In the low-temperature phase below $T_{KT}$, $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ is strongly suppressed as is expected from the analytical result that the leading-order magnon-spin-current is absent, and correspondingly, the relaxation becomes so rapid that $\tau_{s}$ cannot be defined any more. The functional type characterizing the steep increase in $\tau_{s}$ is also the exponential one. By fitting the data at $T\gtrsim T_{KT}$ with the form of \tilde{A}\exp\big{[}\tilde{B}/\sqrt{T/T_{KT}-1}\big{]}, we obtain $\tilde{A}=0.013\pm 0.003$ and $\tilde{B}=2.58\pm 0.07$. The extrapolated $\tau_{s}(T)$ curve is represented by a dashed curve in Fig. 12. One can see that the obtained exponential form well characterizes the numerically-obtained divergent behavior of $\tau_{s}$. As $\sigma^{s}_{xx}$ is related to $\tau_{s}$ and $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ via $\sigma^{s}_{xx}\sim T^{-1}\,\langle|j^{z}_{s,x}(0)|^{2}\rangle\,\tau_{s}$, the divergent behavior in $\sigma^{s}_{xx}$ originates from the divergence of the spin-current-relaxation time $\tau_{s}$ toward $T_{KT}$. Actually, the obtained values of $B\simeq 2.26$ and $\tilde{B}\simeq 2.58$ almost coincide with each other.

Now, we will address the physical interpretation of the above result. In the KT topological transition, the distinct feature above $T_{KT}$ is the existence of an isolated free vortex and its dynamics. Since the vortex interacts with surrounding magnons or spin waves, the vortex motion is not ballistic, but rather diffusive KT-diffusive_Loft_87 ; KT-diffusive_Toyoki_90 ; KT-diffusive_Goldenfeld_90 ; KT-diffusive_Huse_93 ; KT-diffusive_Bray_00 . Thus, the vortex lifetime $\tau_{vtx}$ could be estimated roughly as \tau_{vtx}\propto\xi_{s}^{2}\sim\exp\big{[}2b_{KT}/\sqrt{T/T_{KT}-1}\big{]}, so that $\tau_{vtx}$ should get longer in the exponential form toward $T_{KT}$ with $2b_{XT}\simeq\pi$. Since the two time-scales, $\tau_{s}$ and $\tau_{vtx}$, develop toward $T_{KT}$ in the almost same manner as a function of temperature, it is naturally expected that the vortex excitations play an important role in the spin-current relaxation. Because $\sigma^{s}_{xx}$ is proportional to $\tau_{s}$, we could conclude that the divergent peak at $T_{KT}$ in the $\sigma^{s}_{xx}$ curve is attributed to the topological excitations of the long-life-time vortices.

V.3 Heisenberg-type spin system

In the Heisenberg case of $\Delta=1$, the spin space is isotropic, so that in contrast to the anisotropic cases of $\Delta\neq 1$, not only the $z$ component of the magnetization but also the $x$ and $y$ components are conserved quantities. This enables one to define the spin-currents ${\bf J}_{s}^{x}$ and ${\bf J}_{s}^{y}$ as well as ${\bf J}_{s}^{z}$, where {\bf J}_{s}^{\alpha}=J\sum_{\langle i,j\rangle}\big{(}{\bf r}_{i}-{\bf r}_{j}\big{)}\big{(}{\bf S}_{i}\times{\bf S}_{j}\big{)}^{\alpha} can be derived in the same manner as Eq. (5). Since all the spin currents, ${\bf J}_{s}^{x}$, ${\bf J}_{s}^{y}$, and ${\bf J}_{s}^{z}$, should be equivalent to one another, the associated spin-current conductivities should also be equivalent. Thus, in the Heisenberg case, we calculate the spin-current conductivity averaged over the three spin components

[TABLE]

instead of Eq. (II.2). The spin-current conductivity so obtained is shown in Fig. 7 (c) as a function of temperature. In the Heisenberg case, neither a magnetic transition nor a topological one does not occur, so that the characteristic temperature scale is absent except the exchange interaction $|J|$. In Fig. 7 (c), with decreasing temperature, the longitudinal spin-current conductivity $\sigma^{s}_{xx}$ increases monotonically and a steep increase sets in around $T/|J|\sim 0.8$. As the system size dependence of $\sigma^{s}_{xx}$ becomes considerably larger at lower temperatures, we will extrapolate the low-temperature $\sigma^{s}_{xx}(T)$ curve in the thermodynamic limit.

Figure 13 (a) shows the system-size dependence of $\sigma^{s}_{xx}$ at various temperatures. At lower temperatures, a larger system size is necessary to obtain the thermodynamic-limit value of $\sigma^{s}_{xx}$, suggesting that in contrast to the thermal transport which is a spatially local phenomenon, the spin transport captures the long-length-scale magnetic properties. The extrapolated thermodynamic-limit values of $\sigma^{s}_{xx}$ and the corresponding original finite-size data in Fig. 7 (c) are plotted in Fig. 13 (b) on the semi-logarithmic scale as a function of the inverse temperature $1/T$. As the $L\rightarrow\infty$ data at $T/|J|\lesssim 0.8$ are on a straight line, we fit them by the exponential function of A_{H}\exp\big{[}B_{H}|J|/T\big{]} with $A_{H}$ and $B_{H}$ being fitting parameters. Note that in the Heisenberg model in two dimensions, $\xi_{s}$ increases in the exponential form of \xi_{s}/a\sim\exp\big{[}b_{H}|J|/T\big{]} with $b_{H}\simeq 2\pi$ Heisenberg_Polyakov_75 . The resultant fitting function with the obtained values of $A_{H}=0.0017\pm 0.0003$ and $B_{H}=5.1\pm 0.1$ is represented by a dashed curve in Figs. 13 (b) and 7 (c). Since the obtained value of $B_{H}\simeq 5.1$ is comparable to $b_{H}\simeq 2\pi$, it turns out that $\sigma^{s}_{xx}\propto\xi_{s}$, which is in good agreement with the analytical result in Eq. (58). To get insight into the origin of the rapid increase of $\sigma^{s}_{xx}$, we examine the temperature dependences of $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ and $\tau_{s}$ like in the anisotropic cases of $\Delta\neq 1$.

The temperature dependences of $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ and $\tau_{s}$ are shown in Figs. 14, where $\tau_{s}$ is extracted from $\langle j^{z}_{s,x}(0)\,j^{z}_{s,x}(t)\rangle$ in the same way as before. $\langle|j^{z}_{s,x}(0)|^{2}\rangle$ is size dependent even at the lowest temperature, but its temperature dependence is relatively weak. In LSWT, as shown in Eq. (III.2), $\xi_{s}$ becomes relevant at lower temperatures, so that $L\gg\xi_{s}$ should be satisfied to evaluate the thermodynamic limit of $\langle|j^{z}_{s,x}(0)|^{2}\rangle$. As is suggested from the size dependent data, however, the maximum size of $L=384$ seems to be still small and the expected temperature dependence of $const\,T^{2}+T$ cannot be seen. Compared with $\langle|j^{z}_{s,x}(0)|^{2}\rangle$, the temperature dependence of $\tau_{s}$ is much more remarkable. As one can see from Fig. 14 (b), $\tau_{s}$ gets longer rapidly toward $T=0$, showing the considerably large system-size dependence. We fit the almost size-independent data at $0.54\leq T/|J|\leq 0.74$ with the exponential form of \tilde{A}_{H}\exp\big{[}\tilde{B}_{H}|J|/T\big{]}. The resultant curve with the obtained values of the fitting parameters $\tilde{A}_{H}=0.0006\pm 0.0002$ and $\tilde{B}_{H}=5.8\pm 0.2$ is represented by a dashed curve in Fig. 14 (b). The obtained value of $\tilde{B}_{H}\simeq 5.8$ is close to $b_{H}\simeq 2\pi$ and $B_{H}\simeq 5.1$. As $\sigma^{s}_{xx}$ is estimated roughly as $\sigma^{s}_{xx}\simeq T^{-1}\,\langle|j^{z}_{s,x}(0)|^{2}\rangle\,\tau_{s}$, the origin of the steep increase of $\sigma^{s}_{xx}$ toward $T=0$ is the enhanced relaxation-time $\tau_{s}$ which seems to have a direct association with the rapid growth of the spin correlation length \xi_{s}/a\sim\exp\big{[}b_{H}|J|/T\big{]}.

VI Summary and discussion

We have theoretically investigated transport properties of the classical antiferromagnetic XXZ model on the square lattice in which the anisotropy of the exchange interaction $\Delta\equiv J_{z}/J_{x}$ plays a role to control the universality class of the system. In Ising-type ($\Delta>1$), $XY$-type ($\Delta<1$), and Heisenberg-type ($\Delta=1$) magnets, spins in the low-temperature phase are, respectively, long-range-ordered via a magnetic phase transition, quasi-long-range-ordered via the KT topological transition, and disordered. Based on the linear response theory, we have calculated the thermal conductivity $\kappa_{\mu\nu}$ and the spin-current conductivity $\sigma^{s}_{\mu\nu}$ by means of the hybrid Monte-Carlo and spin-dynamics simulations. It is found that $\sigma^{s}_{\mu\nu}$ reflects the effect of the anisotropy, i.e., the difference in the ordering properties, while $\kappa_{\mu\nu}$ does not with its longitudinal component $\kappa_{xx}$ (=$\kappa_{yy}$) increasing toward $T=0$ as a power function of temperature independently of $\Delta$. For the $XY$-type anisotropy, the longitudinal spin-current conductivity $\sigma^{s}_{xx}$ ($=\sigma^{s}_{yy}$) exhibits a divergence at the Kosterlitz-Thouless (KT) transition temperature $T_{KT}$ obeying the exponential form, \sigma^{s}_{xx}\propto\exp\big{[}B/\sqrt{T/T_{KT}-1}\,\big{]} with $B={\cal O}(1)$, while for the Ising-type anisotropy, the temperature dependence of $\sigma^{s}_{xx}$ is almost monotonic without showing a clear anomaly at the magnetic transition temperature $T_{N}$. In the Heisenberg-type isotropic case, $\sigma^{s}_{xx}$ exhibits a monotonic exponential increase toward $T=0$. By analyzing the time correlation of the spin current at various temperatures, we find that the divergent enhancement of $\sigma^{s}_{xx}$ at $T_{KT}$ is due to the exponential rapid growth of the spin-current-relaxation time toward $T_{KT}$. Such a long spin-current-relaxation time can be interpreted as a manifestation of the topological nature of a vortex whose lifetime is expected to get longer toward $T_{KT}$ since the pair-annihilation of vortices should occur more sporadically with the increase of the inter-free-vortex distance toward $T_{KT}$. This suggests that the topological object of the vortex excitation should be crucial for the spin transport.

Now, we will address possible experimental platforms to investigate the pronounced enhancement of the longitudinal spin-current conductivity $\sigma^{s}_{\mu\mu}$ associated with the KT transition. As the divergent peak in the $\sigma^{s}_{\mu\mu}(T)$ curve toward $T_{KT}$ can commonly be seen in both ferromagnets and antiferromagnets only if an $XY$-type anisotropy exists, good candidate systems are quasi-two-dimensional magnets having the signature of the KT transition such as the $S=1/2$ square-lattice ferromagnet K2CuF4 KCuF_Hirakawa_jpsj_81 ; KCuF_Hirakawa_jpsj_82 ; KCuF_Hirakawa_JAP_82 ; KCuF_Sachs_prb_13 , the $S=1$ honeycomb-lattice antiferromagnets BaNi2X2O8 (X=As, P, V) BaNiXO_Regnault_JMMM_80 ; BaNiXO_Regnault_PhysicaB_86 ; BaNiPO_Regnault_PhysicaB_89 ; BaNiPO_Gaveau_JAP_91 ; BaNiVO_Heinrich_prl_03 ; BaNiVO_Waibel_prb_15 ; BaNiVO_Klyushina_prb_17 , the $S=5/2$ honeycomb-lattice antiferromagnet MnPS3 MnPS_Wildes_JPCM_98 ; MnPS_Ronnow_PhysicaB_00 ; MnPS_Wildes_prb_06 , and the stage-2 NiCl2 NiCl_Karimov_JETP_74 ; NiCl_Karimov_JETP_75 ; NiCl_Ikeda_JPC_81 and CoCl2 CoCl_Matsuura_JMMM_83 ; CoCl_Ikeda_jpsj_85 ; CoCl_Wiesler_Zphys_94 graphite intercalation which are respectively $S=1$ and $S=1/2$ triangular-lattice ferromagnets. In these compounds, a three-dimensional inter-layer coupling is extremely small, so that at first sight, the system may be regarded as a two-dimensional $XY$-type magnet. In reality, however, on approaching $T_{KT}$ at which the spin correlation length $\xi_{s}$ diverges, the effective coupling between neighboring layers grows rapidly as the area of the correlated region $\xi_{s}^{2}$ rapidly increases, eventually leading to a three-dimensional long-range-order as long as such a perturbative coupling is nonzero. Indeed, all the above compounds undergo a phase transition into a long-range-ordered state before reaching $T_{KT}$. Nevertheless, they have a two-dimensional $XY$-like crossover regime just above the magnetic transition, in which the critical phenomena peculiar to the KT transition have been observed. Thus, measurements of the spin-current conductivity in this crossover regime could, in principle, detect the pronounced enhancement of the longitudinal spin-current conductivity toward the virtually existing $T_{KT}$.

In the $XY$ magnets, the true divergence associated with the topological transition cannot be detected because the three-dimensional long-range-order inevitably appears before reaching $T_{KT}$. In Heisenberg magnets, however, such a divergence might be detectable if there exists a magnetic frustration leading to a non-collinear spin-ordering. In such frustrated Heisenberg magnets, a topological defect is the so-called $Z_{2}$ vortex and the KT-type $Z_{2}$-vortex transition is expected to occur at $T_{v}$ Z2_Kawamura_84 ; Z2_Kawamura_10 ; Z2_Kawamura_11 . In contrast to the KT transition, although the inter-free-vortex distance diverges at $T_{v}$, $\xi_{s}$ remains finite at any finite temperature. Thus, a divergent enhancement associated with the $Z_{2}$-vortex transition, if it occurs, is not necessarily masked by a three-dimensional long-range-order in real materials. This may be an interesting issue, but we will leave further detailed analysis for our future work.

As demonstrated in the present paper, the thermal transport is insensitive to the difference in the ordering properties. In extracting the magnetic contribution from the total longitudinal thermal conductivity, great care has to be taken because it contains phonon contribution as well in the temperature range typical for magnetic transitions. In contrast, the spin-current conductivity should be of purely magnetic origin unless a magnon-phonon coupling is strong enough, suggesting that the spin-current measurements may be a promising probe to detect nontrivial magnetic excitations such as vortices.

Acknowledgements.

The authors thank K. Uematsu, S. Furuya, and Y. Niimi for useful discussions. We are thankful to ISSP, the University of Tokyo for providing us with CPU time. This work is supported by JSPS KAKENHI Grant Numbers JP16K17748, JP17H06137.

Appendix A Ordering properties of the classical antiferromagnetic XXZ model on the square lattice

The ordering properties of the classical antiferromagnetic XXZ model (1) on the square lattice can be investigated by MC simulations. Figure 15 shows the temperature dependences of the specific heat $C$ (upper), the order parameter for the two-sublattice antiferromagnetic order (middle), the ratio of the spin-correlation length to the linear system size $L$ (bottom) for $\Delta=1.05$, 0.95, and 1. Here, in the Ising-type ($\Delta=1.05$), $XY$-type ($\Delta=0.95$), and Heisenberg-type ($\Delta=1$) spin systems, the order parameters and the associated spin-correlation lengths are respectively given by $m_{\rm AF}^{z}$ and $\xi_{s}^{z}$, $m_{\rm AF}^{xy}$ and $\xi_{s}^{xy}$, and $m_{\rm AF}$ and $\xi_{s}$ which are defined by

[TABLE]

In our MC simulations, we perform $3\times 10^{5}$ MC sweeps and the first $10^{5}$ sweeps are discarded for thermalization, where one MC sweep consists of the 1 heat-bath sweep and successive 10-30 over-relaxation sweeps. Observations are done in every MC sweep, and the statistical average is taken over 10 independent runs starting from different initial spin configurations.

As one can see from Fig. 15 (a), in the Ising case of $\Delta=1.05$, the specific heat $C$ exhibits a sharp peak associated with the antiferromagnetic transition at $T_{N}/|J|\simeq 0.75$. Correspondingly, $m_{\rm AF}^{z}$ starts growing up at $T_{N}$ and $\xi_{s}^{z}/L$ for different system sizes cross one another at $T_{N}$, which is usually the case for ordinary continuous magnetic phase transitions.

In the $XY$ case of $\Delta=0.95$, the KT transition temperature is estimated to be $T_{KT}/|J|\simeq 0.6$ in Refs. KT_XXZ_Cuccoli_95 ; KT_XXZ_Lee_05 ; KT_XXZ_Pires_96 . Actually, as one can see from Fig. 15 (b), $\xi_{s}^{xy}/L$ for different system sizes merge one another below $T_{KT}$, whereas the specific heat only shows a broad peak slightly above $T_{KT}$ and $m_{\rm AF}^{xy}$ is suppressed with increasing $L$ because of the absence of the true magnetic long-range order.

In the Heisenberg case of $\Delta=1$, the specific heat shows only a broad peak near $T/|J|\simeq 0.7$ and the spin-correlation length $\xi_{s}$ is finite at any finite temperature as is suggested from the fact that in Fig. 15 (c) $\xi_{s}/L$ continues to be suppressed with increasing the system size $L$ at all the temperatures.