The Heisenberg quantum antiferromagnet [@bethe; @anderson] is one of the most important models in solid state physics. A wide variety of physical behavior can be realized, depending on the geometry of the lattice over which the pair interaction $J{\bf S}_i \cdot {\bf S}_j$ is summed. On a 1D linear chain, the quantum fluctuations destroy the long range order, leading to a ground state with gapped (for integer $S$) or critical (for half odd-integer $S$) fluctuations [@luther; @haldane; @nums1]. On a 2D square lattice the ground state is ordered for all $S$[@affleck; @reger]. For $T > 0$ the order is destroyed by thermal fluctuations [@mermin]. In 3D, a second-order phase transition to an ordered state at $T_{\rm c} > 0$ is expected. The quantum fluctuations of the order parameter become negligible (relative to the classical, thermal fluctuations) in the neighborhood of $T_{\rm c}$ and the critical behavior should therefore be that of the classical 3D Heisenberg model. The quantum fluctuations vanish also in the limit $T \to \infty$, and hence their relative strength should have a maximum for $T_{\rm c} < T < \infty$. A universal “quantum critical” regime [@chn; @chubukov; @sandvik1] is in principle possible if the system is sufficiently close to a quantum critical ($T_c=0$) point [@sachdev].

Analytic calculations, employing several different approximations at finite $T$, have recently been carried out for the $S=1/2$ Heisenberg model on a 3D simple cubic lattice [@approx]. Critical temperatures in the range $T_{\rm c}/J = 0.89 - 1.13$ were obtained. Results for thermodynamic quantities depend strongly on the approximation used. High temperature series expansions have been used for the thermodynamics at high temperatures [@rushbrooke], but cannot be expected to be accurate close to $T_{\rm c}$. Quantum Monte Carlo (QMC) methods have been important for non-approximate calculations for 2D quantum antiferromagnets, but so far only limited results have been obtained for 3D systems [@qmc3d]. We here present detailed QMC calculations aimed at reliably extracting $T_{\rm c}$ as well as providing unbiased results for the temperature dependence of thermodynamic quantities. We also discuss the precise nature of the crossover from quantum mechanical to purely classical thermal order parameter fluctuations at the critical point. In particular, we show that the ratio $S({\bf Q})/[T\chi ({\bf Q})]$ of the static order parameter structure factor and susceptibility obeys a finite-size scaling at $T_{\rm c}$ that is [*completely universal*]{}, independent of the exponents of the phase transition or the dimensionality.

For the QMC simulations we have used the Stochastic Series Expansion (SSE) algorithm [@sse] (a generalization of Handscomb’s method [@handscomb]) for $L^3$ lattices with even $L$ up to 16 in the grand canonical ensemble. The method does not introduce any systematic errors in computed quantities. As a significant development allowing us to efficiently obtain results on a dense temperature grid, we have introduced “tempering” within the SSE algorithm. In this approach [@tempering], the configuration space of the simulation is extended to include a range of temperatures between which the system can make transitions, hence allowing calculations for a large number of temperatures much faster than with separate fixed-$T$ simulations. Also crucial in this work is the use of covariance effects to increase the numerical precision. The scheme introduced in allows for reduction of the statistical errors by almost three orders of magnitude for the 3D Heisenberg model.

In the path integral formulation, a quantum system at finite $T$ corresponds to a classical system with an extra “imaginary time” dimension of length $L_\tau \sim 1/T$ [@sondhi]. Standard finite-size scaling arguments [@cardy] can be applied to this system. In order to determine the critical temperature $T_{\rm c}$, we here first consider the $T$ and $L$ dependence of the spin stiffness $\rho_{\rm s}$. With the SSE method $\rho_s$ can be calculated directly from the global winding number fluctuations in periodic systems [@note1], as described in . Hyperscaling theory predicts the finite-size scaling of the stiffness at the critical point [@wallin], $$\rho_{\rm s}(T_{\rm c}) = L^{2-d-z},
\label{rhoscale}$$ where $d$ is the dimensionality and $z$ is the dynamic critical exponent which is zero for a $T_c > 0$ transition. For the 3D Heisenberg model, $L \rho_s$ graphed versus $T$ for different $L$ should hence intersect at $T_{\rm c}$. For $L=4-16$ the curves indeed intersect almost at a single point, as shown in Figure \[figrhos\]. For the smallest systems there is a clear tendency of the intersection points for sizes $L$ and $L+2$ to move to lower $T$ with increasing $L$, but for $L \ge 10$ no shifts can be seen within the statistical errors. From these results we therefore estimate $T_{\rm c} \approx 0.945J$.

=6.5cm

In the thermodynamic limit, the sublattice magnetization $|\langle S^z_i
\rangle |$ (the order parameter) increases continuously from zero for $T < T_{\rm c}$. Consider the static structure factor and susceptibility

$$\begin{aligned}
S({\bf q}) & = & {1\over L^3} \sum\limits_{i,j} 
{\rm e}^{i{\bf q} \cdot ({\bf r}_j - {\bf r}_i)} 
\langle S^z_j S^z_i \rangle,
\label{structq} \\
\chi({\bf q}) & = & {1\over L^3} \sum\limits_{i,j} 
{\rm e}^{i{\bf q} \cdot ({\bf r}_j - {\bf r}_i)} 
\int_0^\beta d\tau \langle S^z_j (\tau) S^z_i (0) \rangle . 
\label{suscq}\end{aligned}$$

\[sqxq\]

At the ordering wave-vector ${\bf Q}=(\pi,\pi,\pi)$, $S$ and $\chi$ should both scale with the system size as $L^{2-\eta}$ at $T=T_{\rm c}$, where $\eta \approx 0.03$ for the 3D Heisenberg transition [@peczak]. $\chi ({\bf Q})$ corresponds to the full space-time integral of the sublattice magnetization of the $3+1$ dimensional classical system and can therefore be expected to be less affected by corrections to scaling. Figure \[figscale\]a shows $\ln{(\chi /L^2)}$ vs. $\ln{(L)}$ at temperatures close to the $T_{\rm c}$ estimated above using the stiffness. Asymptotically we expect the points to fall on a straight line with slope $-\eta \approx -0.03$ if $T=T_{\rm c}$, and eventually diverge upward or downward below and above $T_{\rm c}$, respectively. At $T \approx 0.946J$ the $L \ge 10$ data show the expected scaling. Investigating the behavior for several temperatures in this neighborhood, and using the estimated accuracy of the exponent $\eta$ [@peczak], gives $T_{\rm c}/J = 0.946 
\pm 0.001$, in excellent agreement with the spin stiffness scaling.

With $T_{\rm c}$ determined, the expected scaling for $T > T_{\rm c}$ can be tested. In the thermodynamic limit $\chi$ should diverge as $t^{-\gamma}$, where $t=|T-T_{\rm c}|$ and $\gamma = \nu(2-\eta) 
\approx 1.40$ [@peczak]. Finite-size scaling predicts $\chi_L (t) = \chi_\infty (t) 
f[\xi (t)/L]$, where the correlation length $\xi \sim t^{-\nu}$. Hence, $\chi_L (t)t^\gamma$ graphed versus $Lt^\nu$ for different $L$ should collapse onto a single curve. As shown in Figure \[figscale\]b, the data are indeed collapsed, confirming the expected universality class and the estimated $T_{\rm c}$.

=6.8cm

We next consider the specific heat, $C={(\partial E / \partial T)/L^3}$. Figure \[figc\]a shows the temperature dependence for different system sizes. The position of the maximum of $C$ represents a size dependent definition of the critical point, and, as shown in more detail in Fig. \[figc\]b, seems to indicate that $T_{\rm c}$ is lower than the value estimated above. We believe that the reason for this apparent inconsistency is to be found in the fact that changing $T$ of the quantum system corresponds to changing $L_\tau$ of the corresponding 3+1 dimensional classical system. The calculated $C$ hence corresponds to $\partial E / \partial L_\tau^{-1}$. The critical behavior, governed by the exponent $\alpha$, of the actual classical specific heat and $C$ of the quantum system should nevertheless be the same (except exactly at a $T_{\rm c}=0$ quantum critical point) since the boundary between the ordered and disordered phases in the classical $(T,L_\tau)$ plane is analytic for $T > 0$. However, for small systems the shift of the peak position with $L$ is relatively large, and the associated $L_\tau$ variation can then cause significant deviations from the asymptotic finite-size scaling behavior of $C$. In order to partially compensate for the volume change due to the variable $L_\tau$, we consider the corresponding quantity defined per full space-time volume, i.e., $TC$. Figure \[figc\]c shows the behavior of the peak of $TC$. Compared to the behavior of $C$, we note that the peaks for these small systems are shifted noticeably to the right, and now appear to be consistent with $T_{\rm c} \approx 0.946J$.

=6.9cm

For the 3D classical Heisenberg transition $\alpha$ is negative ($\alpha \approx -0.11$) [@peczak], and hence $C$ should exhibit a cusp-like singularity maximum as $L \to \infty$. Finite-size scaling predicts the size dependence of the maximum value as $C_{\rm max} (L) = C_{\rm max} 
(\infty) - aL^{\alpha/\nu}$ ($\nu \approx 0.71$). The data for $TC$ is indeed well described by this expression for $L > 6$, and gives $T_{\rm c}C_{\rm max}(\infty) \approx 2.4$. The finite-size shift of the peak temperature should scale as $L^{-1/\nu}$. The results do not show this behavior. It is quite likely that this is due to the $L_\tau$ effect discussed above, which is not fully eliminated in $TC$. We note, however, that the expected scaling of $C$ is not easily observed even for the classical 3D Heisenberg model, for which significantly larger systems have been studied [@peczak]. A qualitative difference is that the peak shifts to higher $T$ with increasing $L$ in the classical case [@peczak]. As discussed above, for the quantum model the peak in $C$ initially moves with increasing $L$ from above $T_{\rm c}$ to slightly below $T_{\rm c}$. As $L \to \infty$, $C$ and $TC$ should both peak exactly at $T_{\rm c}$. One can therefore expect a change in the sign of the $C$ peak shift for some large $L$, and then the same asymptotic behavior as for the classical model.

We now discuss the nature of the crossover to purely classical (thermal) fluctuations of the order parameter as $T \to T_{\rm c}$ from above. For this purpose we note the well known sumrules $$\begin{aligned}
S({\bf q}) & = & {1\over \pi}\int_0^\infty d\omega 
(1+{\rm e}^{-\beta\omega}) S({\bf q},\omega), \\
\chi({\bf q}) & = & {2\over \pi}\int_0^\infty d\omega 
\omega^{-1} (1-{\rm e}^{-\beta\omega}) S({\bf q},\omega),\end{aligned}$$ where $S({\bf q},\omega)$ is the dynamic structure factor, and define the ratio $$R({\bf q}) = S({\bf q})/[T\chi({\bf q})].$$ For purely thermal fluctuations of momentum ${\bf q}$, all the spectral weight is in a $\delta$ function at $\omega=0$, and the classical relation $R({\bf q})=1$ is obeyed. Quantum $(|\omega| > 0)$ fluctuations increase $R$ above $1$. For the 3D antiferromagnet the fluctuations at the antiferromagnetic momentum ${\bf Q}$ should be purely classical for $T \le T_{\rm c}$. As $T \to \infty$, $R({\bf q}) \to 1$ for any ${\bf q}$. For $T \to T_{\rm c}$ from above, $L_\tau \ll \xi$ and in evaluating Eqs. (\[sqxq\]) the imaginary-time dependence of the long-distance correlation function can then be expanded to leading order in $\tau/ r$ (or $\tau^{1/z'}/ r$ if the system is not Lorenz invariant, which does not change the result). This gives $R({\bf Q})-1 \sim \xi^{-2} \sim t^{2\nu}$, and the size dependence at $T=T_{\rm c}$ should therefore be (substituting $\xi \to L$): $$R({\bf Q})-1 \sim L^{-2}.
\label{rscale}$$ Note that this scaling does not contain any model-dependent exponents, and remains valid in any number of dimensions (with the $S({\bf Q})$ and $\chi({\bf Q})$ corresponding to the relevant order parameter) as long as $T_{\rm c}>0$ and the order parameter does not commute with the Hamiltonian. (In the commuting case $R({\bf Q})=1$ for all $L$ and $T$). $R({\bf Q})$ should therefore in principle be useful, as an alternative to the stiffness scaling (\[rhoscale\]), for extracting $T_{\rm c}$ in cases where the critical exponents are not known.

Figure \[figratio\]a shows $R({\bf Q})-1$ versus temperature for 3D Heisenberg lattices. The finite-size gap leads to a $T \to 0$ divergence starting at a temperature which decreases with increasing $L$. For the larger systems there is a maximum at $T \approx 1.7J$, and a rapid drop below this temperature. The maximum value $R_{\rm max}-1 \approx 0.05$, which is about half of the value for the 2D Heisenberg model [@chubukov; @sandvik4], reflecting the reduction of quantum fluctuations with increasing dimensionality.

Figure. \[figratio\]b shows the $T$ dependence of $[R({\bf Q})-1]L^2$ for different system sizes. The curves indeed intersect. The intersection between $L$ and $L+2$ moves to lower $T$ with increasing $L$, and is consistent with the estimated $T_{\rm c} 
\approx 0.946J$. The proposed scaling law, Eq. (\[rscale\]), for $R({\bf Q})$ is hence confirmed.

=7.2cm

It is well known that a state of dominant classical order parameter fluctuations can extend to a finite regime well above $T_{\rm c}$ [@chn; @sachdev]. This is not in conflict with the universal scaling behavior of $R({\bf Q})$ noted here, according to which the classical value $R({\bf Q})=1$ always is approached algebraically as $T \to T_{\rm c}$. To see this, consider a system of 2D layers coupled by a very weak interlayer coupling $J_\perp$. As the temperature is lowered, the correlation length within the layers initially grows exponentially [@chn] but there are no significant correlations between the layers. For $T$ sufficiently low, but $T \gg T_{\rm c}$, they are in the “renormalized classical” [@chn] regime and $R({\bf Q}) \to 1$ exponentially [@sandvik4]. However, as $T$ is further lowered to the point, close to $T_{\rm c}$, where the layers start to correlate with each other, there will be a crossover to a $\xi \sim t^{-\nu}$ scaling behavior and the form (\[rscale\]) for $R({\bf Q})$. At this point $R({\bf Q})$ is already exponentially close to $1$, and the further algebraic reduction of $R({\bf Q})$ as $T \to T_{\rm c}$ is exponentially small.

In summary, we have determined the critical temperature $T_{\rm c}/J 
= 0.946 \pm 0.001$ of the 3D Heisenberg antiferromagnet with $S=1/2$, using QMC data of high precision. We have discussed the general nature of the crossover to purely classical thermal order parameter fluctuations as $T \to T_{\rm c}$ from above, and shown that asymptotically the transition is universal in the ratio $S({\bf Q})/[T\chi ({\bf Q})]$, independent of the universality class of the critical point. We have also noted how quantum mechanical effects can affect the finite-size scaling behavior of the specific heat.

I would like to thank D. K. Campbell for stimulating discussions and valuable comments on the manuscript. This work was supported by the NSF under grant No. DMR-97-12765. The numerical simulations were carried out at the NCSA.

H. A. Bethe, Z. Phys. [**21**]{}, 205 (1931).

P. W. Anderson, Phys. Rev. [**86**]{}, 694 (1952).

A. Luther and I. Peschel, Phys. Rev. B [**12**]{}, 3908 (1975).

F. D. M. Haldane, Phys. Lett. [**93A**]{}, 464 (1983); Phys. Rev. Lett. [**50**]{}, 1153 (1983).

K. Nomura, Phys. Rev. B [**40**]{}, 2421 (1989); S. Liang, Phys. Rev. Lett. [**64**]{}, 1597 (1990).

E. J. Neves and J. F. Peres, Phys. Lett. [**114A**]{}, 331 (1986); F. J. Dyson, E. H. Lieb, and B. Simon, J. Stat. Phys. [**18**]{}, 335 (1987); I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Commun. Math. Phys. [**155**]{}, 477 (1988).

J. D. Reger and A. P. Young, Phys. Rev. B [**37**]{}, 5978 (1988).

N. D. Mermin and H. Wagner, Phys. Rev. Lett. [**17**]{}, 1133 (1966).

S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. Lett. [**60**]{}, 1057 (1988); Phys. Rev. B [**39**]{}, 2344 (1989).

A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B [**49**]{}, 11919 (1994).

A. W. Sandvik and D. J. Scalapino, Phys. Rev. Lett. [**72**]{}, 2777 (1994).

S. Sachdev, Phys. Rev. B [**55**]{}, 142 (1997).

A. F. Barabanov and E. V. Zhasinas, Phys. Lett. A [**193**]{}, 191 (1994); R.-J. Liu and T.-L. Chen, Phys. Lett. A [**194**]{}, 137 (1994); A. Fladderjohann, K.-H. Mütter, and P. Wielath, Z. Phys. B [**100**]{}, 277 (1996).

G. S. Rushbrooke, G. A. Baker, and P. J. Wood in: [*Phase Transitions and Critical Phenomena*]{}, Vol. 3, ed. C. Domb and. M. S. Green (Academic, New York, 1974).

O. Nagai, Y. Yamada, and Y. Miyatake, in [*Quantum Monte Carlo Methods in Equilibrium and Nonequilibrium Systems*]{}, ed. M. Suzuki (Springer-Verlag Berlin, 1987); R. A. Sauerwein and M. J. De Oliviera, Mod. Phys. Lett. B [**9**]{}, 619 (1995); M. Troyer, M. E. Zhitomirsky, and K. Ueda, Phys. Rev. B [**55**]{}, R6117 (1997).

A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B [**43**]{}, 5950 (1991); A. W. Sandvik, J. Phys. A [**25**]{}, 3667 (1992).

D. C. Handscomb, Proc. Cambridge Philos. Soc. [**58**]{}, 594 (1962); [**60**]{}, 116 (1964).

E. Marinari and G. Parisi, Europhys. Lett. [**19**]{}, 451 (1992).

A. W. Sandvik, Phys. Rev. B [**54**]{}, 14910 (1996).

For a recent review, see: S. L. Sondhi, S. M. Girvin, J. P. Carini, and S. Shahar, Rev. Mod. Phys. [**69**]{}, 315 (1997).

For a review, see: [*Finite Size Scaling*]{}, ed. J. Cardy, (North-Holland, Amsterdam, 1988).

The Monte Carlo updates changing the winding number become inefficient in SSE simulations for $L > 16$, but for $L \le 16$ studied here they pose no problems.

A. W. Sandvik, Phys. Rev. B [**56**]{}, 11678 (1997).

M. Wallin, E. S. S[ø]{}rensen, S. M. Girvin, and A. P. Young, Phys. Rev. B [**49**]{}, 12115 (1994).

P. Peczak, A. M. Ferrenberg, and D. P. Landau, Phys. Rev. B [**43**]{}, 6087 (1991).

A. W. Sandvik, A. V. Chubukov, and S. Sachdev, Phys. Rev. B [**51**]{}, 16483 (1995); A. W. Sandvik and D. J. Scalapino, Phys. Rev. B [**53**]{}, R526 (1996). J. K. Kim and M. Troyer, Phys. Rev. Lett. [**80**]{}, 2705 (1998).
