Spin-wave thermodynamics of square-lattice antiferromagnets revisited

TL;DR

This paper introduces a novel modified spin-wave theory for square-lattice antiferromagnets that accurately describes thermodynamics across all temperatures without artificial phase transitions, validated by numerical and experimental comparisons.

Contribution

It develops a new spin-wave approach based on Wick decomposition that overcomes previous limitations and accurately reproduces experimental neutron scattering data.

Findings

Reliable thermodynamics from zero to high temperature

Agreement with quantum Monte Carlo and experimental data

Effective inclusion of higher-order interactions like ring exchange

Abstract

Modifying the conventional spin-wave theory in a novel manner based on the Wick decomposition, we present an elaborate thermodynamics of square-lattice quantum antiferromagnets. Our scheme is no longer accompanied by the notorious problem of an artificial transition to the paramagnetic state inherent in modified spin waves in the Hartree-Fock approximation. In the cases of spin $\frac{1}{2}$ and spin $1$, various modified-spin-wave findings for the internal energy, specific heat, static uniform susceptibility, and dynamic structure factor are not only numerically compared with quantum Monte Carlo calculations and Lanczos exact diagonalizations but also analytically expanded into low-temperature series. Modified spin waves interacting via the Wick decomposition provide reliable thermodynamics over the whole temperature range of absolute zero to infinity. Adding higher-order spin couplings such as ring exchange interaction to the naivest Heisenberg Hamiltonian, we precisely reproduce inelastic-neutron-scattering measurements of the high-temperature-superconductor-parent antiferromagnet $\mathrm{La}_2\mathrm{CuO}_4$. Modifying Dyson-Maleev bosons combined with auxiliary pseudofermions also yields thermodynamics of square-lattice antiferromagnets free from thermal breakdown, but it is less precise unless temperature is sufficiently low. Applying all the schemes to layered antiferromagnets as well, we discuss the advantages and disadvantages of modified spin-wave and combined boson-pseudofermion representations.

Tables

Table 1. Table 1: Evaluation of thermodynamics of square-lattice Heisenberg antiferromagnets obtained through various formulations from the points of view of whether it is precise, valid, and correct at low, intermediate, and high temperatures, respectively. We judge the high-temperature asymptotics of the TDD-MHFISW thermodynamics to be incorrect in spite of its correct high-temperature limit, because TDD MHFISWs artificially jump to the paramagnetic phase at a finite temperature.

	Precise low-$T$ analytics	Free from thermal breakdown	Correct high-$T$ asymptotics
TID MLSW	No	Yes	No
TID MWDISW	No	Yes	No
TID MHFISW	No	Yes	No
TDD MLSW	No	Yes	Yes
TDD MWDISW	Yes	Yes	Yes
TDD MHFISW	Yes	No	No
TDD MLSW+PF	No	Yes	No
TDD MWDISW+PF	Yes	Yes	No
TDD MHFISW+PF	Yes	Yes	No

Table 2. Table 2: Pseudotransition temperatures, T N ∗ superscript subscript 𝑇 N T_{\mathrm{N}}^{*} and/or T c ∗ superscript subscript 𝑇 c T_{\mathrm{c}}^{*} , in the TDD-MHFISW and TDD-MHFISW+PF theories for the Hamiltonian ( 126 ) of L → ∞ → 𝐿 L\rightarrow\infty with J ⟂ = 0.3 ​ J subscript 𝐽 perpendicular-to 0.3 𝐽 J_{\perp}=0.3J in the cases of S = 1 2 𝑆 1 2 S=\frac{1}{2} and S = 1 𝑆 1 S=1 . With temperature reaching T N ∗ superscript subscript 𝑇 N T_{\mathrm{N}}^{*} , the interlayer coupling ⟨ 𝒮 𝜹 l = 𝒛 ⟩ T subscript delimited-⟨⟩ subscript 𝒮 subscript 𝜹 𝑙 𝒛 𝑇 \langle\mathcal{S}_{\bm{\delta}_{l}=\bm{z}}\rangle_{T} jumps down to zero and so does the spontaneous magnetization M − z superscript subscript 𝑀 𝑧 M_{-}^{z} , and then at the temperature T c ∗ superscript subscript 𝑇 c T_{\mathrm{c}}^{*} , the intralayer couplings ⟨ 𝒮 𝜹 l = 𝒙 ⟩ T = ⟨ 𝒮 𝜹 l = 𝒚 ⟩ T subscript delimited-⟨⟩ subscript 𝒮 subscript 𝜹 𝑙 𝒙 𝑇 subscript delimited-⟨⟩ subscript 𝒮 subscript 𝜹 𝑙 𝒚 𝑇 \langle\mathcal{S}_{\bm{\delta}_{l}=\bm{x}}\rangle_{T}=\langle\mathcal{S}_{\bm{\delta}_{l}=\bm{y}}\rangle_{T} also jump down to zero to completely lose the correlation energy E 𝐸 E .

	TDD MHFISW			TDD MHFISW+PF
	$k_B{T}_N^{\ast }/JS$	$k_B{T}_c^{\ast }/JS$		$k_B{T}_N^{\ast }/JS$	$k_B{T}_c^{\ast }/JS$
$S=\frac{1}{2}$	$1.614$	$1.822$		$1.726$	—
$S=1$	$2.187$	$2.202$		$2.406$	—

Equations

$$\mathcal{H}=\sum _{i\in \mathrm{A}}\sum _{l=1}^z{J}_{{\mathbf{\delta }}_l}{\mathbf{S}}_{{\mathbf{r}}_i}\cdot {\mathbf{S}}_{{\mathbf{r}}_i+{\mathbf{\delta }}_l}=\sum _{j\in \mathrm{B}}\sum _{l=1}^z{J}_{{\mathbf{\delta }}_l}{\mathbf{S}}_{{\mathbf{r}}_j-{\mathbf{\delta }}_l}\cdot {\mathbf{S}}_{{\mathbf{r}}_j}.$$

$$\displaystyle\left\{\!\!\begin{array}[]{l}\displaystyle S_{\bm{r}_{i}}^{+}=\sqrt{2S}\left(1-\frac{a_{\bm{r}_{i}}^{\dagger}a_{\bm{r}_{i}}}{2S}\right)a_{\bm{r}_{i}}\\ S_{\bm{r}_{i}}^{-}=\sqrt{2S}a_{\bm{r}_{i}}^{\dagger}\\ S_{\bm{r}_{i}}^{z}=S-a_{\bm{r}_{i}}^{\dagger}a_{\bm{r}_{i}}\\ \end{array}\right.\!\!\!,$$

$$\displaystyle\left\{\!\!\begin{array}[]{l}\displaystyle S_{\bm{r}_{j}}^{+}=\sqrt{2S}b_{\bm{r}_{j}}^{\dagger}\left(1-\frac{b_{\bm{r}_{j}}^{\dagger}b_{\bm{r}_{j}}}{2S}\right)\\ S_{\bm{r}_{j}}^{-}=\sqrt{2S}b_{\bm{r}_{j}}\\ S_{\bm{r}_{j}}^{z}=b_{\bm{r}_{j}}^{\dagger}b_{\bm{r}_{j}}-S\\ \end{array}\right.$$

$$\mathcal{H}=\sum _{l=0}^2{\mathcal{H}}^{(l)};{\mathcal{H}}^{(2)}=-NzJ{S}^2,$$

$$\displaystyle\mathcal{H}^{(1)}=JS\sum_{<i,j>}(a_{\bm{r}_{i}}^{\dagger}a_{\bm{r}_{i}}+b_{\bm{r}_{j}}^{\dagger}b_{\bm{r}_{j}}+a_{\bm{r}_{i}}^{\dagger}b_{\bm{r}_{j}}^{\dagger}+a_{\bm{r}_{i}}b_{\bm{r}_{j}}),\$$

$${\mathcal{H}}^{(0)}=-\frac{J}{2}\sum _{<i,j>}{a}_{{\mathbf{r}}_i}^{†}(a_{{\mathbf{r}}_i}+b_{{\mathbf{r}}_j}^{†}{)}^2{b}_{{\mathbf{r}}_j}.$$

$${\mathcal{H}}^{(0)}\simeq NzJ⟨⟨S-\mathcal{S}⟩{⟩}^2-J⟨⟨S-\mathcal{S}⟩⟩$$

$$\times \sum _{<i,j>}(a_{{\mathbf{r}}_i}^{†}{a}_{{\mathbf{r}}_i}+b_{{\mathbf{r}}_j}^{†}{b}_{{\mathbf{r}}_j}+a_{{\mathbf{r}}_i}^{†}{b}_{{\mathbf{r}}_j}^{†}+a_{{\mathbf{r}}_i}{b}_{{\mathbf{r}}_j})\equiv {\mathcal{H}}_{\mathrm{BL}}^{(0)}$$

$$\mathcal{H}\simeq {\mathcal{H}}^{(2)}+{\mathcal{H}}^{(1)}+{\mathcal{H}}_{\mathrm{BL}}^{(0)}\equiv {\mathcal{H}}_{\mathrm{BL}},$$

$$⟨⟨\mathcal{S}⟩⟩\equiv S$$

$$-\frac{1}{2}⟨⟨a_{{\mathbf{r}}_i}^{†}{b}_{{\mathbf{r}}_i+{\mathbf{\delta }}_l}^{†}+a_{{\mathbf{r}}_i}{b}_{{\mathbf{r}}_i+{\mathbf{\delta }}_l}⟩⟩,$$

$${\tilde{\mathcal{H}}}_{\mathrm{BL}}\equiv {\mathcal{H}}_{\mathrm{BL}}+\mu {\mathcal{M}}_{-}^z$$

$$\displaystyle p\equiv\left\{\!\!\begin{array}[]{lr}\displaystyle 1-\frac{\mu}{\sum_{l=1}^{z}J_{\bm{\delta}_{l}}\langle\!\langle\mathcal{S}\rangle\!\rangle}=1-\frac{\mu}{zJ\langle\!\langle\mathcal{S}\rangle\!\rangle}&(\mathrm{TDD})\\ \displaystyle 1&(\mathrm{TID})\\ \end{array}\right.\!\!\!.$$

$${\gamma }_{{\mathbf{k}}_{\nu }}\equiv \frac{{\sum }_{l=1}^z{J}_{{\mathbf{\delta }}_l}⟨⟨\mathcal{S}⟩⟩e^{i{\mathbf{k}}_{\nu }\cdot {\mathbf{\delta }}_l}}{{\sum }_{l=1}^z{J}_{{\mathbf{\delta }}_l}⟨⟨\mathcal{S}⟩⟩}=\frac{1}{z}\sum _{l=1}^z{e}^{i{\mathbf{k}}_{\nu }\cdot {\mathbf{\delta }}_l},$$

$${\omega }_{{\mathbf{k}}_{\nu }}\equiv \frac{{\epsilon }_{{\mathbf{k}}_{\nu }}}{{\sum }_{l=1}^z{J}_{{\mathbf{\delta }}_l}⟨⟨\mathcal{S}⟩⟩}\equiv \sqrt{p^2-{\gamma }_{{\mathbf{k}}_{\nu }}^2},$$

$$\displaystyle\epsilon\equiv\frac{p}{2}-\frac{1}{2N}\sum_{\nu=1}^{N}\omega_{\bm{k}_{\nu}},\$$

$$\tau \equiv \frac{1}{2N}\sum _{\nu =1}^N\frac{p}{{\omega }_{{\mathbf{k}}_{\nu }}}-\frac{1}{2}.$$

$$a_{{\mathbf{k}}_{\nu }}^{†}=\frac{1}{\sqrt{N}}\sum _{i\in \mathrm{A}}{e}^{i{\mathbf{k}}_{\nu }\cdot {\mathbf{r}}_i}{a}_{{\mathbf{r}}_i}^{†},$$

$$b_{{\mathbf{k}}_{\nu }}=\frac{1}{\sqrt{N}}\sum _{j\in \mathrm{B}}{e}^{i{\mathbf{k}}_{\nu }\cdot {\mathbf{r}}_j}{b}_{{\mathbf{r}}_j}$$

$${\alpha }_{{\mathbf{k}}_{\nu }}^{+}=a_{{\mathbf{k}}_{\nu }}^{†}\sinh {\theta }_{{\mathbf{k}}_{\nu }}+b_{{\mathbf{k}}_{\nu }}\cosh {\theta }_{{\mathbf{k}}_{\nu }},$$

$${\alpha }_{{\mathbf{k}}_{\nu }}^{-}=a_{{\mathbf{k}}_{\nu }}\cosh {\theta }_{{\mathbf{k}}_{\nu }}+b_{{\mathbf{k}}_{\nu }}^{†}\sinh {\theta }_{{\mathbf{k}}_{\nu }};$$

$$\cosh 2{\theta }_{{\mathbf{k}}_{\nu }}=\frac{p}{{\omega }_{{\mathbf{k}}_{\nu }}},\sinh 2{\theta }_{{\mathbf{k}}_{\nu }}=\frac{{\gamma }_{{\mathbf{k}}_{\nu }}}{{\omega }_{{\mathbf{k}}_{\nu }}},$$

$${\tilde{\mathcal{H}}}_{\mathrm{BL}}=\sum _{l=0}^2{E}^{(l)}+\sum _{\nu =1}^N{\epsilon }_{{\mathbf{k}}_{\nu }}\sum _{\sigma =\pm }{\alpha }_{{\mathbf{k}}_{\nu }}^{\sigma †}{\alpha }_{{\mathbf{k}}_{\nu }}^{\sigma }+2\mu NS,$$

$$E^{(2)}$$

$$E^{(1)}$$

$$E^{(0)}$$

$$=NzJ\left[(S-⟨⟨\mathcal{S}⟩⟩{)}^2+2(S-⟨⟨\mathcal{S}⟩⟩)ϵ\right].$$

$$⟨{\alpha }_{{\mathbf{k}}_{\nu }}^{\sigma †}{\alpha }_{{\mathbf{k}}_{\nu }}^{\sigma }{⟩}_T$$

$$=\frac{1}{e^{{\epsilon }_{{\mathbf{k}}_{\nu }}/k_{\mathrm{B}}T}-1}\equiv {\bar{n}}_{{\mathbf{k}}_{\nu }}$$

$$⟨\mathcal{S}{⟩}_0^{\prime }$$

$$⟨\mathcal{S}{⟩}_0$$

$$⟨\mathcal{S}{⟩}_T$$

$$I_1\equiv \frac{1}{N}\sum _{\nu =1}^N{\omega }_{{\mathbf{k}}_{\nu }}{\bar{n}}_{{\mathbf{k}}_{\nu }},$$

$$I_2\equiv \frac{1}{N}\sum _{\nu =1}^N\frac{p}{{\omega }_{{\mathbf{k}}_{\nu }}}{\bar{n}}_{{\mathbf{k}}_{\nu }}$$

$$⟨{\mathcal{M}}_{\mathrm{A}}^x{⟩}_T=⟨{\mathcal{M}}_{\mathrm{B}}^x{⟩}_T=0,⟨{\mathcal{M}}_{\mathrm{A}}^y{⟩}_T=⟨{\mathcal{M}}_{\mathrm{B}}^y{⟩}_T=0,$$

$$⟨{\mathcal{M}}_{\mathrm{A}}^z{⟩}_T=-⟨{\mathcal{M}}_{\mathrm{B}}^z{⟩}_T=N(S-\tau -I_2),$$

$$⟨{\mathcal{M}}_{-}^z{⟩}_T=2N(S-\tau -I_2)=0.$$

$$⟨{\mathcal{H}}^{(2)}{⟩}_T=-N{S}^2\sum _{l=1}^z{J}_{{\mathbf{\delta }}_l}=-NzJ{S}^2,$$

$$\displaystyle\langle\mathcal{H}^{(1)}\rangle_{T}=-2NS\sum_{l=1}^{z}J_{\bm{\delta}_{l}}\bigl{(}\langle\mathcal{S}\rangle_{T}-S\bigr{)}$$

$$\displaystyle\qquad\quad\ =-2NzJS\bigl{(}\langle\mathcal{S}\rangle_{T}-S\bigr{)},$$

$$\displaystyle\langle\mathcal{H}^{(0)}\rangle_{T}=-N\sum_{l=1}^{z}J_{\bm{\delta}_{l}}\bigl{(}\langle\mathcal{S}\rangle_{T}-S\bigr{)}^{2}$$

Full text

I Introduction

Some decades ago, the earliest treatment of antiferromagnetic spin waves (SWs) at finite temperatures K568 was modified T1524 ; T2494 ; H4769 ; T5000 in an attempt to formulate thermodynamics of square-lattice Heisenberg antiferromagnets and thereby to interpret neutron scattering measurements on the high-temperature-superconductor-parent compound $\mathrm{La}_{2}\mathrm{CuO}_{4}$ consisting of quantum spins $S=\frac{1}{2}$. S1613 Diagonalizing a bosonic Hamiltonian with its sublattice magnetizations constrained to be zero, Takahashi T1524 ; T2494 gave a precise description of the thermal quantities at sufficiently low temperatures. Hirsch, Tang, and Lazzouni H4769 ; T5000 also demonstrated that such modified SWs (MSWs) well reproduce exact-diagonalization results for small clusters. Since we cannot calculate magnetic susceptibilities, whether uniform or staggered and whether static or dynamic, at finite temperatures within the conventional SW (CSW) theory, their findings opened up a new avenue for the study of thermodynamics. However, Takahashi’s MSW thermodynamics deteriorates with increasing temperature, not only failing to design a Schottky-like peak of the specific heat but even encountering an artificial phase transition of the first order to the trivial paramagnetic solution, T2494 ; N034714 similarly to the Schwinger-boson (SB) mean-field (MF) theory. A617 ; Y064426

In order to avoid thermal breakdown, Ohara and Yosida O2521 ; O3340 proposed another way of modifying CSWs, which consists of diagonalizing the Hamiltonian without any constraint but constructing the free energy under vanishing sublattice magnetizations. This scheme yields a peaked specific heat but spoils the otherwise excellent low-temperature findings. Its high-temperature findings are also unfortunate to deviate from the trivial paramagnetic behavior. The internal energy and uniform susceptibility per spin never approach zero and $(g\mu_{\mathrm{B}})^{2}S(S+1)/3k_{\mathrm{B}}T$, respectively. Takahashi’s MSWs and Auerbach-Arovas’ SBs are accompanied by the discontinuous transition but thereafter stabilized into the correct paramagnetic state, while Ohara-Yosida’s MSWs remain correlated even in the $T\rightarrow\infty$ limit.

In this context, there is a rather different approach to thermodynamics of layered magnets. Combining Dyson-Maleev bosons D1217 ; D1230 and auxiliary pseudofermions B351 to adjust the local Hilbert space dimension to the original spin degrees of freedom, Irkhin, Katanin, and Katsnelson I1082 erased the artificial phase transition without spoiling the successful bosonic description of purely two-dimensional Heisenberg magnets at sufficiently low temperatures. T2494 ; A316 ; S5028 ; Y3733 This formulation unfortunately fails to reproduce the paramagnetic behavior at high temperatures but gives a satisfactory description of layered systems in the truly critical region crossing their magnetic ordering temperatures. It is also noteworthy that a fully convincing description of the two-to-three dimensional crossover is available within the $1/N$ expansion of the $O(N)$ model, I12318 ; I379 i.e. the nonlinear sigma model generalized to $N$-component spins, rather than through the $1/S$ expansion of any SW Hamiltonian.

Under such circumstances, we revisit the SW thermodynamics of square-lattice antiferromagnets to find a better solution with particular emphasis on convenience for practical purposes. Is there something else within a simple spin-wave Hamiltonian that is reliable over the whole temperature range and applicable to various spins and interactions? Since the MSW scheme initiated by Takahashi T1524 ; T2494 and Hirsch et al. H4769 ; T5000 impose a constraint condition of zero staggered magnetization on SWs via a Bogoliubov transformation dependent on temperature, we refer to this way of modifying CSWs as a temperature-dependent-diagonalization (TDD)-MSW scheme. The MSW scheme proposed by Ohara and Yosida O2521 ; O3340 manipulates SWs under the same condition but leaves the CSW Hamiltonian as it is. Then we refer to this way of modifying CSWs as a temperature-independent-diagonalization (TID)-MSW scheme. Modifying sublattice bosons in a TDD manner but bringing them into interaction based on the Wick decomposition (WD) rather than by the Hartree-Fock (HF) approximation, we can retain the excellent low-temperature description and connect it naturally with the paramagnetic behavior without any thermal breakdown.

II Modified Spin-Wave Thermodynamics

We divide the square lattice into two sublattices, referred to as A and B, each containing $N\equiv L/2$ spins of magnitude $S$. We denote a vector connecting nearest neighbors $\bm{r}_{i}\ (i\in\mathrm{A})$ and $\bm{r}_{j}\ (j\in\mathrm{B})$ by $\bm{\delta}_{l}$ with $l$ running from $1$ to $z$ to write the Hamiltonian of our interest as

[TABLE]

$z$ is equal to $4$ and $J_{\bm{\delta}_{l}}$ are all set to $J\,(>0)$ unless otherwise noted. We employ the Dyson-Maleev bosons

[TABLE]

to rewrite the Hamiltonian (1) into

[TABLE]

We decompose the $O(S^{0})$ quartic Hamiltonian $\mathcal{H}^{(0)}$ into quadratic (bilinear in the end) terms

[TABLE]

to have a tractable SW Hamiltonian,

[TABLE]

where we introduce the multivalued double-angle-bracket notation applicable for various approximation schemes

[TABLE]

which we shall read as the quantum average in the Dyson-Maleev-boson vacuum $\langle\mathcal{S}\rangle_{0}^{\prime}$ for the linear SW (LSW) formalism, the quantum average in the magnon vacuum $\langle\mathcal{S}\rangle_{0}$ for the WD-based interacting SW (WDISW) formalism, or the temperature-$T$ thermal average $\langle\mathcal{S}\rangle_{T}$ for the HF-decomposition-based interacting SW (HFISW) formalism. Note that any average $\langle\!\langle\mathcal{S}\rangle\!\rangle$ is independent of the site indices $\bm{r}_{i}$ and $\bm{\delta}_{l}$ by virtue of translation and rotation symmetries.

Let $\mathcal{M}_{\mathrm{A}}^{\lambda}\equiv\sum_{i\in\mathrm{A}}S_{\bm{r}_{i}}^{\lambda}$, $\mathcal{M}_{\mathrm{B}}^{\lambda}\equiv\sum_{j\in\mathrm{B}}S_{\bm{r}_{j}}^{\lambda}$, and $\mathcal{M}_{\pm}^{\lambda}\equiv\mathcal{M}_{\mathrm{A}}^{\lambda}\pm\mathcal{M}_{\mathrm{B}}^{\lambda}$. Then the TDD-MSW theory reads as diagonalizing the effective quadratic Hamiltonian

[TABLE]

with such $\mu$ as to satisfy $\langle\mathcal{M}_{-}^{z}\rangle_{T}=0$. We introduce a key variable $p\equiv\sqrt{q^{2}+1}$ to design various MSWs,

[TABLE]

We further define some functions of $p$,

[TABLE]

Via the Fourier transformation

[TABLE]

and the Bogoliubov transformation

[TABLE]

we can diagonalize the effective Hamiltonian into

[TABLE]

where $E^{(2)}$ is the classical ground-state energy and $E^{(l\leq 1)}$ are its $O(S^{l})$ quantum corrections,

[TABLE]

Every SW thermodynamics can be formulated in terms of $\langle\!\langle\mathcal{S}\rangle\!\rangle$ and $p$, which indicate how the SWs are interacting and modified, respectively. The TDD-MSW thermal distribution function reads

[TABLE]

with $\varepsilon_{\bm{k}_{\nu}}$ containing $\langle\!\langle\mathcal{S}\rangle\!\rangle$. Every time we encounter $\langle\!\langle\mathcal{S}\rangle\!\rangle$, we read it according to the scheme of the time,

[TABLE]

where the sums

[TABLE]

still contain $\langle\!\langle\mathcal{S}\rangle\!\rangle$ to be self-consistently determined. The sublattice magnetizations read

[TABLE]

and therefore, the constraint condition is given as

[TABLE]

In modified HFISW (MHFISW) schemes, we solve the simultaneous equations (29) plus (33) for $\langle\mathcal{S}\rangle_{T}$ and $p$. In modified WDISW (MWDISW) schemes, we substitute Eq. (28) into Eq. (33) to obtain $p$ and then $\langle\mathcal{S}\rangle_{0}$. In modified LSW (MLSW) schemes, $\langle\mathcal{S}\rangle_{0}^{\prime}$ is a constant and therefore we have only to solve Eq. (33) for $p$. Employing the Bloch-De Dominicis theorem B459 to evaluate the thermal average of the quartic Hamiltonian, we can calculate the internal energy,

[TABLE]

Having in mind that

[TABLE]

the dynamic structure factors read

[TABLE]

The static structure factors are available from them,

[TABLE]

and so are the static uniform susceptibilities,

[TABLE]

We compare the TDD-MLSW and TDD-MHFISW calculations of the internal energy $E\equiv\langle\mathcal{H}\rangle_{T}$, specific heat $C\equiv\partial E/\partial T$, and uniform susceptibility $\chi\equiv\sum_{\lambda=x,y,z}\chi^{\lambda\lambda}/3$ with quantum Monte Carlo (QMC) calculations in Figs. 1 ($S=\frac{1}{2}$) and 2 ($S=1$), in the former of which Kim and Troyer’s QMC findings K2705 are also presented ($\times$). While the TDD-MLSW scheme succeeds in designing antiferromagnetic peaks of $C$ and $\chi$, it is far from precise at low temperatures, missing the $T\rightarrow 0$ values of $\chi$ considerably. We can improve these poor low-temperature findings by taking account of the SW interactions $\mathcal{H}^{(0)}$. While the thus-obtained TDD-MHFISW findings are highly precise at sufficiently low temperatures, they completely fail to reproduce the overall temperature dependences. The worst of them is that they are accompanied by an artificial phase transition of the first order to the trivial paramagnetic solution at a certain finite temperature. The specific heat jumps down to zero and the susceptibility switches to that of free spins $L(g\mu_{\mathrm{B}})^{2}S(S+1)/3k_{\mathrm{B}}T$ when $k_{\mathrm{B}}T/J$ reaches $0.9108$ and $2.202$ for $S=\frac{1}{2}$ and $S=1$, respectively, where $p$ diverges, while $\langle\mathcal{S}\rangle_{T}$ vanishes, satisfying $p\langle\mathcal{S}\rangle_{T}=(k_{\mathrm{B}}T/zJ)\mathrm{ln}(1+1/S)$.

The TID-MSW theory starts with diagonalizing the CSW Hamiltonian $\mathcal{H}_{\mathrm{BL}}\equiv\mathrm{lim}_{\mu\rightarrow 0}\widetilde{\mathcal{H}}_{\mathrm{BL}}$. O2521 ; O3340 ; Y14008 ; Y11033 When we denote the probability of $n$ antiferromagnons of mode $\sigma$ with momentum $\bm{k}_{\nu}$ emerging at temperature $T$ by $P_{\bm{k}_{\nu}}(n)$, the free energy reads

[TABLE]

$P_{\bm{k}_{\nu}}(n)$ is determined so as to minimize the effective free energy

[TABLE]

where $\mu_{\nu}$ is obtained from the normalization condition $\sum_{n=0}^{\infty}P_{\bm{k}_{\nu}}(n)=1$, while $\mu$ from the constraint condition (33). Then we have

[TABLE]

with the effective energy

[TABLE]

to yield the TID-MSW optimal distribution

[TABLE]

The TID-MSW calculations of $E$, $C$ and $\chi$ are also shown in Figs. 1 and 2. They are indeed free from thermal breakdown but in poor agreement with QMC calculations through the whole temperature range. The intermediate-temperature peak of $C$ is too high and the low-temperature increase of $\chi$ is too slow. The high-temperature asymptotics is also incorrect (Refer to Appendix).

In order to retain the TDD-MHFISW precise low-temperature findings by all means and connect them naturally with the correct high-temperature asymptotics, we return to the TDD modification scheme but design interacting SWs in a different manner from the HF approximation. A new treatment of the $O(S^{0})$ quartic Hamiltonian $\mathcal{H}^{(0)}$ consists of applying the Wick theorem based on the magnon operators $\alpha_{\bm{k}_{\nu}}^{\sigma\dagger}$ and $\alpha_{\bm{k}_{\nu}}^{\sigma}$ to it and neglecting the residual normal-ordered interaction $:\!\mathcal{H}^{(0)}\!:$. Then we have the bilinear Hamiltonian (11) with $\langle\!\langle\mathcal{S}\rangle\!\rangle$ read as the SW ground-state expectation value $\langle\mathcal{S}\rangle_{0}$ given by Eq. (28). We show in Figs. 1 and 2 the MWDISW calculations as well. Both TDD and TID MWDISWs describe the Schottky-like peak of $C$ much better than MHFISWs and their descriptions of $E$ and $C$ are exactly the same at sufficiently low temperatures [cf. Eq. (84)]. However, a significant difference is detected between the TDD and TID modification schemes in describing $\chi$. The low-temperature increase of $\chi$ is imprecisely described by TID MSWs but precisely reproduced by TDD MWDISWs as well as by TDD MHFISWs [cf. Eq. (96)]. TDD MHFISWs encounter an artificial phase transition to the paramagnetic solution at an intermediate temperature, while TDD MWDISWs are free from thermal breakdown. Unlike TID MSWs, TDD MSWs inherently hit the correct high-temperature limit (Refer to Appendix). TDD MHFISWs are unfortunate to artificially jump to the paramagnetic solution, but TDD MLSWs and TDD MWDISWs are so successful as to give correct high-temperature asymptotics (See Fig. 8). The TDD-MWDISW thermodynamics is precise at both low and high temperatures and free from any thermal breakdown. Let us inquire further into low-temperature MSW findings analytically.

III Low-Temperature Series Expansion

In order to convert the $\nu$ summations (30), (31), and (38) into $x$ integrations and thereby expand them into low-temperature series, we define a state-density function,

[TABLE]

For the square lattice in the thermodynamic limit, this reads

[TABLE]

with leading coefficients explicitly given as

[TABLE]

We define integral functions as

[TABLE]

If we put

[TABLE]

and assume that $t\ll p-q$, $f_{\xi}(v)$ for TDD MSWs become Bose-Einstein integral functions and are therefore expanded in powers of $v$ as

[TABLE]

for $\xi=1,2,3,\cdots$, while $f_{\xi}(v)$ for TID MSWs are similarly expanded as

[TABLE]

for $\xi=1,3,5,\cdots$.

Now we can expand the integrals (30), (31), and (38) in powers of $t$ and $v$ as

[TABLE]

Having in mind that $p=\sqrt{1+q^{2}}=\sqrt{1+(vt)^{2}}$ for TDD MSWs and $p=\sqrt{1+q^{2}}=1$ for TID MSWs, we solve Eq. (33) for $v$ in an iterative manner to have

[TABLE]

in both cases. Hence Eqs. (58) and (62) read

[TABLE]

Considering that $p=1+t^{2}O(e^{-1/t})$, $\tau=\tau_{p=1}+tO(e^{-1/t})$, $\epsilon=\epsilon_{p=1}+t^{2}O(e^{-1/t})$, and therefore

[TABLE]

the nonlinearity of $t$ as a function of $T$ is weak, if any,

[TABLE]

We eventually have

[TABLE]

There is little difference of $O(e^{-1/T})$ between the TDD-MSW and TID-MSW low-temperature series expansions of $E$. Within the first four terms, they are exactly the same not only as each other but also as the CSW one. As far as $E$ and therefore $C$ at low temperatures are concerned, it does not matter whether and how SWs are modified but does whether and how they are interacting. Note that the chemical potential $v$ has little effect of $O(e^{-1/T})$ on $I_{1}$. On the other hand, there is a serious difference of $O(T)$ between the TDD-MSW and TID-MSW low-temperature series expansions of $\chi$. While they converge to the same $T\rightarrow 0$ limit, the TID-MSW scheme underestimates the initial slope by a factor of two. The TDD-MWDISW and TDD-MHFISW findings are precise and exactly the same as each other within the first two terms. With further increasing temperature, the latter deviate from the former and end in the artificial phase transition to the paramagnetic solution. We should be reminded that CSWs can reproduce nothing about $\chi$. $I_{3}$ expanded in powers of $t$ and $v$ contains a term $\propto-t^{2}\mathrm{ln}v$, which diverges in the $v\rightarrow 0$ limit, i.e. within CSW theories, but stays finite by virtue of the constraint condition (33) in MSW theories.

IV Dynamic Structure Factor

In order to further demonstrate the quality and reliability of the TDD-MWDISW thermodynamics, we show in Fig. 3 its findings for the dynamic structure factor $S(\bm{q},\omega)\equiv\sum_{\lambda=x,y,z}S^{\lambda\lambda}(\bm{q},\omega)/3$ in comparison with exact calculations of the expression (36) expanded as a continued fraction. G2999 ; G11766 TDD MWDISWs give an excellent description especially of the $\delta$-function peaks at the one-magnon frequency $\varepsilon_{\bm{k}_{\nu}}$. There is no visible difference between the MSW and exact calculations of them at all in the case of $S=1$. There are two facts noteworthy in this context. One is that CSWs cannot reproduce anything about $S(\bm{q},\omega)$ even at $T=0$ unless $L\rightarrow\infty$ and the other is that the TDD-MWDISW and TDD-MHFISW schemes are equivalent at $T=0$. In CSW theories, $p=1$ and therefore the excitation energy $\varepsilon_{\bm{k}_{\nu}}\equiv zJ\langle\!\langle\mathcal{S}\rangle\!\rangle\omega_{\bm{k}_{\nu}}$ goes to zero at $\bm{k}_{\nu}=\bm{0}$. Then we cannot calculate the $\nu$ summation in Eq. (36), because $\mathrm{sinh}2\theta_{\bm{k}_{\nu}}$ and $\mathrm{cosh}2\theta_{\bm{k}_{\nu}}$ as well as $\bar{n}_{\bm{k}_{\nu}}$ are divergent at $\bm{k}_{\nu}=\bm{0}$. If we identify quantum averages at absolute zero, $\langle\ \rangle_{T=0}$, with those in the magnon vacuum, $\langle\ \rangle_{0}$, to set $\bar{n}_{\bm{k}}$’s all equal to zero and take the thermodynamic limit $L\rightarrow\infty$ to convert the intractable $\nu$ summation into a convergent integral, the ground-state structure factors are available within CSW theories, however the situation is still the same that $\chi\propto S(\bm{0})/T$ is not available there. The thus-calculated CSW findings for the ground-state dynamic structure factor are exactly the same as the $T=0$ MSW calculations at each approximation level, i.e. LSW, WDISW, or HFISW. Still, the fact remains that we can calculate no magnetic correlation of finite-sized square-lattice Heisenberg antiferromagnets in terms of SWs without modifying them in a TDD manner. Since TID MSWs share the same Bogoliubov transformation as CSWs, they are also of little use in this context.

$S(\bm{q},\omega)$ of the Hamiltonian (1) with $S=\frac{1}{2}$ has ever been calculated in terms of SBs at a MF level C239 and MHFISWs. T487 They yield exactly the same findings at every temperature. Their findings at $T=0$ are therefore exactly the same as the MWDISW calculations shown in Fig. 3. However, as was fully demonstrated in Figs. 1 and 2, MHFISWs are fragile and easy to break down with increasing temperature, while MWDISWs are robust and free from thermal breakdown.

V Comparison with Experiments

It is interesting to analyze inelastic-neutron-scattering (INS) measurements on the spin-$\frac{1}{2}$ square-lattice antiferromagnet $\mathrm{La}_{2}\mathrm{CuO}_{4}$ in terms of our MWDISW theory. This material is poorly fitted for the naivest Heisenberg model (1) but well describable with a higher-order spin Hamiltonian C5377 based on a strongly correlated single-band Hubbard model at half filling, T1289 ; M9753 ; D235130

[TABLE]

Allowing electrons to directly hop only between nearest-neighbor Cu sites and then denoting the hopping energy and on-site interaction by $t$ and $U$, respectively, we have

[TABLE]

In the expression (97), we regard $\mathcal{H}/U$ as second order in $t/U$ so as to reproduce the Heisenberg Hamiltonian with a nearest-neighbor-only coupling (1) and therefore $\mathcal{V}/U$ as fourth order in $t/U$.

Let $\bm{q}\equiv\bm{k}_{\mathrm{in}}-\bm{k}_{\mathrm{sc}}$ and $\hbar\omega\equiv E_{\mathrm{in}}-E_{\mathrm{sc}}$ with $\bm{k}_{\mathrm{in}}$ ($\bm{k}_{\mathrm{sc}}$) and $E_{\mathrm{in}}$ ($E_{\mathrm{sc}}$) being the initial (final) wavevector and energy of neutrons, respectively. The double-differential INS cross section, defining the probability of scattering an incident neutron beam into a particular energy range $dE$ and direction perpendicular to the surface area subtending the solid angle $d\varOmega$, is directly related to the imaginary parts of the dynamic susceptibilities $\chi^{\lambda\lambda^{\prime}}(\bm{q},\omega)$,

[TABLE]

where $(\gamma r_{\mathrm{e}}/2\mu_{\mathrm{B}})^{2}$ measures $72.65\times 10^{-3}\,\mathrm{barn}/\mu_{\mathrm{B}}^{2}$, $F(\bm{q})$ is the magnetic form factor defined as $\int e^{i\bm{q}\cdot\bm{r}}|\phi(\bm{r})|^{2}d\bm{r}$ with $\phi(\bm{r})$ being a Wannier orbital, and $e^{-2W(\bm{q})}$ is the Debye-Waller factor. We rewrite $\mathcal{D}_{\bm{\alpha}}$ and $\mathcal{Q}_{\bm{\alpha},\bm{\beta},\bm{\gamma}}$ in terms of the Dyson-Maleev bosons (9) and denote their $O(S^{l})$ terms by $\mathcal{D}_{\bm{\alpha}}^{(l)}$ and $\mathcal{Q}_{\bm{\alpha},\bm{\beta},\bm{\gamma}}^{(l)}$, respectively, similarly to Eq. (10). Decomposing the $O(S^{0})$ quartic, $O(S^{2})$ quartic, $O(S^{1})$ sextic, and $O(S^{0})$ octic interactions $\mathcal{D}_{\bm{\alpha}}^{(0)}$, $\mathcal{Q}_{\bm{\alpha},\bm{\beta},\bm{\gamma}}^{(2)}$, $\mathcal{Q}_{\bm{\alpha},\bm{\beta},\bm{\gamma}}^{(1)}$, and $\mathcal{Q}_{\bm{\alpha},\bm{\beta},\bm{\gamma}}^{(0)}$ all into quadratic terms through the Wick theorem, we evaluate Eq. (99) in terms of MWDISWs and interpret separate INS experiments on different samples of $\mathrm{La}_{2}\mathrm{CuO}_{4}$ performed by Headings et al. at $T=10\,\mathrm{K}$ H247001 and by Coldea et al. at $T=295\,\mathrm{K}$. C5377 While the Landé $g$ factor may depend on direction in a solid, W1082 here we take an isotropic $g$ and set it to $2$ for simplicity. L224511

Figure 4 shows the experimental and theoretical findings for the imaginary part of the dynamic susceptibility $\chi(\bm{q},\omega)\equiv\sum_{\lambda=x,y,z}\chi^{\lambda\lambda}(\bm{q},\omega)$ at $T=10\,\mathrm{K}$. Supposing the ring exchange interaction is considerably strong, $J_{\square}=0.41J$, and making a manual correction to the conventional LSW (CLSW) dispersion relation, Headings et al. H247001 gave a good guide to the one-magnon cross section [Fig. 4(a)]. Taking account of the $O(S^{0})$ quartic terms $\mathcal{D}_{\bm{\alpha}}^{(0)}$ in the Heisenberg interactions but discarding the $O(S^{1})$ sextic terms $\mathcal{Q}_{\bm{\alpha},\bm{\beta},\bm{\gamma}}^{(1)}$ and $O(S^{0})$ octic terms $\mathcal{Q}_{\bm{\alpha},\bm{\beta},\bm{\gamma}}^{(0)}$ in the ring exchange interactions, Katanin and Kampf K100403 demonstrated that conventional HFISWs (CHFISWs) can indeed reproduce Headings’ guiding line with $J=151.9\,\textrm{meV}$ and $J^{\prime}=J^{\prime\prime}=J_{\square}/9.6=3.80\,\textrm{meV}$. Their estimate sounds more convincing with a moderate ring exchange interaction, $J_{\square}=0.24J$, claiming that any orbital other than Cu $3d_{x^{2}-y^{2}}$ should further contribute to magnetic interactions. Our full calculation, including all the up-to-$O(S^{0})$ terms, can also yield a moderate ring exchange interaction, $J_{\square}=20J^{\prime}=20J^{\prime\prime}=0.22J$, within the fourth-order perturbation theory (98).

Figure 5 shows the experimental findings for the INS cross section $\sigma$ at room temperature in comparison with our theoretical findings for $\sigma$ and $\chi(\bm{q},\omega)$ at various temperatures. Coldea et al. C5377 report that the energy dispersion of magnetic excitations along the high-symmetry directions $\mathrm{X}\,(\pi/a,0)$ to $\Lambda\,(\pi/2a,\pi/2a)$ becomes less pronounced upon heating. While Figs. 4 and 5 look consistent with such observations, we should be reminded that they are separate observations of different samples. The two samples have a small but non-negligible difference of magnetic interaction and that is mainly why they exhibit a visibly different zone-boundary dispersion. In this context, we further note that different mechanisms may be responsible for the zone-boundary dispersion. Higher-order expansions in $t/U$ and $1/S$ have competing effects on the zone-boundary one-magnon energies. Within the LSW description of the naivest Heisenberg Hamiltonian (1), there is no dispersion along the zone boundary between $\mathrm{X}$ and $\Lambda$. Higher-order spin couplings such as ring exchange interaction contribute to raising the one-magnon energies at around $\mathrm{X}$ with respect to those at around $\Lambda$, K100403 as was demonstrated in Figs. 4 and 5. Without any such higher-order exchange coupling, higher-order perturbation corrections within the Hamiltonian (1) also make the one-magnon energies along the zone boundary dispersive, raising those at around $\Lambda$ with respect to those at around $\mathrm{X}$. S216003 ; U282 Such observations are indeed obtained by INS experiments on another spin-$\frac{1}{2}$ square-lattice antiferromagnet, $\mathrm{Cu}(\mathrm{DCOO})_{2}\cdot 4\mathrm{D}_{2}\mathrm{O}$, R037202 ; C15264 whose $t/U$ is relatively small with respect to that of $\mathrm{La}_{2}\mathrm{CuO}_{4}$ and hence its nearest-neighbor-only coupling of $J=6.3\,\mathrm{meV}$. C561 There are some other frustrated square-lattice antiferromagnets with dispersive one-magnon energies along the zone boundary. K852 ; T197201 Dependences of the zone-boundary dispersion on exchange coupling, spin quantum number, and temperature remain to be investigated from both experimental and theoretical points of view.

Magnetization measurements on carrier-free $\mathrm{La}_{2}\mathrm{CuO}_{4}$ (Ref. K7430, ) and $\mathrm{Cu}(\mathrm{DCOO})_{2}\cdot 4\mathrm{D}_{2}\mathrm{O}$ (Refs. C561, ,B883, ) reveal their Néel transitions at $325\,\mathrm{K}$ and $16.5\,\mathrm{K}$, respectively. Hence it follows that the experimental observations in Fig. 5 describe the characteristics of SWs in the vicinity of the Néel temperature $T_{\mathrm{N}}$. WDISWs overestimate $T_{\mathrm{N}}$ of layered antiferromagnets, as will be shown in Sec. VII. We can make higher-order perturbation corrections to WDISWs in an attempt to reduce their overestimation of $T_{\mathrm{N}}$ and further interpret various experimental observations of layered antiferromagnets in the truly critical temperature region near $T_{\mathrm{N}}$. Such fluctuation corrections, together with auxiliary pseudofermions, indeed improve the HFISW calculations of sublattice magnetizations in a spin-$1$ layered perovskite with easy-axis single-ion anisotropy, $\mathrm{K}_{2}\mathrm{NiF}_{4}$. I1082 In three dimensions, however, the present MSW theories all reduce to a CSW formulation below $T_{\mathrm{N}}$ with their chemical potential vanishing, $\mu\rightarrow-0\,(T\rightarrow T_{\mathrm{N}}+0)$. Therefore, we take more interest in developing an efficient MSW theory in lower dimensions. The TDD-MWDISW thermodynamics of square-lattice antiferromagnets is precise and analytic at low temperatures and remains reliable at high temperatures. We expect that fluctuation corrections will further improve it at intermediate temperatures rather than attempt to refine that of layered antiferromagnets.

VI Modified Spin Waves Combined with Pseudofermions

A combined boson-pseudofermion representation of spin operators can also give a reasonable description of thermodynamic properties. Tuning this tool in various aspects, Irkhin, Katanin, and Katsnelson I1082 reproduced magnetization measurements on a layered antiferromagnet with particular emphasis on the dimensional-crossover temperature region. It must be of benefit to our future study to compare MSW theories of current interest with what they call self-consistent SW (SCSW) theories. Irkhin-Katanin-Katsnelson’s SCSWs within one-particle picture are obtained by combining TDD MHFISWs with pseudofermions. Besides them, various MSWs combined with pseudofermions (MSWs+PFs) are available to formulate thermodynamics. Before concluding our study, we investigate boson-pseudofermion mixed languages in detail. Indeed combining TDD MSWs with pseudofermions results in preventing them from thermal breakdown, but the resultant findings are not necessarily superior to those of pure TDD MSWs. TDD MLSWs, MHFISWs, and MWDISWs combined with pseudofermions (MLSWs+PFs, MHFISWs+PFs, and MWDISWs+PFs) all fail to reproduce the high-temperature paramagnetic behavior correctly. They remain correlated to underestimate the total spin degrees of freedom even at sufficiently high temperatures. Irkhin-Katanin-Katsnelson’s SCSWs work better in the vicinity of magnetic ordering than elsewhere and are therefore suitable for describing two-dimensional magnets with interlayer coupling and/or magnetic anisotropy, I1082 ; G024427 which we shall demonstrate in the final section, in comparison with purely bosonic TDD-MSW calculations.

The Bar’yakhtar-Krivoruchko-Yablonskiĭ representation of spin operators B351 reads

[TABLE]

where auxiliary pseudofermions, created on sites $\bm{r}_{i}\ (i\in\mathrm{A})$ and $\bm{r}_{j}\ (j\in\mathrm{B})$ by $c_{\bm{r}_{i}}^{\dagger}$ and $d_{\bm{r}_{j}}^{\dagger}$, respectively, are employed to bring Dyson-Maleev bosons, created by $a_{\bm{r}_{i}}^{\dagger}$ and $b_{\bm{r}_{j}}^{\dagger}$, into kinematic interaction. The local Hilbert space in which the Bar’yakhtar-Krivoruchko-Yablonskiĭ transforms (107) operate is spanned by the basis vectors

[TABLE]

Applying $S_{\bm{r}_{l}}^{z}$ to each basis set shows that $|n_{l},0\rangle_{\bm{r}_{l}}$ with $n_{l}$ larger than $2S$ and $|n_{l},1\rangle_{\bm{r}_{l}}$ are nonphysical states. The Bar’yakhtar-Krivoruchko-Yablonskiĭ bosons combined with pseudofermions (107) rewrite the Hamiltonian (1) into

[TABLE]

Similarly to Eq. (11), we decompose the quartic terms to have a bilinear Hamiltonian,

[TABLE]

where the multivalued double-angle-bracket notation is applied to a boson-pseudofermion operator

[TABLE]

In an attempt to formulate thermodynamics in terms of MSWs+PFs, we again introduce the effective quadratic Hamiltonian $\widetilde{\mathcal{H}}_{\mathrm{BL}}\equiv\mathcal{H}_{\mathrm{BL}}+\mu\mathcal{M}_{-}^{z}$, similarly to Eq. (14), where $\mu$ is determined so as to satisfy

[TABLE]

In taking every Bar’yakhtar-Krivoruchko-Yablonskiĭ thermal average, the “pseudoprojection” operator $\mathcal{P}$ serves to cancel nonphysical-state contributions. The MSW+PF effective Hamiltonian is diagonalized into

[TABLE]

where the MSW energies $\varepsilon_{\bm{k}_{\nu}}$ look the same as those in Eq. (24), while pseudofermions, occupying only nonphysical states and therefore making no physical sense, are immobile to have flat bands above the MSW dispersive bands,

[TABLE]

In terms of the Bar’yakhtar-Krivoruchko-Yablonskiĭ thermal averages of quasiparticles

[TABLE]

the order parameter $\langle\!\langle\mathcal{S}\rangle\!\rangle$ is explicitly written as

[TABLE]

The sublattice magnetizations read

[TABLE]

and therefore, the constraint condition is given as

[TABLE]

We choose one from Eqs. (119)–(121) and solve it simultaneously with Eq. (123) for $\langle\!\langle\mathcal{S}\rangle\!\rangle$ and/or $p$, where every time we encounter $\langle\!\langle\mathcal{S}\rangle\!\rangle$, we read it as one of Eqs. (119)–(121) according to the scheme of the time. The internal energy and uniform susceptibilities are expressed as

[TABLE]

We compare the thus-calculated $E$, $C$, and $\chi$ of the purely two-dimensional square-lattice Heisenberg antiferromagnets (1) with the TDD-MSW calculations in Fig. 6. MHFISWs and MWDISWs, whether combined with pseudofermions or not, bring exactly the same results at sufficiently low temperatures, namely, within the first three (up to $T^{5}$) terms of $E$ and the first two (up to $T^{2}$) terms of $\chi T$. With further increasing temperature, there occurs a difference between their calculations. The MHFISW findings increase too rapidly as functions of temperature and encounter an artificial transition to the paramagnetic solution without auxiliary pseudofermions. Their first-order transition is fictitious indeed, but they hit the correct high-temperature limit. From this point of view, MHFISWs+PFs unfortunately mishit the high-temperature paramagnetic behavior (Refer to Appendix), even though they are free from thermal breakdown. In the TDD-MHFISW+PF formulation, every thermodynamic quantity is overestimated and underestimated at intermediate and high temperatures, respectively. The internal energy and uniform susceptibility per spin never approach zero and $(g\mu_{\mathrm{B}})^{2}S(S+1)/3k_{\mathrm{B}}T$, respectively, even when temperature rises. This behavior is reminiscent of the TID-MSW findings.

Now we are in a position to have a bird’s-eye view of various SW theories designed for thermodynamics of square-lattice antiferromagnets. Whether they are useful will depend on how we use them. If we take interest only in the low-temperature properties, TDD MSWs, no matter how they are interacting and whether or not they are combined with PFs, may serve our purpose. If we are under the necessity of reproducing the overall temperature dependences of experimental findings, TDD MHFISWs are no longer useful at all but TDD MWDISWs remain reliable. Indeed TDD MHFISWs are inferior in describing the thermodynamic properties of low-dimensional antiferromagnets, but they are potentially superior in predicting the ordering temperature of three-dimensional layered antiferromagnets (cf. Fig. 7 and see Ref. I1082, as well). While TID MSWs appear to have no particular merit, there is an example of their successful detection Y074703 of novel nuclear spin relaxation mediated by multimagnons P398 ; B359 in an intertwining double-chain ferrimagnet. It is hard to say in a word which scheme is the best. When we investigate two-dimensional antiferromagnets with the aim to quantitatively analyze their thermodynamic properties as functions of temperature, Table 1 serves to evaluate various MSW representations.

VII Summary and Discussion

TDD MWDISWs serve as a sophisticated language for extensive quantum antiferromagnets in two dimensions. Their thermodynamics is highly precise at low temperatures, sufficiently robust against thermalization, and well applicable to frustrated antiferromagnets as well, covering both static and dynamic properties. The spin-$\frac{1}{2}$ square-lattice Heisenberg antiferromagnet (1) frustrated by diagonal nearest-neighbor couplings have also ever been studied in terms of MSWs by numerous authors with particular emphasis on its nature at absolute zero. Chandra and Doucot C9335 discussed within the CLSW scheme that a spin liquid phase may exist in the sense of the conventional antiferromagnetic long-range order disappearing. Hirsch and Tang H2887 employed the TDD-MLSW scheme in an attempt to reveal a transition to a disordered phase. Nishimori and Saika N4454 followed to introduce the TDD-MHFISW scheme and claim that a Néel-type long-range order consisting of either two or four sublattices should survive the strong frustration. While there are analytic I8206 and numerical Y2384 ; W1350021 MSW calculations of static properties at finite temperatures as well, the TDD-MWDISW thermodynamics of this model is never yet revealed. Such calculation is not only intriguing in itself but may also give a piece of information otherwise unavailable.

We take a particular interest in developing a SW thermodynamics unattainable within the CSW theory. We therefore devote much effort to describing low-dimensional magnets with no long-range order at finite temperatures. In three dimensions, a long-range order survives thermalization more or less, during which every MSW scheme reduces to a CSW formulation with its chemical potential unavailable from the constraint condition (33). R119 In lower dimensions, SWs may be modified over the whole temperature range to have their otherwise unavailable thermodynamics. However, real layered magnets encounter an order-disorder phase transition at a finite temperature, which we shall denote by $T_{\mathrm{N}}$, due to weak interlayer coupling and/or magnetic anisotropy. In an attempt to interpret the spontaneous magnetization as an anomalous function of temperature in a three-dimensional ferrimagnet, Karchev K325219 ; K216003 proposed modifying CSWs with a temperature- and sublattice-dependent constraint condition. Focusing on a narrow critical region near $T_{\mathrm{N}}$, Irkhin, Katanin, and Katsnelson I1082 applied a MSW theory refined with auxiliary pseudofermions to quasi-two-dimensional ferromagnets and antiferromagnets. There is also a possibility of further developing our scheme along this line.

In order to gain better understanding of SW thermodynamics, we finally compare the TDD-MSW and TDD-MSW+PF findings for the order-disorder phase transitions of layered Heisenberg antiferromagnets describable within the Hamiltonian (1) of $z=6$. When we set $J_{\bm{\delta}_{l}}$ to $J\,(>0)$ and $J_{\perp}\,(>0)$ for $1\leq l\leq 4$ and $5\leq l\leq 6$, respectively, the Hamiltonian of current interest reads

[TABLE]

In terms of the Dyson-Maleev bosons (9) and Bar’yakhtar-Krivoruchko-Yablonskiĭ bosons combined with pseudofermions (107), this Hamiltonian is also rewritten into (10) and (109) and then decomposed into the quadratic forms (12) and (110), respectively. Under the constraint condition of zero staggered magnetization, the MSW and MSW+PF effective quadratic Hamiltonians $\widetilde{\mathcal{H}}_{\mathrm{BL}}\equiv\mathcal{H}_{\mathrm{BL}}+\mu\mathcal{M}_{-}^{z}$ are diagonalized into (24) and (113), respectively. The bond-dependent order parameters $\langle\!\langle\mathcal{S}\rangle\!\rangle$ are explicitly written as

[TABLE]

with the bond-dependent sums and their averages over bonds

[TABLE]

When we set $J$ and $J_{\perp}$ equal, Eqs. (127)–(129) reduce to Eqs. (27)–(29) and Eqs. (119)–(121) with $\eta^{\prime}$ and $H^{\prime}$ coinciding with $\eta$ and $H$, respectively.

We compare the TDD-MSW and TDD-MSW+PF calculations of the staggered magnetization $M_{-}^{z}\equiv\langle\mathcal{M}_{-}^{z}\rangle_{T}$ and internal energy $E\equiv\langle\mathcal{H}\rangle_{T}$ in Fig. 7. Similarly to the case with single-layer antiferromagnets (Fig. 6), the TDD-MHFISW calculations of layered antiferromagnets are accompanied by an artificial phase transition of the first order to the trivial paramagnetic solution. Every TDD-MHFISW formulation inevitably suffers a pseudotransition from correlated to uncorrelated spins at a finite temperature, which we shall denote by $T_{\mathrm{c}}^{*}$ distinguishably from the real scenario $T_{\mathrm{c}}=\infty$. However, the present findings require more careful observation and need to be distinguished from what we have observed in Fig. 6. There are two artificial phase transitions of the first order in the TDD-MHFISW thermodynamics of layered antiferromagnets (cf. Table 2), unless $J=J_{\perp}$. Preceding the total breakdown of the system into free spins at $T=T_{\mathrm{c}}^{*}$, there occurs another fictitious phase transition of the first order at a temperature below $T_{\mathrm{c}}^{*}$ which we shall denote by $T_{\mathrm{N}}^{*}$, where the spontaneous magnetization $\langle\mathcal{M}_{-}^{z}\rangle_{T}$ disappears, while the exchange interaction $\langle\mathcal{H}\rangle_{T}$ stays finite, losing the interlayer coupling $\langle\mathcal{S}_{\bm{\delta}_{l}=\bm{z}}\rangle_{T}$ but maintaining the intralayer correlation $\langle\mathcal{S}_{\bm{\delta}_{l}=\bm{x}}\rangle_{T}=\langle\mathcal{S}_{\bm{\delta}_{l}=\bm{y}}\rangle_{T}$. First arises a three-to-two-dimensional phase transition at $T=T_{\mathrm{N}}^{*}$ and then an antiferromagnetic-to-paramagnetic phase transition at $T=T_{\mathrm{c}}^{*}$ on TDD MHFISWs. That is why we find exactly the same discontinuous transition to occur at exactly the same temperature in Figs. 6 and 7. A layered antiferromagnet, once its interlayer coupling is lost, degenerates into planar antiferromagnets merely stacked up.

While TDD MHFISWs get free from a total breakdown with the help of auxiliary pseudofermions, yet they suffer a partial breakdown (cf. Table 2). Even in the TDD-MHFISW+PF scheme, a layered antiferromagnet is still misled to discontinuously decouple into planar antiferromagnets at $T=T_{\mathrm{N}}^{*}$ with its spontaneous magnetization $\langle\mathcal{M}_{-}^{z}\rangle_{T}$ jumping down to zero. In order to overcome this difficulty, Irkhin, Katanin, and Katsnelson I1082 considered replacing the bond-dependent short-range order parameter (129) by the averaged one

[TABLE]

The thus-tuned MHFISWs+PFs are free from any fictitious transition but misread the ground-state properties (cf. insets in Fig. 7). With $\langle\mathcal{S}\rangle_{T}$ averaged over $\bm{\delta}_{l}$’s, MHFISWs+PFs have the same Bogoliubov transformation (23) as MLSWs at $T=0$ and therefore give the same ground-state energy and magnetization as MLSWs. Averaging $\langle\mathcal{S}\rangle_{T}$ spoils the otherwise precise low-temperature findings of MHFISWs fully demonstrated in Figs. 1 and 2. In addition, every MSW+PF formulation, whether it employs $\langle\mathcal{S}\rangle_{T}$ as they are or averages them, misreads the high-temperature paramagnetic behavior with everlasting two-dimensional correlation $\langle\mathcal{S}_{\bm{\delta}_{l}=\bm{x}}\rangle_{T}=\langle\mathcal{S}_{\bm{\delta}_{l}=\bm{y}}\rangle_{T}$.

Indeed the MHFISW+PF thermodynamics formulated with the effective order parameter (134) is not so precise as the MWDISW one at temperatures far below and above the ordering temperature $T_{\mathrm{N}}$, but it is highly successful in the truly critical temperature region near $T_{\mathrm{N}}$. It reproduces the emergent magnetization $\langle\mathcal{M}_{-}^{z}\rangle_{T}$ much more reasonably than the MWDISW thermodynamics and can be further improved by taking account of fluctuation corrections, within a random-phase approximation for instance. I1082 Irkhin and Katanin I12318 ; I379 further demonstrated that the $O(N)$ model can more suitably describe the truly critical temperature region. While the MWDISW thermodynamics is less precise there, yet there may be a possibility of applying it to layered ferrimagnets, where an anomalous temperature dependence of the spontaneous magnetization far below $T_{\mathrm{N}}$ is also an interesting topic and remains to be solved. K325219 ; K216003 Though we have here focused our attention on planar antiferromagnets aiming to solve a longstanding problem lying in MSW theories for them, the MWDISW scheme can be applied to many other intriguing magnetic systems including even “zero-dimensional” clusters Y157603 ; C280 to reveal their overall thermodynamic properties possibly with precise expressions at low temperatures.

Acknowledgements.

We are grateful to J. Ohara for useful comments. This study was supported by the Ministry of Education, Culture, Sports, Science, and Technology of Japan.

APPENDIX A HIGH-TEMPERATURE ASYMPTOTICS

We investigate high-temperature asymptotic thermodynamics of square-lattice Heisenberg antiferromagnets obtained through various MSW formulations. We employ the effective temperature $t\equiv{k_{\mathrm{B}}T}/zJ\langle\!\langle\mathcal{S}\rangle\!\rangle$ (55).

First, we consider TDD MSWs. Having in mind that $-1\leq\gamma_{\bm{k}_{\nu}}\leq 1$, we find that the thermal distribution function (26) no longer depends on the wavevector in the high-temperature limit,

[TABLE]

With converging $\bar{n}_{\bm{k}_{\nu}}$, i.e. diverging $p$, in the $t\rightarrow\infty$ limit, we have

[TABLE]

Then the constraint condition (33) becomes

[TABLE]

to consistently yield diverging $p$ in the $t\rightarrow\infty$ limit,

[TABLE]

Equations (29), (34), and (38) read

[TABLE]

We learn that TDD MSWs hit the correct high-temperature limit. TDD MHFISWs end up with an artificial transition to the paramagnetic solution, but TDD MLSWs and TDD MWDISWs really give correct high-temperature asymptotics, as is demonstrated in Fig. 8.

Next, we consider TID MSWs. Having in mind that $p=1$ as well as $-1\leq\gamma_{\bm{k}_{\nu}}\leq 1$, we find that the thermal distribution function (43) still depends on the wavevector in the high-temperature limit,

[TABLE]

making it difficult to analytically calculate the high-temperature asymptotics. $\epsilon$, $\tau$, and $\eta$ are not dependent on temperature. With converging $\bar{n}_{\bm{k}_{\nu}}$ in the $t\rightarrow\infty$ limit, we have

[TABLE]

We numerically solve the constraint condition (33) in the high-temperature limit

[TABLE]

to find $\lim_{t\rightarrow\infty}(-\mu/tzJ\langle\!\langle\mathcal{S}\rangle\!\rangle)$ to be nonzero and thus $\lim_{t\rightarrow\infty}\bar{n}_{\bm{k}_{\nu}}$ to be finite. Equations (29), (34), and (38) read

[TABLE]

We find $\lim_{t\rightarrow\infty}(-\mu/tzJ\langle\!\langle\mathcal{S}\rangle\!\rangle)$ to be $1.31972709$ and $0.74682827$, $\lim_{t\rightarrow\infty}E/LJ$ to be $-0.21475012$ and $-0.37814423$, and $\lim_{t\rightarrow\infty}\chi k_{\mathrm{B}}T/L(g\mu_{\mathrm{B}})^{2}$ to be $0.11970015$ and $0.42202662$ for $S=\frac{1}{2}$ and $S=1$, respectively (See Fig. 8). We learn that TID MSWs cannot give correct high-temperature asymptotics.

Finally, we consider TDD MSWs combined with PFs. The Bar’yakhtar-Krivoruchko-Yablonskiĭ thermal average of MSWs (115) in the high-temperature limit looks the same as (135). That of PFs (118) is free from wavevector dependence at every temperature,

[TABLE]

If the calculations of the high-temperature limit (136)–(142) remain unchanged, the constraint condition (123) becomes

[TABLE]

with the Brillouin function

[TABLE]

Since $B_{S}(x)$ is a monotonically increasing function, $\lim_{t\rightarrow\infty}p/t=0$ is one and only real solution of Eq. (156), and therefore, the Bar’yakhtar-Krivoruchko-Yablonskiĭ thermal averages (135) and (155) both diverge in the $t\rightarrow\infty$ limit. Then, the calculations (139)–(142) are no longer reliable, because every thermal-distribution-weighted integration over $\bm{k}$ is not necessarily exchangeable with the $t\rightarrow\infty$ operation. Yet applying Eqs. (139)–(142) to Eqs. (121), (124), and (125) yields

[TABLE]

which are indeed inconsistent with the correct calculations demonstrated in Fig. 8.

In order to analyze the high-temperature asymptotics of the TDD-MSW+PF thermodynamics correctly, we expand every thermal quantity into high-temperature series. Employing the state-density function (44) and its explicit expression (45) again, the thermal-distribution-weighted sums read

[TABLE]

At low temperatures, having in mind that low-energy excitations, i.e. $x$ much smaller than $p-q$, make major contributions to the integrals (160)–(162), we expand $w(x)$ in powers of $x$ as Eqs. (45) and (46), while at high temperatures of current interest, having in mind that

[TABLE]

we expand thermal distribution functions as

[TABLE]

Note further that

[TABLE]

Supposing $p$ to be of order $\sqrt{t}$, we consider such high temperatures as to satisfy $1\ll p\ll t$. Via the expansions (164) with $\beta\equiv(p/t)\sqrt{1-(y/p)^{2}}$, the integrals (160)–(162) become

[TABLE]

where the $p$ dependences of $\epsilon$ and $\tau$ read

[TABLE]

Putting $\beta\equiv(2S+1)p/t$ in Eq. (164) yields the high-temperature series expansions of $-\bar{f}$ and $-\bar{f}(1-\bar{f})$. The constraint condition (123) becomes

[TABLE]

to reveal that $p$ indeed diverges as $\sqrt{t}$ at high temperatures,

[TABLE]

Equations (121) and (125) read

[TABLE]

Having in mind that

[TABLE]

we eventually have

[TABLE]

which are now consistent with the numerical calculations demonstrated in Fig. 8. We learn that TDD MSWs+PFs cannot give correct high-temperature asymptotics.