\section{Introduction}
\label{sec:introduction}

In the research field of spintronics, various types of magnetoresistance such as giant magnetoresistance~\cite{Baibich1988,Binasch1989,Fert2008} and tunneling magnetoresistance~\cite{Julliere1975,Miyazaki1995,Moodera1995,Yuasa2004,Parkin2004} have been applied to spintronics devices for sensors, memories, and data storage.
Recently, a novel type of magnetoresistance, which called spin Hall magnetoresistance (SMR), has been attracting much attention.
The SMR has first been observed in a normal metal (NM)/ferromagnetic insulator (FI) bilayer\cite{Nakayama2013a,Chen2013a,Hahn2013,Vlietstra2013,Althammer2013,Meyer2014,Marmion2014,Cho2015,Kim2016,Chen2016,Sterk2019,Tolle2019} and subsequently in a NM/antiferromagnetic insulator (AFI) bilayer\cite{Hou2017,Lin2017,Cheng2019,Hoogeboom2017a,Fischer2018,Ji2018,Lebrun2019}. 
Because the SMR reflects information on a FI magnetization (or AFI N\'{e}el vector), it can be utilized for detection of the orientation of the ordered spins in magnetic materials.

The SMR also has a potential possibility as a useful probe for two-dimensional magnetic materials such as van der Waals atomic layers.
It is remarkable that several two-dimensional atomic layer compounds showing ferromagnetism\cite{2d_fi} and antiferromagnetism\cite{2d_afi} have been synthesized recently.
Since the SMR measurement for such a magnetic atomic layer will be possible in a near future, it is a urgent task to construct the theory of the SMR which is applicable to two-dimensional quantum magnets.

\begin{figure}[tb]
 \begin{center}
  \includegraphics[clip,width=7.0cm]{fig/smr_fig.pdf}
    \caption{Setup of spin Hall Magnetoresistance. Spin absorption at the interface changes according to the orientation of the Neel vector of the antiferromagnetic insulator, and the magnitude of the magnetoresistance changes.}
    \label{fig:SMR_setup}
  \end{center}
\end{figure}

The existing theory\cite{Nakayama2013a,Chen2013a} explains the SMR as follows (see Fig.~\ref{fig:SMR_setup}).
When an in-plane charge current flows in the NM layer, there occurs spin accumulation near the NM/FI(AFI) interface due to the spin Hall effect\cite{Sinova2015}.
Then, a backflow spin current, which is induced by spin diffusive in the NM, is converted into the charge current again by the inverse spin Hall effect (ISHE), resulting in the longitudinal magnetoresistance.
This longitudinal magnetoresistance depends on the orientation of FI magnetization (or AFI N\'{e}el vector),
because the amount of spin accumulation at the interface is changed by it.
We note that the strength of SMR is proportional to square of the spin Hall angle of the NM.

Although this semiclassical theory\cite{Nakayama2013a,Chen2013a} succeeded in explaining qualitative features of the SMR measurement, it could not explain temperature dependence of the SMR signal.
Recently, two of the authors have constructed a microscopic theory using Green's function~\cite{Kato2020}. 
In this microscopic theory, the SMR is described in terms of local spin susceptibilities of the FI(AFI).
In particular, this microscopic theory can include dynamic processes such as magnon absorption and emission, which are neglected in the semiclassical theory.

In Ref.~\onlinecite{Kato2020}, the SMR has been calculated analytically by employing the spin-wave approximation.
However, the spin-wave approximation cannot be applied near the transition temperature.
Furthermore, it becomes inaccurate when the magnitude of the localized spin in the FI(AFI) is small.
Therefore, to obtain the features of the SMR in a wide range of the temperature and arbitrary magnitude of the localized spin, $S$, we need to calculate it without employing the spin-wave approximation.

In our study, we numerically calculate the SMR using quantum Monte Carlo method based on the formulation of Ref.~\onlinecite{Kato2020}.
We utilize a new method for accurately obtaining an integral of the local spin susceptibilities from the QMC data without numerical analytic continuation.
We evaluate detailed temperature dependence of the SMR for the $S=1/2$ and $S=1$ spin systems on a two-dimensional square lattice and a quasi-two-dimensional cubic lattice with a finite number of layers.
Based on the numerical results for these models, we discuss qualitative features of the SMR.

This work is organized as follows. In Sec.\ref{sec:formulation}, we formulate SMR by the microscopic theory based on a microscopic Hamiltonian of NM, AFI, and exchange interaction at interface. In Sec.\ref{sec:Numerical}, we propose numerical calculation method for SMR with high accuracy. 
In Sec.\ref{sec:Result}, we show numerical results of the SMR obtained from quantum Monte Carlo simulation.
Finally, we discuss experimental relevance of our work in Sec.~\ref{sec:discussion} and summarize our results in Sec.~\ref{sec:summary}.
We provide the details of the derivation and comparison with theoretical calculation in Appendixes.

\section{Formulation}
\label{sec:formulation}

\subsection{Spin conductance}
\label{sec:SpinConductance}

We first summarize a theoretical framework for SMR~\cite{Chen2013a,Kato2020}, which is employed in this work.
When an electric field is applied to the metal side in the $x$-direction, a charge current $j_x$ induces a spin current $j_s^{\rm SH}=\theta_{\rm SH} j_x$ in the $y$-direction due to the spin Hall effect, where $\theta_{\rm SH}$ is the spin Hall angle. 
As a result, the electron spins parallel to the $z$ axis are accumulated near the interface.
This spin accumulation is described by the spin chemical potential $\mu_{s}(y)=\mu_{\uparrow}(y)-\mu_{\downarrow}(y)$, where $\mu_{\uparrow}(y)$ and $\mu_{\downarrow}(y)$ are the chemical potentials of two different spins. 
The spatial gradient of $\mu_{s}(y)$ generates a backflow $j_s^{\rm B}$ in the $y$ direction due to spin diffusion.
Finally, this backflow is converted into a charge current in the $x$ direction due to the inverse spin Hall effect, which results in the SMR. 

The amount of the spin accumulation near the interface is affected by the spin loss rate at the interface $I_{S}$, which is assumed to be
proportional to the spin chemical potential $\mu_{s}(0)$, where the position of the interface is set as $y=0$.
Following Ref.~\onlinecite{Kato2020}, we introduce the spin conductance $G_s$ as
\begin{align}
G_s &= \lim_{\mu_{s}(0)\rightarrow 0} \frac{I_{S}}{\mu_{s}(0)}.
\label{eq:spinconductancedef}
\end{align}
We note that $G_s$ depends on the orientation of the N\'{e}el vector of the AFI.
Combining the spin diffusion theory and the spin conductance, the SMR ratio can be derived as
\begin{align}
\frac{\Delta \rho}{\rho} &= \theta_{\rm SH} \gamma \tanh (d/2\lambda), \\
\gamma &= \frac{4e^{2}}{\hbar} \frac{G_s}{S\sigma / \lambda},
\end{align}
where $e$ is elementary charge $e>0$, $S$ is surface area of interface, $\lambda$ is spin diffusion length, $d$ is a thickness of the NM, and $\gamma$ is a normalized spin conductance (we have assumed that $\gamma \ll 1$).
Thus, calculation of the SMR signal is attributed to that of the spin conductance $G_s$.

In the remaining part of Sec.~\ref{sec:formulation}, we formulate the spin current $I_s$ and the spin conductance $G_s$ as a function of the N\'{e}el vector of the AFI. 

\subsection{Normal metal}
\label{sec:NM}

We describe a normal metal by noninteracting electron system whose Hamiltonian is given as
\begin{align}
H_{\rm{NM}}=\sum_{\bm{k}\sigma} \epsilon_{\bm{k}} c^{\dag}_{\bm{k}\sigma}c_{\bm{k}\sigma} .
\label{eq:Ham}
\end{align}
where $\epsilon_{\bm{k}}$ is a kinetic energy and $c_{{\bm k}\sigma}$ is an annihilation operator of conduction electrons with a wavenumber ${\bm k}$ and a spin ${\bm \sigma}$.
The spin accumulation at the interface is modeled by quasi-equilibrium distribution that is described by the effective Hamiltonian
\begin{align}
\mathcal{H}_{\rm{NM}}&= H_{\rm NM} - \sum_{\sigma} \mu_\sigma N_{\sigma} ,
\end{align}
where $\mu_\sigma$ is a spin-dependent chemical potential and $N_\sigma$ is a number operator of conduction electrons with a spin $\sigma \equiv  \pm 1$.
This effective Hamiltonian is rewritten as
\begin{align}
\mathcal{H}_{\rm{NM}}=
\sum_{\bm{k}}(\xi_{\bm{k}} - \sigma
\delta \mu_s /2)c^{\dag}_{\bm{k}\sigma}c_{\bm{k}\sigma}.
\label{eq:calHam}
\end{align}
where $\xi_{\bm{k}}=\epsilon_{\bm{k}}-\mu$ is a kinetic energy measured from the averaged chemical potential $\mu=(\mu_\uparrow + \mu_\downarrow)/2$ and $\delta \mu = \mu_\uparrow - \mu_\downarrow$ is a spin chemical potential near the interface.

As shown in Sec.~\ref{sec:SpinCond2}, the spin conductance is written in terms of a local spin susceptibility defined as
\begin{align}
\chi^{R}_{\rm{loc}}(\omega) &={\frac{1}{{N}_{\rm{NM}}}}\sum_{\bm{q}}
\int dt \, \chi^{R}_{+-}(\bm{q},t) e^{i\omega t},
\label{eq:3_spin_sus_loc} \\
    \chi^{R}_{+-}(\bm{q},t) &= {\frac{i}{{{N}_{\rm{NM}}} \hbar}}\theta(t)
    \langle [s^{+}_{\bm{q}}(t),s^{-}_{\bm{q}}(0)] \rangle , 
\end{align}
where $\langle \cdots \rangle = {\rm Tr}\, (e^{-\beta {\cal H}_{\rm NM}}\cdots)/{\rm Tr}\, e^{-\beta {\cal H}_{\rm NM}}$ indicates a average with respect to the quasi-equilibrium state, ${N}_{\rm{NM}}$ is the number of unit cells in normal metal, $s^{\pm}_{\bm{q}}$ are spin ladder operators defined as
\begin{align}
    s^{+}_{\bm{q}}&=(s^{-}_{\bm{q}})^{\dagger} = \sum_{\bm{k}} c^{\dag}_{\bm{k}+{\bm q}\uparrow}c_{\bm{k}\downarrow} ,
\end{align}
and $s^+_{\bm q}(t) = e^{iH_{\rm NM}t/\hbar} s^+_{\bm q}e^{-iH_{\rm NM}t/\hbar}$.
For the present model described by Eqs.~(\ref{eq:Ham}) and (\ref{eq:calHam}), the imaginary part of the local spin susceptibility is calculated as~\cite{Kato2020}
\begin{align}
    {\rm Im} \, \chi^{R}_{\rm{loc}}(\omega) 
    = \pi N(0)^2 (\hbar \omega + \delta \mu_s),
    \label{eq:chiloc}
\end{align}
where $N(0)$ is a density of states near the Fermi energy.

\subsection{Antiferromagnetic insulator}
\label{sec:AFI}

\begin{figure}[tb]
 \begin{center}
  \includegraphics[clip,width=5.0cm]{fig/axis_change.pdf}
    \caption{Coordinate transformation from the laboratory system $O$-$xyz$ to the coordinate system $O$-$x'y'z'$ with the $z'$ axis fixed to the magnetization orientation.}
    \label{fig:smr_axis}
  \end{center}
\end{figure}

Let us denote an azimuth angle of the N\'{e}el vector of the AFI measured from the $z$ axis with $\theta$.
We introduce a new magnetization-fixed coordinate $O$-$x'y'z'$, which is obtained by rotating the laboratory coordinate $O$-$xyz$ by $\theta$ around the $y$-axis as shown in Fig.~\ref{fig:smr_axis}.
When the components of the spin operator ${\bm S}$ in these two coordinates are denoted with $(S^x, S^y, S^z) $ and $(S^{x'}, S^{y'}, S^{z'})$, the transformation between them is given as
\begin{align}
\left( \begin{array}{c} S^{x'} \\ S^{y'} \\ S^{z'} \end{array}\right)
= \left(\begin{array}{ccc}
\cos \theta & 0 & - \sin \theta \\
0 & 1 & 0 \\
\sin \theta \, & 0 & \cos \theta \, \end{array} \right) 
\left( \begin{array}{c} S^{x} \\ S^{y} \\ S^{z} \end{array}\right)
.
\label{spin_coordinate_transformation}
\end{align}

For the AFI, we consider the quantum Heisenberg model in which localized spins are located on a cubic lattice.
We divide the lattice into A and B sublattices and define the spin operator at a site $l$ on the sublattice $\nu$ ($=$A,B) as ${\bm S}_{\nu,l}$.
Then, the Hamiltonian of the AFI is described using the spin operators in the magnetization-fixed coordinate as
\begin{align}
H_{\rm AFI} = J \sum_{\langle {l,l'} \rangle}
\Biggl[
{\frac{1}{2}} \left( S^{+'}_{{\rm A},l}S^{-'}_{{\rm B},l'} +S^{-'}_{{\rm A},l}S^{+'}_{{\rm B},l'} \right) +S^{z'}_{{\rm A},l}S^{z'}_{{\rm B},l'} \Biggl] \nonumber \\
+ D \sum_{l} \Biggl[ (S^{z'}_{{\rm A},l})^2 + (S^{z'}_{{\rm B},l})^2 \Biggl]+ h \sum_{l} \Biggl[ S^{z'}_{{\rm A},l} - S^{z'}_{{\rm B},l} \Biggl],
\label{eq:HAFI}
\end{align}
where $\langle l,l'\rangle$ indicates a pair of the nearest-neighbor sites, $S_{\nu,l}^{\pm'}=S_{\nu,l}^{x'}\pm i S_{\nu,l}^{y'}$ are spin ladder operators, and $J$ is a magnitude of the exchange interaction.
In order to fix the N\'{e}el vector in the quantum Monte Carlo simulation, we consider a small anisotropy of magnetization and an alternating magnetic field, whose amplitudes are denoted with $D$ and $h$, respectively.
In the following discussion, we set $J$ as a unit of the energy.

The spin conductance is written in terms of a local spin correlation functions defined as
\begin{align}
    G^{R,(a)}_{{\rm loc},\nu\nu'}(\omega) &= \frac{1}{N_{\rm AFI}} \sum_{\bm q} \int dt \, G^{R,(a)}_{\nu\nu'}(\bm{q},t) e^{i\omega t},
\end{align}
where $N_{\rm AFI}$ is a number of unit cells of the AFI and $G^{R,(a)}_{\nu\nu'}(\bm{q},t)$ ($a=1,2,3$) indicate three kinds of the spin correlation functions defined as
\begin{align}
G^{R,(1)}_{\nu\nu'}(\bm{q},t) &= -\frac{i}{\hbar} \theta(t)\langle{[S^{z'}_{\nu,\bm{q}}(t),S^{z'}_{\nu',{\bm{q}}}(0)]} \rangle,\\
G^{R,(2)}_{\nu\nu'}(\bm{q},t) &= -\frac{i}{\hbar} \theta(t)\langle{[S^{+'}_{\nu,\bm{q}}(t),S^{-'}_{\nu',\bm{q}}(0)]} \rangle, \\
G^{R,(3)}_{\nu\nu'}(\bm{q},t) &= -\frac{i}{\hbar}
    \theta(t)\langle{[S^{-'}_{\nu,\bm{q}}(t),S^{+'}_{\nu',\bm{q}}(0)]} \rangle.
\end{align}
Here, $\langle \cdots \rangle$ indicates a thermal average in the AFI and 
the Fourier transformations of the spin operators are defined as
\begin{align}
S_{\nu,{\bm q}}^{+'} &= (S_{\nu,{\bm q}}^{-'})^\dagger =\sum_{l} S_{\nu,l}^{+'} e^{-i{\bm q}\cdot {\bm R}_{\nu,l}}, \\
S_{\nu,{\bm q}}^{z'} &= \sum_{l} S_{\nu,l}^{z'} e^{-i{\bm q}\cdot {\bm R}_{\nu,l}}, 
\end{align}
where ${\bm R}_{\nu,l}$ indicates a position of the site $l$ on the sublattice $\nu$.

\subsection{Microscopic description of spin conductance}
\label{sec:SpinCond2}

We consider an exchange coupling at the interface between the NM and AFI, whose Hamiltonian is given as
\begin{align}
    H_{\rm{ex}}=\sum_{\bm{k},\bm{q},\nu}[{\cal T}^\nu_{\bm{k},\bm{q}}{\mathcal{S}}^{+}_{\nu\bm{k}}s^{-}_{\bm{q}}+({\cal T}^\nu_{\bm{k},\bm{q}})^* {\mathcal{S}}^{-}_{\nu\bm{k}}s^{+}_{\bm{q}}],
 \label{eq:Hex}
\end{align}
using the laboratory coordinate, where ${\cal T}^\nu_{\bm{k},\bm{q}}$ is magnitude of exchange interaction at the interface.
In this work, we set ${\cal T}^\nu_{\bm{k},\bm{q}} = {\cal T}$ assuming that ${\cal T}^\nu_{\bm{k},\bm{q}}$ is independent of ${\bm k}$, ${\bm q}$, and $\nu$.
By using the transformation (\ref{spin_coordinate_transformation}),
we can rewrite Hamiltonian in the magnetization-fixed coordinate as
\begin{align}
H_{\rm ex} &= \sum_{a=1}^3 H_{\rm{ex}}^{(a)}, 
\label{eq:Hex1} \\
    H_{\rm{ex}}^{(a)} &=g_{a}(\theta) \sum_{\bm{k},\bm{q},\nu}[{\cal T} {S}^{(a)}_{\nu\bm{k}}s^{-}_{\bm{q}}+{\cal T}^* ({S}^{(a)}_{\nu\bm{k}})^\dagger s^{+}_{\bm{q}}], \label{eq:Hex2}
\end{align}
where $S_{\nu{\bm k}}^{(a)}$ and $g_a(\theta)$ ($a=1,2,3$) are defined as
\begin{align}
    & {S}^{(1)}_{\nu\bm{k}}={S}^{z'}_{\nu\bm{k}},
    \quad g_{1}(\theta)=-\sin{\theta}, \label{eq:Hex3} \\
    & {S}^{(2)}_{\nu\bm{k}}={S}^{+'}_{\nu\bm{k}},
    \quad g_{2}(\theta)=\cos^2(\theta/2), \label{eq:Hex4} \\
    & {S}^{(3)}_{\nu\bm{k}}={S}^{-'}_{\nu\bm{k}},
    \quad g_{3}(\theta)=-\sin^2 (\theta/2). \label{eq:Hex5}
\end{align}
A detailed derivation is given in Appendix~\ref{app:SpinOperator}.

By the second-order perturbation with respect to $H_{\rm ex}$, the spin conductance is calculated as~\cite{Kato2020}
\begin{align}
G_s &= G_0 g_1(\theta)^2 \langle S_{\nu,l}^{z'} \rangle^2 \nonumber \\
&+\sum_{a=1}^{3} \sum_{\nu}
2 G_0 g_{a}(\theta)^2 \int \frac{d\epsilon}{2 \pi}\, \left. {\rm Im} \, \chi_{\rm loc}^R(\epsilon/\hbar) \right|_{\delta \mu_s=0}
    \nonumber \\ 
& \hspace{10mm} \times (-{\rm{Im}}\, G^{R,(a)}_{\nu\nu,\rm{loc}}(\epsilon/\hbar)) \left( - \frac{\partial f}{\partial \epsilon} \right),
\label{eq:Gs}
\end{align}
where $f(\epsilon)=(e^{\beta\epsilon}-1)^{-1}$ is the Bose distribution function and $G_0=2\pi |{\cal T}|^2 N_{\rm NM}^2 N_{\rm AFI} N(0)^2$ is a dimensionless parameter which represents a strength of the interfacial exchange coupling.
A detailed derivation is given in Appendix~\ref{app:SpinCurrent}.

The amplitude of the SMR signal is proportional to\footnote{In Ref.~\onlinecite{Hou2017}, the orientation of the N\'{e}el vector in the AFI is controlled by magnetization of an adjacent YIG.
The negative SMR observed in this experiment can be reasonably understood if the N\'{e}el vector is perpendicular to the magnetization of YIG.
In our work, we also employ this assumption.}
\begin{align}
    \Delta G &\equiv G(\theta = 0)-G(\theta = \pi/2).
\end{align}
For convenience of discussion, we express $\Delta G$ by a sum of the two contributions as 
\begin{align}
& \Delta G = \Delta G_{z} + \Delta G_{xy}, 
\label{eq:Gsum} \\
& \Delta G_z/G_0 = - \langle S^{z'}_{\nu,l} \rangle^2 - 2J_1 , \\
& \Delta G_{xy}/G_0 = \frac{3}{2}J_2 -\frac{1}{2} J_3 , \\
& J_a = \sum_\nu \int_{-\infty}^\infty \frac{d\epsilon}{2\pi}
    \left[
    (-{\rm{Im}}\, G^{R,(a)}_{\rm{loc},\nu\nu}(\epsilon/\hbar))\frac{\beta \epsilon}{\sinh^2(\beta \epsilon/2)}
    \right], \label{eq:Jadef}
\end{align}
where $\Delta G_z$ and $\Delta G_{xy}$ are contributions from the spin-spin correlation in the $z$- and $xy$-directions, respectively.
We will show later that $\Delta G_z$ ($\Delta G_{xy}$) gives a negative (positive) contribution to the SMR signal. 
We note that the integrals, $J_a$ ($a=1,2,3$), vanish at zero temperature.
Therefore, the negative SMR signal $\Delta G/G_0 = -\langle S^{z'}_{\nu,l} \rangle^2$ is obtained at zero temperature in consistent with the experiment~\cite{Hou2017}.
Thus, the SMR signal is formulated in terms of the spin correlation function $G_{\nu \nu,{\rm loc}}^{R,(a)}(\omega)$ and the staggered magnetization $|\langle S^{z'}_{\nu,l} \rangle|$.
In the next section, we explain a numerical method to calculate them.

\section{Numerical method}
\label{sec:Numerical}

We employ the continuous-time quantum Monte Carlo (QMC) method to calculate the spin correlation functions and the staggered magnetization for a finite temperature~\cite{Evertz2003,Kawashima2004,Gubernatis2016}.
We use the program package, Discrete Space Quantum Systems Solver (DSQSS)~\cite{Motoyama2021,DSQSSWeb}, which implements a continuous-time path-integral QMC method based on a directed-loop algorithm~\cite{Syljuasen2002,Syljuasen2003} and the Suwa-Todo algorithm without detailed balance condition~\cite{Suwa2010}.
To calculate the spin correlation functions $G_{\nu \nu,{\rm loc}}^{R,(a)}(\omega)$, we need to perform an analytic continuation from the Matsubara frequency $i\omega_n$ to the real frequency $\omega$.
However, the numerical analytic continuation is usually unstable and inaccurate for obtaining reliable results.

\begin{figure}[tb]
 \begin{center}
    \includegraphics[clip,width=8.0cm]{fig/path_change.pdf}
    \caption{Change of integral path.
    (a) is before change integral path and needs to calculate all $G^{R}_{\rm{loc}}(\omega)$.  (b) is after change path and can calculate only Matsubara frequency $G^{R}_{\rm{loc}}(i\omega_{n})$ and $G^{\prime R,(a)}_{\rm{loc}}(i\omega_{n})$.}
    \label{fig:integral_change}
  \end{center}
\end{figure}

In this work, we employ a new method without performing numerical analytic continuation and directly evaluate the integrals, $J_a$, defined by Eq.~(\ref{eq:Jadef}).
These integrals shares a function $\beta \epsilon /\sinh^2(\beta \epsilon/2)$, which have double poles at $\epsilon = 2 i\pi n/\beta \equiv i\epsilon_n$ where $n$ is an integer
\footnote{We note that $G^{R,(a)}_{\rm{loc}}(\epsilon/\hbar)$ vanishes at $\epsilon=0$ and therefore the integrand has no pole there.}.
Therefore, by modifying the integration path in the complex plane of $\epsilon$ as in Fig.~\ref{fig:integral_change}, $J_a$ can be expressed by a sum of the residues of the integrand at $\epsilon = i\epsilon_n$ ($n\ge 1$).
Since the residues are evaluated as
\begin{align}
    & R^{(a)}_{\nu,n} \equiv \underset{\epsilon = i\epsilon_{n}}{{\rm{Res}}}\frac{\beta \epsilon G^{R,(a)}_{\rm{loc},\nu\nu}(\epsilon/\hbar) e^{\beta  \epsilon}}{(e^{ \beta\epsilon}-1)^2}\nonumber \\
    & = \frac{1}{\beta} \left[ G^{R,(a)}_{\rm{loc},\nu\nu}(i\epsilon_{n}/\hbar) +  \epsilon_{n} \frac{dG_{\rm{loc},\nu\nu}^{R,(a)}(i\epsilon_n/\hbar)}{d\epsilon_n} \right] ,
\end{align}
we can express the integrals as

\begin{align}
    J_a = - \sum_{n=1}^\infty \sum_{\nu} {\rm Re}\, R^{(a)}_{\nu,n} .
    \label{eq:sumformula}
\end{align}
Thus, the integrals can be evaluated only from the information of the imaginary-times spin correlation functions.
The derivative of $G^{R,(a)}_{\rm{loc},\nu\nu}(i\epsilon_n/\hbar)$ with respect to the Matsubara frequency can be obtained by numerically differentiation using the Pade approximation.
Here, the Pade approximation is simply used for interpolation on the imaginary axis, and therefore the accuracy of the simulation is greatly improved compared with direct numerical analytic continuation.

In our Monte Carlo simulation, we typically used $10^6$ Monte Carlo samples for each point.
To improve the accuracy, the mesh number of the histogram along the imaginary time was increased at low temperatures.
In performing the sum in Eq.~(\ref{eq:sumformula}), we used numerical data of the spin correlation function for the Matsubara frequencies below a cutoff frequency and extrapolated it in the form $C/\epsilon_n^2$ ($C$: a constant) for higher Matsubara frequencies. 

\section{Result}
\label{sec:Result}

In this section, we show the numerical results of the amplitude of the SMR signal, $\Delta G/G_0$, for the $S=1/2$ and $S=1$ spin systems on a two-dimensional square lattice and a quasi-two-dimensional cubic lattice with a finite number of layers.
We also show the two contributions, $\Delta G_z/G_0$ and $\Delta G_{xy}/G_0$, separately (see Eqs.~\eqref{eq:Gsum}-\eqref{eq:Jadef}).
Finally, we discuss the effect of randomness using a model with disordered exchange interactions.

We note that the quasi-two-dimensional quantum Heisenberg model has no long-range order at finite temperatures because of the Mermin–Wagner theorem.
However, the staggered magnetization grows rapidly below a specific temperature, at which the correlation length exceeds the system size at low temperatures.
Since this rapid growth of the staggered magnetization is expected to simulate the actual behavior of thin antiferromagnet layers, our numerical results are used for discussing qualitative features of the SMR signal (for a detailed discussion, see Sec.~\ref{sec:discussion}).

\subsection{Case of $S=1/2$}
\label{sec:s1/2}

\begin{figure}[tb]
 \begin{center}
    \includegraphics[clip,width=8.0cm]{fig_result/m1.pdf}
    \caption{Temperature dependence of the spin conductance $\Delta G/G_0 (= \Delta G_{xy}/G_0+\Delta G_z/G_0)$ for the $S=1/2$ quantum Heisenberg model with the
    alternating magnetic field $h=0.01J$ on (a) a $L\times L$ square lattice ($L=24,32,48,64$) and (b) a $16\times 16 \times W$ cubic lattice ($W=2,4,6$). The two contributions, $\Delta G_{xy}/G_0$ and $\Delta G_{z}/G_0$, are also plotted separately.
    }
    \label{fig:smr_m1_2d_l}
  \end{center}
\end{figure}

We first show the SMR signal for the $S = 1/2$ quantum Heisenberg model on a $L\times L$ square lattice ($L=24, 32, 48, 64$) in Fig.~\ref{fig:smr_m1_2d_l}~(a).
The black plots indicate the total SMR signal, $\Delta G/G_0$, whereas the red and blue plots indicate $\Delta G_z/G_0$ and $\Delta G_{xy}/G_0$, respectively.
To fix the N\'{e}el vector to the $z$ direction, we introduced a weak alternating magnetic field, $h=0.01J$ (see also Eq.~\eqref{eq:HAFI}).
We find that the size dependence of the SMR signal is weak for $L\ge 24$.
The solid curves in the inset of Fig.~\ref{fig:smr_m1_2d_l} shows the staggered magnetization as a function of the temperature for $L=16,32,64$. 
We can see that its size dependence becomes weak for $L\ge 32$.
For rough estimate of the temperature at which the staggered magnetization starts to grow, we perform fitting in the form of $|\langle S_i \rangle| \propto 1-(T/T_c)^\alpha$ in the range of
$|\langle S_i \rangle /S_0| > 0.2$ for $L=64$ as indicated by the red dashed line in the inset of Fig.~\ref{fig:smr_m1_2d_l}).
For convenience, this characteristic temperature is called the ordering temperature hereafter.
The estimated ordering temperature is indicated by the green arrow in the main graph of Fig.~\ref{fig:smr_m1_2d_l}.
Above the transition temperature, the SMR signal vanishes because $\Delta G_z/G_0$ and $\Delta G_{xy}/G_0$ cancels each other out.
This feature can be understood analytically by the high-temperature expansion (see Appendix~\ref{app:high-temp}).
In contrast, we find that the SMR signal becomes finite below the ordering temperature.
As the temperature is further lowered below the transition temperature, the SMR signal once increases, and then decreases toward zero temperature; the spin conductance has a peak roughly at three-fifths of the transition temperature.
Although the SMR signal is predicted to be negative at low temperatures as predicted by the spin-wave approximation~\cite{Kato2020}, such sign change is not found in Fig.~\ref{fig:smr_m1_2d_l}~(a) down to the lowest simulation temperature ($k_{\rm B}T/J=0.125$).

Fig.~\ref{fig:smr_m1_2d_l}~(b) shows the result for a quasi-two-dimensional system with a finite thickness, i.e., a $16 \times 16 \times W$ cubic lattice ($W=2,4,6$).
As the thickness $W$ increases, the ordering temperature (indicated by the green arrows) is enhanced and the spin conductance induced below the transition temperature becomes large.
We find that the peak of the spin conductance is largely enhanced compared with that for a single-layer spin system ($W=1$).
Although the maximum of the spin conductance increases with increasing the thickness $W$, the difference between the results for $W=4$ and $6$ is much smaller than that for $W=2$ and $4$.

\subsection{Case of $S=1$}
\label{sec:s1}

Next, we consider the $S = 1$ quasi-two-dimensional quantum Heisenberg model with a finite thickness.
We introduce a small axial anisotropy $D=-0.1J$ for fixing the N\'{e}el vector.
The SMR signal for a $L\times L\times 6$ cubic lattice ($L=8, 12, 16$) is shown in Fig.~\ref{fig:smr_m2_3d_thi}~(a).
We find that the size dependence between $L=12$ and $L=16$ is weak enough for examining qualitative features.
The SMR signal for a $16\times 16\times W$ cubic lattice ($W=2, 4, 6$) is shown in Fig.~\ref{fig:smr_m2_3d_thi}~(b).
The ordering temperature (indicated by the green arrows) is enhanced for $S=1$ compared with the case of $S=1/2$ (see Fig.~\ref{fig:smr_m1_2d_l}~(b)) because quantum fluctuations are suppressed.
As the thickness $W$ increases, the transition temperature is raised and the SMR signal becomes large.
As in the case of $S=1/2$, the difference between the results for $W=4$ and $6$ is much smaller than that for $W=2$ and $4$.

Above the transition temperature, the SMR signal becomes small but remains a finite value in contrast to the case of $S=1/2$.
This feature can be understood by the high-temperature expansion (see Appendix~\ref{app:high-temp}).
As the temperature is lowered below the ordering temperature, the SMR signal once increases, takes a maximum roughly at three-fifths of the ordering temperature, and then decreases.
The maximum value of the SMR signal is rather enhanced compared with the case of $S=1/2$.
At sufficiently low temperatures, the SMR signal becomes negative as predicted by the spin-wave approximation (for a detailed comparison, see Appendix~\ref{app:sw}).

\begin{figure}[tb]
 \begin{center}
    \includegraphics[clip,width=8.0cm]{fig_result/m2.pdf}
    \caption{Temperature dependence of the spin conductance $\Delta G/G_0 (= \Delta G_{xy}/G_0+\Delta G_z/G_0)$ for the $S=1$ quantum Heisenberg model with an axial anisotropy $D=-0.1J$ on (a) a $L\times L\times 6$ cubic lattice ($L=8,12,16$) and (b) a $16\times 16 \times W$ cubic lattice ($W=2,4,6$). The two contributions, $\Delta G_{xy}/G_0$ and $\Delta G_{z}/G_0$, are also plotted separately.
    }
    \label{fig:smr_m2_3d_thi}
  \end{center}
\end{figure}

\subsection{Disordered case for $S=1$}
\label{sec:disordered}

Finally, let us consider effect of randomness by introducing disordered exchange interactions.
We consider disordered exchange interactions whose probability density function~\cite{Laflorencie_random_2006} is given by
\begin{align}
P(J) =  J^{-1+\delta^{-1}} \delta^{-1} \Theta (J) \Theta (1-J),
\label{eq6:dis}
\end{align}
where $\Theta(x)$ is the Heaviside step function and $\delta$ is a parameter
which is chosen as 1.5 in our simulation.

\begin{figure}[tb]
 \begin{center}
    \includegraphics[clip,width=8.0cm]{fig_result/m2_dis.pdf}
\caption{Temperature dependence of $\Delta G/G_0$ for the disordered $S=1$ quantum Heisenberg model on a $16\times 16 \times 6$ cubic lattice.
The parameter $\delta$ in the probability distribution (\ref{eq6:dis}) is chosen as $1.5$.}
    \label{fig:smr_m2_3d_dis}
  \end{center}
\end{figure}

We show the SMR signal for the disordered $S=1$ quantum Heisenberg model on a $16\times 16 \times 6$ cubic lattice.
The two legends indicate the results for the uniform and disordered cases, respectively.
We note that the value of the exchange interaction $J$ is taken as constant for the uniform system.
In the disorder case, the ordering temperature is suppressed and the SMR signal becomes small.
However, the qualitative features are common as the uniform case; the SMR signal becomes small above the transition temperature whereas it once increases and decreases as the temperature is lowered below the ordering temperature.

\section{Experimental Relevance}
\label{sec:discussion}

In this section, we discuss how our simulation is related to the SMR measurements.
As already stated, our calculation has several limitations in comparing it with experimental results.
Because of the Mermin-Wagner theorem~\cite{Mermin66}, quasi-two-dimensional quantum spin systems never have a finite-temperature phase transition whereas its correlation length diverges toward zero temperature.
To capture this feature of quasi-two-dimensional quantum systems, we need to perform large-scale QMC simulation for large system sizes and to carefully analyse numerical data by the finite-size scaling.
Furthermore, we introduced additional parameters, i.e., a staggered magnetic field and an anisotropy in our model to fix the direction of the AFI N\'{e}el vector.
These additional terms in the Hamiltonian, as well as finite-size effect, smear out rapid rise of the staggered magnetization near the ordering temperature.
On the other hand, the AFI/NM bilayer used in experiments has, more or less, a finite size.
By experimental methods such as a neutron scattering, we usually observe that the staggered magnetization of the AFI layer grows rapidly below a specific finite temperature, which is frequency called the transition temperature.
In this sense, our simulation, in which the antiferromagnetic ordering grows below a specific temperature due to the finite-size effect, resembles the behavior in real experiments.
Therefore, qualitative features of our simulation are expected to be observed experimentally.

In our simulation, characteristic features are summarized as follows.
(a) When the magnitude of the localized spin, $S$, decreases, the SMR signal becomes small.
In particular, the SMR signal is rather suppressed for the case of $S=1/2$.
(b) The SMR signal becomes large below the transition temperature.
The peak of the SMR signal {\it does not} correspond to the transition temperature.
(c) The SMR signal becomes negative at sufficiently low temperatures for the case of $S=1$.
This feature is consistent with the previous theory based on the spin-wave approximation~\cite{Kato2020}.
Although we expect that the sign change of SMR occurs at low temperatures also for the case of $S=1/2$, we need to perform QMC simulation at lower temperatures to check it.
(d) Although randomness of exchange interactions reduces both the transition ordering temperature and the SMR signal, it does not change qualitative features in the temperature dependence of the SMR.

The features obtained in our simulation are consistent with the experiment for the NiO thin layer~\cite{Hou2017}, in which a $S=1$ spin exists per Ni site.
However, there are also differences; The SMR never becomes zero even at high temperatures and the curves of the SMR signal in our simulation does not completely fit to the experimental result.
We attribute these differences to
additional factors such as domain formation, roughness of interfaces, proximity effect by adjacent magnetic materials, and defects in a real AFI layer.
More realistic simulation which considers these factors is left as a future study.

\section{Summary}
\label{sec:summary}

We studied the temperature dependence of the amplitude of the SMR using the quantum Monte Carlo simulation.
Following Ref.~\onlinecite{Kato2020},
we formulated the spin conductance at the AFI/NM interface in terms of the retarded component of the local spin susceptibility and calculated it by the quantum Monte Carlo method employing an improved method which does not need direct numerical analytic continuation.
We showed that the SMR starts to grow below a specific temperature at which the correlation length reaches the system size.
As the temperature is lowered, the SMR once increases, takes a maximum value roughly at three-fifths of the ordering temperature, and then decreases.
We discussed the dependence of the magnitude of the spin and thickness of the AFI layer as well as effect of disordered exchange interactions.
These qualitative features are expected to be observed experimentally.

Our simulation for a simplified model is expected to provide a useful starting point for understanding of the SMR in quasi-two-dimensional quantum spin systems.
Detailed comparison with experiments using more realistic models, as well as detailed finite-size analysis with large-scale QMC simulation, is left for a future study.

\section*{Acknowledgements}
T. K. acknowledges support from the Japan Society for the Promotion of Science (JSPS KAKENHI Grants No. 20K03831). M. M. is financially supported by a Grant-in-Aid for Scientific Research B (20H01863) from MEXT, Japan. 
T. K. and T. I acknowledge K. Yoshimi and Y. Motoyama for giving us useful information on detailed implementation of QMC in the DSQSS package.

