Quantum paramagnetism and helimagnetic orders in the Heisenberg model on the body centered cubic lattice

TL;DR

This study explores the quantum and classical magnetic phases of the Heisenberg model on a BCC lattice, revealing quantum fluctuations significantly alter classical phase boundaries and suppress nonmagnetic phases for spin-1/2.

Contribution

It provides a comprehensive analysis of quantum effects on complex magnetic orders in a BCC lattice, highlighting the impact of quantum fluctuations on phase stability and ordering vectors.

Findings

Quantum fluctuations shift phase boundaries and helimagnetic pitches.

Nonmagnetic phase exists for S=1/2 but disappears at S=1.

Calculated ordering temperatures agree with experimental data.

Abstract

We investigate the spin $S=1/2$ Heisenberg model on the body centered cubic lattice in the presence of ferromagnetic and antiferromagnetic nearest-neighbor $J_{1}$, second-neighbor $J_{2}$, and third-neighbor $J_{3}$ exchange interactions. The classical ground state phase diagram obtained by a Luttinger-Tisza analysis is shown to host six different (noncollinear) helimagnetic orders in addition to ferromagnetic, N\'eel, stripe and planar antiferromagnetic orders. Employing the pseudofermion functional renormalization group (PFFRG) method for quantum spins ($S=1/2$) we find an extended nonmagnetic region, and significant shifts to the classical phase boundaries and helimagnetic pitch vectors caused by quantum fluctuations while no new long-range dipolar magnetic orders are stabilized. The nonmagnetic phase is found to disappear for $S=1$. We calculate the magnetic ordering temperatures from PFFRG and quantum Monte Carlo methods, and make comparisons to available data

Tables

Table 1. Table 1: The pitch vectors of the helimagnetic orders in the classical J 1 subscript 𝐽 1 J_{1} – J 2 subscript 𝐽 2 J_{2} – J 3 subscript 𝐽 3 J_{3} Heisenberg model on the BCC lattice together with their degeneracy in the first Brillouin zone and the energy per spin in the corresponding ground state.

Pitch vector ($𝐪$)	Component $q$	Degeneracy of $𝐪$ in the first Brillouin zone	Energy $\frac{E}{N{S}^2}$
$(q,0,0)$	$2{\cos }^{-1}\left(\frac{\left|J_1\right|}{J_2+4J_3}\right)$	6–fold	$J_2-2J_3-\frac{2{\left(J_1\right)}^2}{J_2+4J_3}$
$(2\pi -q,0,0)$
$(q,q,0)$	${\cos }^{-1}\left(\frac{\left|J_1\right|-J_2-2J_3}{2J_3}\right)$	12–fold	$-J_2-2J_3-\frac{{\left(-\left|J_1\right|+J_2\right)}^2}{2J_3}$
$(2\pi -q,q,0)$		24–fold
$(q,q,q)$	$2{\cos }^{-1}\left(\frac{\left|J_1\right|+\sqrt{{\left(J_1\right)}^2-32J_3\left(J_2-4J_3\right)}}{16J_3}\right)$	8–fold	$-\frac{{\left(J_1\right)}^2\left(u+16J_3\left(4J_3+J_2\right)\right)+\left|J_1\right|u^{\frac{3}{2}}+384{\left(J_2\right)}^2{\left(J_3\right)}^2}{1024{\left(J_3\right)}^3}$
$(2\pi -q,q,q)$		24–fold	$u={\left(J_1\right)}^2+32J_3\left(4J_3-J_2\right)$
$(0,0,0)$		1–fold	$-4\left|J_1\right|+3J_2+6J_3$
$(2\pi ,0,0)$
$(\pi ,\pi ,0)$		6–fold	$-J_2-2J_3$
$(\pi ,\pi ,\pi )$		1–fold	$-3J_2+6J_3$

Table 2. Table 2: For the classical phase diagram of the J 1 subscript 𝐽 1 J_{1} – J 2 subscript 𝐽 2 J_{2} – J 3 subscript 𝐽 3 J_{3} Heisenberg model on the BCC lattice, we provide the analytical expression for the different phase boundaries between the members of column Phase I and Phase II, and the order of the corresponding phase transitions.

Phase I	Phase II	Equation for the phase boundary	Type of Phase Transition
$(0,0,0)$	$(q,0,0)$	$J_2=\left|J_1\right|-4J_3$	2nd Order
$(2\pi ,0,0)$	$(2\pi -q,0,0)$
$(0,0,0)$	$(q,q,q)$	$J_2=\left|J_1\right|-4J_3$	2nd Order
$(2\pi ,0,0)$	$(2\pi -q,q,q)$
$(0,0,0)$	$(\pi ,\pi ,\pi )$	$J_2=\frac{2}{3}\left|J_1\right|$	1st Order
$(2\pi ,0,0)$
$(q,0,0)$	$(q,q,0)$	$J_2=0$	1st Order
$(2\pi -q,0,0)$	$(2\pi -q,q,0)$
$(q,q,0)$	$(\pi ,\pi ,0)$	$J_2=\left|J_1\right|$	2nd Order
$(2\pi -q,q,0)$
$(q,q,0)$	$(q,q,q)$	$J_2=\frac{56\left|J_1\right|{\left(J_3\right)}^2-11{\left(J_1\right)}^2{J}_3-32{\left(J_3\right)}^3-\sqrt{\left|J_1\right|{\left(J_3\right)}^2{\left(5\left|J_1\right|-16J_3\right)}^3}}{8{\left(J_3\right)}^2}$	1st Order
$(2\pi -q,q,0)$	$(2\pi -q,q,q)$
$(q,q,q)$	$(\pi ,\pi ,0)$	$J_2=\frac{{\left(J_1\right)}^4+4416{\left(J_1\right)}^2{\left(J_3\right)}^2+\left({\left(J_1\right)}^2+48{\left(J_3\right)}^2\right)\sqrt[3]{u}+{\left(u\right)}^{2/3}+36864{\left(J_3\right)}^4}{108J_3\sqrt[3]{u}}$	1st Order
$(2\pi -q,q,q)$
$(q,q,q)$	$(\pi ,\pi ,\pi )$	$J_2=\frac{{\left(J_1\right)}^2+144{\left(J_3\right)}^2}{36J_3}$	1st Order
$(2\pi -q,q,q)$
with $u={\left(J_1\right)}^6-10872{\left(J_1\right)}^4{\left(J_3\right)}^2-2709504{\left(J_1\right)}^2{\left(J_3\right)}^4+108\sqrt{3}{J}_3\sqrt{-{\left(J_1\right)}^2\left(12{\left(J_3\right)}^2-{\left(J_1\right)}^2\right){\left(512{\left(J_3\right)}^2-{\left(J_1\right)}^2\right)}^3}+7077888{\left(J_3\right)}^6$

Table 3. Table 3: The critical value J 2 c / | J 1 | superscript subscript 𝐽 2 𝑐 subscript 𝐽 1 J_{2}^{c}/|J_{1}| of the transition between the FM/Néel and the stripe order obtained from PFFRG and compared to different methods for the S = 1 / 2 𝑆 1 2 S=1/2 J 1 subscript 𝐽 1 J_{1} – J 2 subscript 𝐽 2 J_{2} Heisenberg model on the BCC lattice with J 3 = 0 subscript 𝐽 3 0 J_{3}=0 . The results marked with an asterisk are from the present study.

	Phase I	Phase II	Method	$J_2^c/\left|J_1\right|$
$J_1$–Ferromagnetic	$(0,0,0)$	$(\pi ,\pi ,\pi )$	PFFRG∗	$0.56\left(2\right)$
			Exact Diagonalization∗	$0.568$
			Coupled Cluster Method∗	$0.579$
			Rotation-invariant Green’s function method Müller, Patrick et al. (2015)	$0.68$
			Random phase approximation Tahir-Kheli and Jarrett (1964)	$0.6799$
$J_1$–Antiferromagnetic	$(2\pi ,0,0)$	$(\pi ,\pi ,\pi )$	PFFRG∗	$0.70\left(2\right)$
			Coupled Cluster Method Farnell et al. (2016)	$0.704$
			Exact Diagonalization Schmidt et al. (2002)	$0.7$
			Non-linear spin-wave theory Majumdar and Datta (2009)	$0.705$
			Random phase approximation Pantić et al. (2014)	$0.72$
			Linked Cluster Series expansions Oitmaa and Zheng (2004a)	$0.705\left(5\right)$

Table 4. Table 4: Critical temperatures T c / | J 1 | subscript 𝑇 𝑐 subscript 𝐽 1 T_{c}/|J_{1}| of the S = 1 / 2 𝑆 1 2 S=1/2 J 1 subscript 𝐽 1 J_{1} – J 2 subscript 𝐽 2 J_{2} BCC Heisenberg model for various values of the frustrating AF J 2 subscript 𝐽 2 J_{2} coupling, obtained from different methods for both FM as well as AF J 1 subscript 𝐽 1 J_{1} . The results marked with an asterisk are from the present study. For plotted data see Fig.  12 .

	Method	$\frac{J_2}{\left|J_1\right|}=0$	$\frac{J_2}{\left|J_1\right|}=0.1$	$\frac{J_2}{\left|J_1\right|}=0.2$	$\frac{J_2}{\left|J_1\right|}=0.3$	$\frac{J_2}{\left|J_1\right|}=0.4$	$\frac{J_2}{\left|J_1\right|}=0.5$	$\frac{J_2}{\left|J_1\right|}=0.6$
$J_1$–FM	PFFRG (one loop)∗	$1.45\left(1\right)$	$1.33\left(1\right)$	$1.18\left(1\right)$	$1.04\left(1\right)$	$0.88\left(1\right)$	$0.66\left(1\right)$
	PFFRG (two loop)∗	$1.37\left(1\right)$	$1.26\left(3\right)$	$1.11\left(1\right)$	$0.97\left(1\right)$	$0.84\left(3\right)$	$0.62\left(1\right)$
	QMC∗	$1.260\left(1\right)$
	HTE [${\chi }_u$] Rushbrooke et al. (1964); Oitmaa and Bornilla (1996); Oitmaa and Zheng (2004b)	$1.2602\left(5\right)$
	HTE [${\chi }_u^{\mathrm{quot}}$] Müller, Patrick et al. (2015)	$1.253\left(8\right)$	$1.125\left(15\right)$		$1.06\left(32\right)$
	HTE [${\chi }_u^{\mathrm{DA}}$] Müller, Patrick et al. (2015)	$1.268\left(8\right)$	$1.16\left(3\right)$		$0.89\left(7\right)$
	HTE [$S\left(𝑸\right)$] Richter et al. (2015)	$1.273\left(8\right)$	$1.17\left(3\right)$	$1.09\left(10\right)$	$0.89\left(7\right)$	$0.72\left(11\right)$	$0.45\left(18\right)$
	GFA Müller, Patrick et al. (2015)	$1.359$	$1.247$	$1.136$	$1.022$	$0.903$	$0.771$
$J_1$–AF	PFFRG (one loop)∗	$1.63\left(1\right)$	$1.52\left(1\right)$	$1.36\left(1\right)$	$1.26\left(2\right)$	$1.10\left(1\right)$	$0.90\left(1\right)$	$0.71\left(1\right)$
	PFFRG (two loop)∗	$1.50\left(1\right)$	$1.42\left(1\right)$	$1.32\left(1\right)$	$1.16\left(1\right)$	$1.06\left(1\right)$	$0.89\left(1\right)$	$0.75\left(1\right)$
	QMC∗	$1.377\left(2\right)$
	HTE [${\chi }_s$] Oitmaa and Zheng (2004b)	$1.376\left(4\right)$
	HTE $\left[S\left(𝑸\right)\right]$ Richter et al. (2015)	$1.50\left(8\right)$	$1.36\left(10\right)$	$1.26\left(13\right)$	$1.09\left(13\right)$	$0.96\left(7\right)$	$0.75\left(6\right)$	$0.61\left(10\right)$
	HTE [${\chi }_s$] Oitmaa and Zheng (2004a)	$1.38\left(2\right)$	$1.26\left(2\right)$	$1.13\left(1\right)$	$1.00\left(2\right)$	$0.86\left(2\right)$	$0.71\left(2\right)$	$0.57\left(3\right)$
	GFA Richter et al. (2015)	$1.530$	$1.369$	$1.195$	$1.004$	$0.786$	$0.520$

Table 5. Table 5: For the nearest-neighbor Heisenberg FM and the AF on the BCC and simple cubic (SC) lattices, we provide for S = 1 / 2 𝑆 1 2 S=1/2 , the Curie T C S = 1 / 2 / [ | J 1 | ​ S ​ ( S + 1 ) ] superscript subscript 𝑇 𝐶 𝑆 1 2 delimited-[] subscript 𝐽 1 𝑆 𝑆 1 T_{C}^{S=1/2}/[|J_{1}|S(S+1)] temperature (column 2) and the Néel temperature T N S = 1 / 2 / [ J 1 ​ S ​ ( S + 1 ) ] superscript subscript 𝑇 𝑁 𝑆 1 2 delimited-[] subscript 𝐽 1 𝑆 𝑆 1 T_{N}^{S=1/2}/[J_{1}S(S+1)] (column 3) as obtained from Quantum Monte Carlo. We also provide the ordering temperature T C / N S → ∞ / [ | J 1 | ​ S ​ ( S + 1 ) ] superscript subscript 𝑇 𝐶 𝑁 → 𝑆 delimited-[] subscript 𝐽 1 𝑆 𝑆 1 T_{C/N}^{S\to\infty}/[|J_{1}|S(S+1)] for the corresponding classical Heisenberg model as obtained from Classical Monte Carlo. In the classical model, the equality of the ordering temperatures for FM exchange and AF exchange is a consequence of the fact that the free energy is an even function of the coupling J 1 subscript 𝐽 1 J_{1}   Oitmaa and Zheng ( 2004b ) .

Lattice	$T_C^{S=1/2}/\left[\left|J_1\right|S\left(S+1\right)\right]$	$T_N^{S=1/2}/\left[J_{1}S\left(S+1\right)\right]$	$T_{C/N}^{S\to \infty }/\left[\left|J_1\right|S\left(S+1\right)\right]$	$\frac{T_C^{S=1/2}}{T_{C/N}^{S\to \infty }}$	$\frac{T_N^{S=1/2}}{T_{C/N}^{S\to \infty }}$
BCC	$1.680\left(1\right)$	$1.836\left(3\right)$	$\mathrm{2.054\hspace{0.17em}241}\left(52\right)$ Chen et al. (1993)	$0.818\left(1\right)$	$0.894\left(3\right)$
SC	$1.119\left(1\right)$ Troyer et al. (2004); Wessel (2010)	$1.261\left(1\right)$ Sandvik (1998)	$\mathrm{1.442\hspace{0.17em}929}\left(77\right)$ Chen et al. (1993)	$0.776\left(1\right)$	$0.874\left(1\right)$

Table 6. Table 6: Results of fits to Eq. ( 11 ). An absent k max subscript 𝑘 max k_{\rm max} indicates a fit without including scaling corrections. In the fits indicated with ∗ we fix ω = 0.8 𝜔 0.8 \omega=0.8 .

	$R=U_4$,  $U_4^{\mathrm{AF}}$							$R=\xi /L$,  ${\xi }^{\mathrm{AF}}/L$
	$L_{\min }$	$m_{\max }$	$k_{\max }$	$R^{\ast }$	$T_c$	${\chi }^2/\text{DOF}$		$L_{\min }$	$m_{\max }$	$k_{\max }$	$R^{\ast }$	$T_c$	${\chi }^2/\text{DOF}$
	$8$	$1$		$0.4984\left(1\right)$	$1.2571\left(1\right)$	$310/48$		$8$	$1$		$0.6199\left(4\right)$	$1.26094\left(7\right)$	$101/48$
	$12$	$1$		$0.4960\left(2\right)$	$1.2583\left(1\right)$	$96/35$		$12$	$1$		$0.6228\left(6\right)$	$1.26058\left(9\right)$	$49/35$
	$16$	$1$		$0.4940\left(4\right)$	$1.2592\left(2\right)$	$44/22$		$16$	$1$		$0.623\left(1\right)$	$1.2606\left(2\right)$	$22/22$
	$8$	$2$		$0.4986\left(2\right)$	$1.2575\left(1\right)$	$190/47$		$8$	$1$	$0^{\ast }$	$0.636\left(3\right)$	$1.2599\left(2\right)$	$64/47$
FM	$12$	$2$		$0.4965\left(2\right)$	$1.2583\left(2\right)$	$40/34$
	$16$	$2$		$0.4951\left(4\right)$	$1.2589\left(2\right)$	$10/21$
	$8$	$2$	$0^{}$	$0.4912\left(6\right)$	$1.2598\left(2\right)$	$30/45$
	$8$	$2$	$0^{\ast }$	$0.486\left(1\right)$	$1.2605\left(3\right)$	$30/46$
	$8$	$1$		$0.4965\left(2\right)$	$1.3752\left(1\right)$	$156/43$		$8$	$1$		$0.5843\left(5\right)$	$1.38212\left(8\right)$	$1067/43$
	$12$	$1$		$0.4949\left(2\right)$	$1.3760\left(2\right)$	$64/32$		$12$	$1$		$0.5986\left(6\right)$	$1.3801\left(1\right)$	$122/32$
	$16$	$1$		$0.4940\left(4\right)$	$1.3764\left(2\right)$	$44/21$		$16$	$1$		$0.607\left(1\right)$	$1.3791\left(2\right)$	$19/21$
	$8$	$2$		$0.4966\left(2\right)$	$1.3754\left(1\right)$	$91/42$		$8$	$1$	$0^{}$	$0.68\left(7\right)$	$1.375\left(2\right)$	$52/41$
AF	$12$	$2$		$0.4953\left(2\right)$	$1.3760\left(2\right)$	$27/31$		$8$	$1$	$0^{\ast }$	$0.672\left(3\right)$	$1.3756\left(2\right)$	$52/42$
	$16$	$2$		$0.4949\left(4\right)$	$1.3762\left(2\right)$	$17/20$
	$8$	$2$	$0^{}$	$0.4947\left(3\right)$	$1.3762\left(2\right)$	$31/40$
	$8$	$2$	$0^{\ast }$	$0.489\left(1\right)$	$1.3773\left(3\right)$	$34/41$

Equations

$$\hat{\mathcal{H}}=J_1\sum _{⟨i,j{⟩}_1}{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j+J_2\sum _{⟨i,j{⟩}_2}{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j+J_3\sum _{⟨i,j{⟩}_3}{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j,$$

$$\mathcal{H}=\sum _{i,j}J(|{\mathbf{R}}_i-{\mathbf{R}}_j|){\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$\mathcal{H}=\sum _{\mathbf{k}}J(\mathbf{k})\mathbf{S}(\mathbf{k})\cdot \mathbf{S}(-\mathbf{k}).$$

$$\mathbf{S}(\mathbf{k})=\frac{1}{\sqrt{N}}\sum _i{e}^{-ı\mathbf{k}\cdot {\mathbf{R}}_i}{\mathbf{S}}_i.$$

$${\mathbf{S}}_i=R(\pm \cos (\mathbf{k}\cdot {\mathbf{R}}_i),\sin (\mathbf{k}\cdot {\mathbf{R}}_i),0),$$

$$\displaystyle J(\mathbf{k})=8J_{1}\cos\bigg{(}\frac{k_{x}}{2}\bigg{)}\cos\bigg{(}\frac{k_{y}}{2}\bigg{)}\cos\bigg{(}\frac{k_{z}}{2}\bigg{)}+$$

$$+$$

$$C_{ij}=C_{ij}^{zz}={\int }_0^{\infty }d\tau ⟨{\hat{S}}_i^z(\tau ){\hat{S}}_j^z(0)⟩,$$

$$\begin{array}{rl}{\mathbf{S}}_i=S\{& \mathbf{a}\cos [(\pi ,\pi ,0)\cdot {\mathbf{R}}_i]+\mathbf{b}\cos [(\pi ,0,\pi )\cdot {\mathbf{R}}_i]+\mathbf{c}\cos [(0,\pi ,\pi )\cdot {\mathbf{R}}_i]\\ & +\mathbf{d}\cos [(\pi ,-\pi ,0)\cdot {\mathbf{R}}_i]+\mathbf{e}\cos [(-\pi ,0,\pi )\cdot {\mathbf{R}}_i]+\mathbf{f}\cos [(0,-\pi ,\pi )\cdot {\mathbf{R}}_i]\}\end{array}$$

$$\begin{array}{rl}{\mathbf{a}}^2+{\mathbf{b}}^2+{\mathbf{c}}^2+{\mathbf{d}}^2+{\mathbf{e}}^2+{\mathbf{f}}^2& =1\\ \mathbf{a}\cdot \mathbf{b}+\mathbf{d}\cdot \mathbf{e}=\mathbf{a}\cdot \mathbf{c}+\mathbf{d}\cdot \mathbf{f}=\mathbf{a}\cdot \mathbf{e}+\mathbf{b}\cdot \mathbf{d}=\mathbf{a}\cdot \mathbf{f}+\mathbf{c}\cdot \mathbf{d}=\mathbf{b}\cdot \mathbf{c}+\mathbf{e}\cdot \mathbf{f}& =\mathbf{c}\cdot \mathbf{e}+\mathbf{b}\cdot \mathbf{f}=\mathbf{b}\cdot \mathbf{e}+\mathbf{a}\cdot \mathbf{d}+\mathbf{c}\cdot \mathbf{f}=0\end{array}$$

$$R(T,L)=R^{\ast }+\sum _{m=1}^{m_{\max }}{a}_m(T-T_c{)}^m{L}^{m/\nu },$$

$$\begin{array}{rl}R(T,L)=R^{\ast }& +\sum _{m=1}^{m_{\max }}{a}_m(T-T_c{)}^m{L}^{m/\nu }\\ & +L^{-\omega }\sum _{k=0}^{k_{\max }}{b}_k(T-T_c{)}^k{L}^{k/\nu },\end{array}$$

$$U_4^{\ast }=U_4^{\mathrm{AF}\ast }=0.491(5),$$

$$(\xi /L{)}^{\ast }=({\xi }^{\mathrm{AF}}/L{)}^{\ast }=0.62(2),$$

Full text

I Introduction

The long-range ferromagnetic (FM) or antiferromagnetic (AF) order of spins pinned to the sites of a bipartite crystal lattice becomes frustrated in the presence of longer range AF interactions, a scenario called parametric frustration. For Heisenberg spins in the classical limit, i.e., spin $S\to\infty$, these competing interactions provide a promising route towards realizing (noncollinear) helimagnetic orders, i.e., spiral spin structures Villain (1959); Yoshimori (1959); Nagamiya (1968); Rastelli et al. (1979). On the square lattice, the FM or AF ordering of spins when frustrated via AF second- and third neighbor Heisenberg interactions is known to stabilize one- and two-dimensional helimagnetic orders Rastelli et al. (1979). When the reciprocal spin $1/S$ becomes nonzero, quantum fluctuations enter the picture, and their amplitude increases with increasing reciprocal spin $1/S$. In fact, the quantity $1/S$ plays the same role for quantum fluctuations as the temperature does for classical fluctuations Kaganov and Chubukov (1987), although they may act differently as has been suggested in the kagome Heisenberg AF Harris et al. (1992); Reimers and Berlinsky (1993); Huse and Rutenberg (1992); Henley (2009); Korshunov (2002); Chernyshev and Zhitomirsky (2014); Götze and Richter (2015). In the semiclassical ($1/S\ll 1$) regime, it is known that for collinear phases the quantum corrections to the ground state and the spin wave spectrum are modest Kaganov and Chubukov (1987). On the other hand, in helimagnets, owing to the delicate interplay of competing interactions the impact of quantum fluctuations is likely to be of significance. It was shown by Chubukov Chubukov (1984) that quantum fluctuations lead to a shift of the spiral pitch vector value, but keep the two Goldstone modes ($\mathbf{k}=\mathbf{0}$ and $\mathbf{k}=\pm\mathbf{Q}$) intact, thus preserving the general structure of the magnon spectrum. However, in the small spin–$S$ limit where strong quantum fluctuations are at play, the fate of the helimagnetic ground states remains largely not studied. In particular, it is of interest to investigate whether in the extreme quantum limit of $S=1/2$, quantum fluctuations could melt the helimagnetic structures Iqbal et al. (2017, 2018, 2019) and potentially realize a quantum paramagnetic ground state Anderson (1973); Balents (2010); Savary and Balents (2016). In this context, low-dimensional quantum spin systems have traditionally attracted much attention due to the significant increase in the role played by quantum fluctuation effects. On the square lattice, for $S=1/2$, the helimagnetic orders give way to a quantum paramagnet over an appreciable region in parameter space Sindzingre et al. (2009, 2010); Iqbal et al. (2016a), however, similar scenarios in three-dimensional lattices remain largely unexplored.

In this paper, employing the pseudofermion functional renormalization group (PFFRG) method Reuther and Wölfle (2010), we address the question as to which degree the impact of quantum fluctuations is mellowed down with an increase in dimensionality of the lattice to three spatial dimensions (3D). In order to accommodate, without frustration, both the two-sublattice Néel $(\mathbf{k}=(\pi,\pi))$ and the stripe $(\mathbf{k}=(\pi,0))$ orders of the square lattice in a 3D lattice, we require a bipartite lattice which itself is composed of two interpenetrating bipartite lattices, i.e., it is a bi-bipartite lattice. The body centered cubic (BCC) lattice [Fig. 1] has precisely this property; it is a Bravais lattice which is composed of two interpenetrating, identical simple cubic sublattices, and thus serves as a natural analogue of the square lattice in 3D Schmidt et al. (2002). We investigate the classical and $S=1/2$ Heisenberg model on the BCC lattice in the presence of nearest-neighbor $J_{1}$, second-neighbor $J_{2}$ and third-neighbor $J_{3}$ exchange couplings,

[TABLE]

where the $\mathbf{\hat{S}}_{i}$ are the $S=1/2$ Heisenberg spin operators on site $i$. In the classical limit ($S\to\infty$), the $\mathbf{\hat{S}}_{i}$ reduce to three-component vectors. The symbols $\langle i,j\rangle_{1}$, $\langle i,j\rangle_{2}$ and $\langle i,j\rangle_{3}$ denote sums over nearest-neighbor, second-neighbor, and third-neighbor pairs of sites, respectively. The $J_{1}$, $J_{2}$ and $J_{3}$ are allowed to be both FM and AF, and thus we will consider all possible combinations of the signs of the couplings in Eq. (1). Early interest and investigations into the model (with additional four spin interaction terms) stemmed from its relevance to the description of the BCC phase of solid ${}^{3}{\rm He}$ at low-temperatures Utsumi and Izuyama (1977); Okada and Ishikawa (1978); Yosida (1980).

The classical ground states of the BCC $J_{1}$–$J_{2}$ model are Utsumi and Izuyama (1977); Shender (1982): (i) For $J_{2}/|J_{1}|<2/3$, a FM state [Fig. 2(a)] with $\mathbf{k}=(0,0,0)$ for FM $J_{1}$ or a two-sublattice Néel state [Fig. 2(b)] with $\mathbf{k}=(2\pi,0,0)$ for AF $J_{1}$. (ii) For $J_{2}/|J_{1}|>2/3$, a stripe antiferromagnet [Fig. 2(c)] with $\mathbf{k}=(\pi,\pi,\pi)$ is stabilized in both cases, a FM or AF $J_{1}$. This is because the transition point depends only on the coordination number at nearest-neighbor $z_{1}$ ($=8$) and second-neighbor $z_{2}$ ($=6$) distances, with the critical $(J_{2}/|J_{1}|)_{c}=z_{1}/2z_{2}$, hence $(J_{2}/|J_{1}|)_{c}=2/3$ Utsumi and Izuyama (1977); Schmidt et al. (2002). For the corresponding $S=1/2$ BCC $J_{1}$–$J_{2}$ model, all previous studies suggest a single direct phase transition from the FM or Néel state to the stripe ordered state Tahir-Kheli and Jarrett (1964); Müller, Patrick et al. (2015); Schmidt et al. (2002); Oitmaa and Zheng (2004a); Majumdar and Datta (2009); Pantić et al. (2014); Farnell et al. (2016). Thus, in contrast to the square lattice $S=1/2$ $J_{1}$–$J_{2}$ Heisenberg model Schulz and Ziman (1992); H.J. Schulz et al. (1996); Shannon et al. (2006); Iqbal et al. (2016a), there is an absence of an intermediate quantum paramagnetic phase, a manifestation of the weakening of quantum fluctuations in 3D. Note however, that for the square lattice model with FM $J_{1}$, the very existence of quantum paramagnetic phase is not very clear yet Shindou et al. (2011); Richter et al. (2010); Iqbal et al. (2016a). On the BCC lattice, the role of a further neighbor frustrating AF $J_{3}$ coupling in Eq. (1) has not yet been investigated, neither at the classical or semiclassical level nor in the limit of small spin–$S$. At the classical level, our Luttinger-Tisza analysis shows that the inclusion of an AF $J_{3}$ coupling stabilizes a plethora of helimagnetic structures, and a planar AF order Utsumi and Izuyama (1977). In particular, for a model with FM $J_{1}$ we find three incommensurate spiral orders, namely, a 1D spiral with $\mathbf{q}=(q,0,0)$ [Fig. 2(e)], a 2D spiral with $\mathbf{q}=(q,q,0)$ [Fig. 2(f)], and a 3D spiral with $\mathbf{q}=(q,q,q)$ [Fig. 2(g)]. Similarly, in the case of AF $J_{1}$ we find three corresponding incommensurate spiral orders, namely, a 1D spiral with $\mathbf{q}=(2\pi-q,0,0)$ [Fig. 2(h)], a 2D spiral with $\mathbf{q}=(2\pi-q,q,0)$ [Fig. 2(i)], and a 3D spiral with $\mathbf{q}=(2\pi-q,q,q)$ [Fig. 2(j)]. In addition, for both FM and AF $J_{1}$, a planar AF order with $\mathbf{q}=(\pi,\pi,0)$ [Fig. 2(d)] is stabilized at large $J_{2}$ and $J_{3}$. The global classical phase diagram is presented in Fig. 4(a) and Fig. 4(c) together with the pitch vectors of these incommensurate spirals given in Table 1.

For the quantum $S=1/2$ $J_{1}$–$J_{2}$–$J_{3}$ model for both FM and AF $J_{1}$, our PFFRG analysis reveals that the most salient manifestation of quantum fluctuations is the realization of a quantum paramagnetic (PM) phase centered at the tricritical point of the 2D spiral, 3D spiral and planar AF orders [Fig. 4(b) and Fig. 4(d)]. This PM phase has an extended span in parameter space and is stabilized principally at the expense of the 2D and 3D spiral orders, and lesser so at the cost of the stripe and planar AF orders. The phase boundaries and the pitch vectors of the helimagnetic orders are found to be strongly renormalized compared to their classical values, however, no new magnetic orders are found to be stabilized by quantum fluctuations. We estimate the critical magnetic ordering temperature $T_{c}$ of the Néel and FM orders in the $J_{1}$–$J_{2}$ Heisenberg model and compare our findings against Quantum Monte Carlo estimates for nonfrustrated case of nearest-neighbor FM and AF couplings only, and previously obtained high-temperature series expansion estimates in the frustrated regime.

The paper is organized as follows: In Sec. II, we briefly introduce the main methods which are employed in this paper. These are the Luttinger-Tisza approach (Sec. II.1) which is used to determine the classical phase diagram and the PFFRG method (Sec. II.2) which is employed to map out the quantum phase diagram for $S=1/2$. The following Sec. III presents the results of this study, wherein Sec. III.1 (Sec. III.2) discuss the classical (quantum) phase diagrams. Sec. III.3 is devoted to an analysis of critical magnetic ordering temperatures. Finally, we summarize our findings and present an outlook for future studies in Sec. IV. In Appendix A we provide a brief illustration of the coupled-cluster method (CCM) that we use to complement our calculations. In Appendix B we report some details on the quantum Monte Carlo simulations of the finite-temperature behavior in the unfrustrated regime.

II Methods

II.1 Luttinger-Tisza method

The classical limit of a quantum spin model can be obtained by replacing all spin operators on the lattice by unit vectors. For general Heisenberg-type interactions on the BCC lattice, the classical spin Hamiltonian is thus

[TABLE]

where the Heisenberg spin operators $\mathbf{\hat{S}}_{i}$, such as in Eq. (1), are reduced to standard three-component unit vectors $\mathbf{S}_{i}$. Here, $\mathbf{R}_{i}$ is the position of the lattice site $i$. Formally, this limit can be understood as normalization of the spin operators by the total angular momentum $\sqrt{S(S+1)}$, and subsequently taking the limit $S\to\infty$ Millard and Leff (1971); Lieb (1973). The Luttinger-Tisza method Luttinger and Tisza (1946); Luttinger (1951); Kaplan and Menyuk (2007) is an approach to find the approximate groundstate of the classical model in Eq. (2) by replacing the unit length constraint for each classical spin vector with a global constraint $\sum_{i}|\mathbf{S}_{i}^{2}|=S^{2}N$, called a weak constraint. Here, $N$ denotes the total number of spins in the system. This approximation, in principle, allows for local fluctuations in the spin length as only the average of the local moments is fixed. On Bravais lattices, such as the BCC lattice, however, the ground state subject to the weak constraint automatically fixes $|\mathbf{S}_{i}^{2}|=S^{2}\ \forall~{}i$, which renders the Luttinger-Tisza method exact on Bravais lattices.

To understand this and also solve the weakly constrained problem we switch to reciprocal space where the classical Heisenberg Hamiltonian Eq. (2) reads

[TABLE]

Here, we have used the Fourier transform of the spin configuration

[TABLE]

An analogous expression also gives the Fourier transform of the interaction $J(\mathbf{k})$.

On a Bravais lattice, the normalized Fourier modes in real space are planar spin spirals given by

[TABLE]

where $R\in O(3)$ is an arbitrary rotation and reflection matrix allowed by the symmetry of the Heisenberg model. The choice of sign in the first component of $\mathbf{S}_{i}$ reflects the chirality of the spiral and is not fixed within the Heisenberg model alone and instead requires the presence of terms anisotropic in spin space such as dipolar interactions. The classical ground state is obtained by minimizing $J(\mathbf{k})$ with respect to $\mathbf{k}$ within the first Brillouin zone, and the real space spin configuration is then given by the corresponding spin spiral. For the case of the $J_{1}$-$J_{2}$-$J_{3}$ Heisenberg model on the BCC lattice in Eq. (1) we obtain

[TABLE]

where, $k_{x}$, $k_{y}$ and $k_{z}$ are the three components of the wavevector $\mathbf{k}$. It is possible to analytically carry out the minimization of $J(\mathbf{k})$, and the wave vector corresponding to the minima in $\mathbf{k}$-space is termed as the ordering wavevector, and subsequently denoted by $\mathbf{q}$. The ordering wave vectors can be unique or degenerate for a particular ordered state, but are always distinct for two different magnetic orders. Hence, an ordered state can be uniquely specified by its $\mathbf{q}$ vector(s).

Employing this scheme, we obtain all the different classical ground states stabilized in the $J_{1}$–$J_{2}$–$J_{3}$ BCC Heisenberg model. The analytical minimizations are performed only along the high symmetry lines of the first Brillouin zone of the BCC lattice, where all the $\mathbf{q}$’s for the different ground states are found to be located. Away from these lines numerical minimizations are done only for completeness. Mapping out the respective ground states in the $J_{1}$–$J_{2}$–$J_{3}$ parameter space permits us to build the complete analytical phase diagram of the Heisenberg model on the BCC lattice. We also obtain the analytical expressions for the phase boundaries between different magnetically ordered ground states, and the order of the phase transitions.

II.2 Pseudofermion functional renormalization group method

We now briefly introduce the PFFRG method which is used to calculate the quantum ($S=1/2$) phase diagram of the system. The first key step of this approach Reuther and Wölfle (2010); Reuther and Thomale (2011); Reuther et al. (2011a, b, c); Reuther and Thomale (2014) is to re-express the spin Hamiltonian, e.g., Eq. (1), in terms of Abrikosov pseudofermions using $\mathbf{\hat{S}}_{i}=\frac{1}{2}\sum_{\alpha,\beta}\hat{c}^{\dagger}_{i,\alpha}\boldsymbol{\sigma}_{\alpha\beta}\hat{c}_{i,\beta}$ Abrikosov (1965), where $\alpha$, $\beta=\uparrow$ or $\downarrow$, and $\hat{c}^{\dagger}_{i,\alpha}$ ($\hat{c}_{i,\alpha}$) are the pseudofermion creation (annhilation) operators, and $\boldsymbol{\sigma}$ is the Pauli matrix vector. The introduction of pseudofermions leads to an artificial enlargement of the Hilbert space, which, apart from the local physical spin states |$\uparrow\rangle=|1,0\rangle$ and |$\downarrow\rangle=|0,1\rangle$ also contains empty and doubly occupied states, $|0,0\rangle$ and $|1,1\rangle$, respectively, carrying zero spin (here, the notation $|n_{\uparrow},n_{\downarrow}\rangle$ indicates the occupations of the $\uparrow$ and $\downarrow$ modes). The problem of possible spurious contributions from unphysical $S=0$ states in the PFFRG can be cured by adding level repulsion terms $-A\sum_{i}{\bf S}_{i}^{2}$ to the Hamiltonian which, upon choosing $A$ sufficiently large (and positive), energetically separate the Hilbert spaces of the |$\uparrow\rangle$, |$\downarrow\rangle$ and $|0,0\rangle$, $|1,1\rangle$ states, respectively Baez and Reuther (2017). For most practical purposes (such as an application to the models considered here), it is sufficient to set $A=0$ since the single occupancy constraint is already naturally fulfilled in the ground state Baez and Reuther (2017). This is because an unphysical occupation may be viewed as a vacancy in the spin lattice associated with a finite excitation energy.

The resulting fermionic theory is then equipped with a step-like infrared frequency cutoff $\Lambda$ suppressing the bare fermion propagator between $\omega=-\Lambda$ and $\omega=\Lambda$ on the Matsubara axis. This manipulation generates a $\Lambda$-dependence of all $m$-particle vertex functions which may be formulated as an exact but infinite hierarchy of coupled differential equations Metzner et al. (2012); Platt et al. (2013), where the ones for the self energy $\Sigma$ and for the two-particle vertex, $\Gamma_{2}$, are illustrated in Fig. 3(a) and Fig. 3(b). To be amenable to numerical solutions, this hierarchy of equations needs to be truncated, which, for the results presented below, amounts to approximating the three-particle vertex, $\Gamma_{3}$, using two different schemes:

(i) One-loop plus Katanin scheme: This approach has been widely used and has proven suitable to capture the right balance between ordering tendencies and quantum fluctuations. It will be applied in most calculations presented below. Within this scheme, the contracted three-particle vertex in Fig. 3(b) takes the form of Fig. 3(c) which corresponds to a self-energy insertion in the interaction channels of the two-particle vertex flow. The crucial benefit of this approximation as compared to neglecting $\Gamma_{3}$ completely is that it guarantees the full feedback of the self-energy into the two-particle vertex flow, hence, leading to a fully self-consistent RG scheme. From a different perspective, it can be shown that the one-loop plus Katanin scheme exactly sums up all diagrammatic contributions separately in the large $S$ limit Baez and Reuther (2017) and in the large $N$ limit (where in the latter case the spins’ symmetry group is promoted from SU(2) to SU($N$)) Buessen et al. (2018a). This ensures that magnetically ordered phases (as typically encountered for $S\rightarrow\infty$) as well as disordered phases (obtained for $N\rightarrow\infty$) may both be faithfully described. The three-particle terms which are neglected within the one-loop plus Katanin scheme can be shown to be subleading in both $1/S$ and $1/N$.

(ii) Two-loop scheme: This approach adds further corrections to the three-particle term such as those shown in Fig. 3(d). In this diagram, different two-particle interaction channels are inserted into each other resulting in effective two-loop contributions. It should also be noted that in similarity to the Katanin scheme in (i), self-consistency again requires the full feedback of self-energy into such nested diagrams which even generates certain three-loop contributions (see Rück and Reuther (2018) for details). All these corrections ensure that the aforementioned subleading terms in $1/S$ and $1/N$ are better approximated which allows for a more accurate investigation of quantum critical parameter regions where the detailed interplay between magnetic ordering and quantum fluctuations becomes crucial. While this may generally lead to shifted phase boundaries compared to the scheme in (i) it has been shown in Ref. Rück and Reuther (2018) that such shifts turn out to be rather small. Another benefit of the two-loop scheme is that it determines Néel/Curie temperatures more accurately Rück and Reuther (2018) which is also the context in which it will be applied below.

Up to fermionic contractions, the two-particle vertex [either calculated via (i) or (ii)] is the diagrammatic representation of the static (i.e., imaginary time-integrated) spin correlator given by

[TABLE]

where $\hat{S}_{i}^{\mu}(\tau)=e^{\tau\hat{\mathcal{H}}}\hat{S}_{i}^{\mu}e^{-\tau\hat{\mathcal{H}}}$. Since the Heisenberg model is spin-rotation invariant all diagonal components $C_{ij}^{\mu\mu}$ of the spin correlator are identical. Without loss of generality we have chosen the $zz$-component here. Within PFFRG, the thermodynamic limit is approached by calculating the correlators $C_{ij}$ only up to a maximal distance between sites $i$ and $j$. Fourier-transforming $C_{ij}$ into momentum space, we then obtain the static susceptibility $\chi^{\Lambda}({\bf k})$ as a function of $\Lambda$ which represents the central physical outcome of this approach. While the Fourier-transform generally allows to access a continuous set of wave vectors k within the Brillouin zone, the restriction to a finite set of correlators $C_{ij}$ limits the number of harmonics in the Fourier-sums and, therefore, smoothens sudden changes in the susceptibility. In the present study of the $J_{1}$–$J_{2}$–$J_{3}$ model on the BCC lattice, we have set the maximal length of spin correlators to be equal to 10 nearest-neighbor lattice spacings, which incorporates a total of $2331$ correlated sites, producing well converged results with a proper k-space resolution. We, furthermore, approximate the frequency dependences of the vertex functions by discrete grids containing $64$ or $100$ points for each frequency variable. The number of coupled differential equations for the above given system size and 64 frequencies is (i) $611\,057\,728$—without using any point group symmetry and (ii) $1\,537\,568$—upon exploiting the complete $O_{h}$ point group Landau and Lifshitz (1977) symmetries, while for a calculation with $100$ frequencies, we have (i) $2\,331\,000\,100$—without symmetries and (ii) $5\,801\,300$—with symmetries. When a system spontaneously develops magnetic order, the susceptibility $\chi^{\Lambda}({\bf k})$ shows a sharp increase at the corresponding wave vector q upon decreasing $\Lambda$ and eventually the flow becomes unstable. The $\Lambda$ values at which this breakdown takes place can be associated with the critical ordering temperature via the relation $T_{c}=(\frac{\pi}{2})\Lambda_{c}$ for $S=1/2$ Iqbal et al. (2016b, 2019). In contrast, a smooth flow of the susceptibility down to $\Lambda\rightarrow 0$ indicates a magnetically disordered state. After the initial applications of the PFFRG in two dimensions, it has subsequently been applied with much success to three-dimensional systems Iqbal et al. (2016b); Balz et al. (2016); Iqbal et al. (2019); Buessen and Trebst (2016); Iqbal et al. (2017); Chillal et al. (2017); Iqbal et al. (2018); Buessen et al. (2018b). For further details about the PFFRG procedure, its subsequent refinements, and expansions to handle a larger class of magnetic Hamiltonians we refer the reader to Refs. Reuther and Wölfle (2010); Reuther and Thomale (2011); Reuther et al. (2011b); Iqbal et al. (2016b, c); Hering and Reuther (2017); Keleş and Zhao (2018); Buessen et al. (2018c); Roscher et al. (2018); Iqbal et al. (2015).

III Results

III.1 Classical Phase diagram

\floatsetup

[table]capposition=bottom

We begin by presenting the classical ground state phase diagram, i.e., at $T=0$ and a detailed analysis of the magnetic orders in the classical $J_{1}$–$J_{2}$–$J_{3}$ Heisenberg model on a BCC lattice as obtained from the Luttinger-Tisza method [Sec. II.1].

III.1.1 Ferromagnetic $J_{1}$

The classical phase diagram in the $J_{1}$–$J_{2}$–$J_{3}$ parameter space with FM $J_{1}$ is shown in Fig. 4(a). It is host to six different types of magnetic orders; three incommensurate coplanar spiral structures and three collinear orders. Starting with both FM $J_{1}$ and $J_{2}$ we trivially find a FM ground state [Fig. 2(a)]. The inclusion of an AF $J_{3}$ coupling above a critical value $J_{3}^{c}=\frac{1}{4}(|J_{1}|+|J_{2}|)$ destabilizes the FM state Kaplan (1959), via a 2nd order phase transition, into an incommensurate 1D spiral [Fig. 2(e)] with a pitch vector $\mathbf{q}=(q,0,0)$ with $q$ given in Table 1. This pitch vector is 6–fold degenerate within the first Brillouin zone. It is important to emphasize that the spiral state is only governed by one of these symmetry equivalent pitch vectors, i.e., superpositions are not possible as they would violate the classical length constraint, and hence the degeneracy remains discrete. This 1D spiral structure is stabilized purely by a FM $J_{2}$ interaction. Indeed, along the line $J_{2}=0$, there is a 1st order phase transition to a 2D incommensurate spiral [Fig. 2(f)] with a pitch vector $\mathbf{q}=(q,q,0)$ with $q$ given in Table 1. Similar to the 1D spiral, the ground state in this phase is determined only by one of the $12$ symmetry equivalent $(q,q,0)$-type pitch vectors. Upon increasing $J_{2}$, we observe that the value of $q$ continuously evolves towards $\pi$. At the line $J_{2}/|J_{1}|=1$ and above a critical $J_{3}^{c}/|J_{1}|\approx 0.29$ [see Table 2 for an analytical expression of the phase boundaries], there is a 2nd order phase transition to a planar AF order [Fig. 2(d)] with $\mathbf{q}=(\pi,\pi,0)$ Utsumi and Izuyama (1977). In contrast to the incommensurate orders discussed above, the pitch vector of the planar AF is half of a reciprocal lattice vector, i.e., $2\mathbf{q}\equiv\mathbf{0}$. As pointed out by Villain Villain, J. (1977), this characteristic allows the ground state to be composed of all six symmetry equivalent pitch vectors, namely, $(\pi,\pi,0)$, $(\pi,0,\pi)$, $(0,\pi,\pi)$, $(\pi,-\pi,0)$, $(-\pi,0,\pi)$, and $(0,-\pi,\pi)$. All six pitch vectors satisfy the property $\sin(\mathbf{q}\cdot\mathbf{R}_{i})=0$ at every lattice site. Therefore, the general ground state can be written as

[TABLE]

where $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$, $\mathbf{d}$, $\mathbf{e}$, $\mathbf{f}$ are arbitrary vectors constrained by Eq. (9) which normalizes the spin length at each site. The arbitrary vectors defining the spin configuration have in total $18$ continuous degrees of freedom, and are subject to $8$ constraints. Accounting for the global spin rotation invariance (2 degrees of freedom) of the Heisenberg model, the continuous ground state manifold of the planar AF order is $8$-dimensional. This is in contrast to the other magnetic orders in the $J_{1}$–$J_{2}$–$J_{3}$ parameter space which feature only a $n$-fold discrete degeneracy [see Table 1]. A similar enhancement for the available degrees of freedom is also found for the stripe AF phase on the square lattice Danu et al. (2016) and other cubic lattice systems Ignatenko and Irkhin (2016).

Upon decreasing the value of $J_{3}$, we find that the interplay between AF $J_{2}$ and $J_{3}$ couplings, leads to the appearance of an incommensurate 3D spiral [Fig. 2(g)] in a sliver of parameter space. This state also continuously evolves from the FM state via a 2nd order phase transition, however, its transition into the $(q,q,0)$ and $(\pi,\pi,0)$ states is of 1st order. Its pitch vector $\mathbf{q}=(q,q,q)$ is 8–fold degenerate, but again, only one of them is present in any given ground state. Lowering the AF $J_{3}$ coupling even further, and for $J_{2}/|J_{1}|>2/3$, the collinear stripe order [Fig. 2(c)] with a wave vector $\mathbf{q}=(\pi,\pi,\pi)$ is stabilized. This state is composed of two interpenetrating simple cubic lattices which are Néel ordered. For any spin on a given sublattice, all its nearest-neighbor spins residing on the other sublattice add up to zero. Hence, the energy of the $(\pi,\pi,\pi)$ state is independent of the relative orientation of the two sublattices Shender (1982), which is thus not determined within the $J_{1}$–$J_{2}$–$J_{3}$ Heisenberg model. The pitch vector of this state resides at the corners of the first Brillouin zone, and is therefore unique. Since the third neighbor spins in this state are FM ordered, changing $J_{3}$ to FM only enhances its stability, and thus this state occupies the entire parameter space for $J_{2}/|J_{1}|>2/3$ and FM $J_{3}$. The phase boundary between the FM and the stripe collinear order is determined solely by the coordination number at nearest-neighbor and second-neighbor distances, and is given by $J_{2}/|J_{1}|=z_{1}/2z_{2}=2/3$.

III.1.2 Antiferromagnetic $J_{1}$

A change in the $J_{1}$ coupling from FM to AF is found not to alter the phase boundaries and the order of the phase transitions in the $J_{2}$–$J_{3}$ parameter space, as observed in the corresponding phase diagram [Fig. 4(c)] and Table 2. This feature is most easily understood by viewing the BCC lattice as being composed of two interpenetrating simple cubic sublattices with one being positioned at the body center of the other. As $J_{2}$ and $J_{3}$ couple sites only within the sublattices, it is only the $J_{1}$ coupling which connects the two simple cubic lattices. Hence, a sign reversal of $J_{1}$ can be undone by flipping the spins on one of the sublattices without affecting the $J_{2}$ and $J_{3}$ couplings. Therefore, the phase boundaries remain unchanged, however, in reciprocal space this flipping amounts to a shift of the wave vector $\mathbf{q}\to(2\pi,2\pi,2\pi)-\mathbf{q}$. After folding back this wave vector into the first Brillouin zone this amounts to a shift of any one of the wave vector components $q\to 2\pi-q$. Consequently, the FM state is replaced by a Néel AF [Fig. 2(b)] with wavevector $\mathbf{q}=(2\pi,0,0)$ 111In literature, the Néel AF state is also specified by the wave vector $\mathbf{q}=(2\pi,2\pi,2\pi)$. This wave vector does not reside in the first Brillouin zone, however, it is equivalent to $\mathbf{q}=(2\pi,0,0)$. Similarly, the incommensurate 1D spiral is now characterized by the pitch vector $\mathbf{q}=(2\pi-q,0,0)$ [Fig. 2(h)] with $q$ and the degeneracy of the ground state being the same as in the FM case [see Table 1]. Along the same lines, the pitch vector of the incommensurate 2D spiral in the AF $J_{1}$ case is given by $\mathbf{q}=(2\pi-q,q,0)$, which is now $24$–fold degenerate in contrast to the $12$–fold degeneracy present in the ferromagnetic $J_{1}$ case. Nonetheless, the ground state is still composed of only one of these pitch vectors. This state also evolves continuously to the planar AF upon increasing $J_{2}$, and the planar AF remains unchanged compared to the FM $J_{1}$ case. The incommensurate 3D spiral now has a $24$–fold degenerate pitch vector $\mathbf{q}=(2\pi-q,q,q)$. Finally, the stripe AF $\mathbf{q}=(\pi,\pi,\pi)$-phase remains unchanged.

The complete analytical expressions for the pitch vectors of the spiral orders together with their degeneracies are to be found in Table 1. Also, in Table 2 we provide the analytical expressions of all phase boundaries together with the order of the phase transitions.

III.2 Quantum Phase diagram

We now investigate the effects of quantum fluctuations on the classical phase diagram employing one-loop PFFRG. As found in Ref. Rück and Reuther (2018), the one-loop formulation is mostly sufficient to correctly determine the phase boundaries for spin $S=1/2$. The more numerically-intensive two-loop scheme will only be applied for a select set of coupling parameters in order to calculate critical magnetic ordering temperatures more accurately, the results of which are presented in Sec. III.3. As described in Sec. II.2, at each point in parameter space we track the evolution of the susceptibility $\chi^{\Lambda}({\bf k})$ as a function of $\Lambda$ for all $\mathbf{k}$ in the first Brillouin zone. The $\mathbf{k}$-vector which yields the dominant susceptibility at the point of breakdown of the RG flow then determines the nature of the magnetically ordered ground state. On the other hand, the absence of a breakdown in the limit $\Lambda\to 0$ signals the absence of long-range dipolar magnetic order, and points to a paramagnetic ground state. We find that all the magnetic orders in the classical phase diagram are to be found in the quantum phase diagram [Fig. 4(b) and Fig. 4(d)], and that no new types of long-range dipolar magnetic orders are stabilized by quantum fluctuations. In Fig. 5, we show the representative RG flows for all the quantum phases with their momentum resolved susceptibility profiles shown in Fig. 6 and Fig. 7.

For $S=1/2$, the most salient effect of quantum fluctuations is the appearance of a PM phase over an extended region in the $J_{2}$–$J_{3}$ parameter space for FM as well as AF $J_{1}$ [see Fig. 4(b) and Fig. 4(d)]. The PM phase extends over a larger region in the case of FM $J_{1}$ as compared to AF $J_{1}$. It engulfs a substantial portion of parameter space occupied classically by the 2D and 3D incommensurate spiral orders, and to a lesser degree cuts into the classical domain of the $\mathbf{q}=(\pi,\pi,0)$ planar AF state. For FM $J_{1}$, we find that a tiny region of the classical FM phase is destabilized into a PM phase by quantum fluctuations. On the other hand, for AF $J_{1}$, the Néel order does not succumb at all to quantum fluctuations. Interestingly, we find that for FM as well as AF $J_{1}$ quantum fluctuations do not destabilize the $\mathbf{q}=(\pi,\pi,\pi)$ stripe AF order into a PM phase, and consequently, the PM phase does not occupy any portion of the parameter space which classically hosts the stripe AF. In Fig. 8, we present the momentum resolved susceptibilities for the PM phase. For completeness we have also considered larger spin magnitudes $S>1/2$, using a modification of the one-loop PFFRG as described in Ref. Baez and Reuther (2017). Interestingly, we find that already at $S=1$ the PM phase disappears entirely owing to the weakening of quantum fluctuations.

In the following, we investigate in more detail, the effects of quantum fluctuations for $S=1/2$ on the different magnetically ordered phases. Although, our principal method to investigate the quantum phase diagram is the PFFRG, for comparison we also add here results using the exact diagonalization (ED) and the coupled-cluster method (CCM). Both methods have been used previously to study the $S=1/2$ BCC $J_{1}$–$J_{2}$ model with AF $J_{1}$ Schmidt et al. (2002); Farnell et al. (2016), but so far no results for FM $J_{1}$ are available. The ED is a well-established method, see, e.g., H.J. Schulz et al. (1996); Schmidt et al. (2002); Richter et al. (2010). Here, we use J. Schulenburg’s spinpack Schulenburg (2018); Richter and Schulenburg (2010). The CCM is a universal many-body method  Bishop (1991); Zeng et al. (1998); Bishop (1998); Farnell and Bishop (2004) that has been successfully applied on frustrated quantum spin systems, see, e.g. Darradi et al. (2005, 2008); Farnell et al. (2009); Götze et al. (2011); Li et al. (2015); Bishop and Li (2017); Farnell et al. (2018). We will give a brief illustration of the CCM in Appendix A. Moreover, we mention that ED and CCM calculations used here follow closely Ref. Schmidt et al. (2002) and Ref. Farnell et al. (2016), respectively. The main ED and CCM results are summarized in Fig. 9.

\floatsetup

[table]capposition=bottom

Firstly, the only impact of quantum fluctuations on the stripe AF state is to shift its phase boundaries compared to the classical ones. In particular, we find that for large enough $J_{2}$ (i) the stripe AF replaces the classical incommensurate 3D spirals, whose existence is thus reduced to a tiny sliver in the $J_{2}$–$J_{3}$ plane and (ii) the stripe AF cuts into the classical domain of the planar AF with which it now shares a phase boundary hitherto absent in the classical phase diagram. This is similar to the findings on the square lattice wherein quantum fluctuations are found to favor the $(\pi,\pi)$ state over the $(\pi,0)$ state Danu et al. (2016). In contrast, quantum fluctuations act differently on the other phase boundary of the stripe AF with the FM or Néel orders depending on whether $J_{1}$ is FM or AF, respectively. For FM $J_{1}$, we find that the phase boundary shifts to a smaller value of $J_{2}$, whereas for AF $J_{1}$ the phase boundary is shifted to a larger value as also observed on the simple cubic lattice Iqbal et al. (2016b). In particular, there is no intermediate PM phase in between the FM/Néel and stripe AF orders in the $J_{1}$–$J_{2}$ model in agreement with previous studies Tahir-Kheli and Jarrett (1964); Müller, Patrick et al. (2015); Farnell et al. (2016); Schmidt et al. (2002); Majumdar and Datta (2009); Pantić et al. (2014); Oitmaa and Zheng (2004a). This is in contrast to the findings on the square lattice, and can be attributed to the diminished quantum fluctuations in 3D. In Table 3, we provide numerical estimates of the phase boundary for $J_{3}=0$, i.e., along the $J_{1}$–$J_{2}$ line obtained by PFFRG and other numerical approaches. The observation that for FM $J_{1}$ the stripe AF order extends at the expense of FM order follows from the fact that quantum fluctuations, in general, do not alter the FM state (including its ground state energy) as its an eigenstate of the Heisenberg exchange Hamiltonian and thus free of macroscopic zero-point vibrations Nagaev (1984); Kaganov and Chubukov (1987), however, act on AF orders, e.g., by lowering their ground state energies. Hence, compared to the classical case, phase boundaries between FM and AF orders are typically shifted towards the FM side. One further observes that the shift in the phase boundary from the classical value of $J_{2}^{c}/|J_{1}|=2/3$ is stronger for FM $J_{1}$ compared to AF $J_{1}$. We mention that our PFFRG findings for the BCC $J_{1}$–$J_{2}$ model are in excellent agreement with ED and CCM results, cf. Table 3, where we provide a comparison of the critical values $J_{2}^{c}/|J_{1}|$ of the transition between the FM/Néel and the stripe orders. Moreover, the ED and CCM results for the spin-spin correlation functions and the order parameter shown in Fig. 9 clearly demonstrate an absence of an intermediate quantum paramagnetic phase, cf. the PFFRG phase diagram in Fig. 4(b).

A plot of the susceptibility along a path in reciprocal space for different $J_{2}/|J_{1}|$ within the FM and Néel ordered phases is presented in Fig. 10. The maxima of the susceptibility at the magnetic wave vectors of the respective orders are seen to be clearly resolved. The frustrating effect of a $J_{2}$ coupling on the FM/Néel orders leads to a reduction in the dominant susceptibility peaks and to the development of a peak at the incipient stripe AF order at $\mathbf{q}=(\pi,\pi,\pi)$, similar to the findings by high temperature series expansion Richter et al. (2015).

Another phase boundary which is significantly shifted by quantum fluctuations is the one between the 2D spiral and planar AF orders [see Fig. 4]. For FM as well as AF $J_{1}$, the classical phase boundary at $J_{2}/|J_{1}|=1$ is shifted to a smaller value, implying that quantum effects enhance the stability of the planar AF state. In particular, for FM $J_{1}$ the phase boundary shifts to $J_{2}/|J_{1}|\approx 0.82$, while for AF $J_{1}$ the shift is comparatively weaker, and the phase boundary is found to be located at $J_{2}/|J_{1}|\approx 0.93$. This strong effect of quantum fluctuations can be explained by the fact that the $\mathbf{q}=(\pi,\pi,0)$ planar AF has a continuous $8$–dimensional classical ground state degeneracy whereas the 2D spiral only features an $8$–fold discrete degeneracy. Consequently, quantum fluctuations play a more prominent role on top of the classical planar AF state in comparison to the 2D spiral.

Similar to the shifts observed in the phase boundary between the FM/Néel and stripe AF order, we find that the boundaries between FM/Néel orders to the 1D and 3D spiral phases behave differently depending on whether $J_{1}$ is FM or AF. For FM $J_{1}$, the 1D and 3D spirals enhance their domain of stability beyond the classically allowed region of their existence, and thus the domain of the FM phase shrinks compared to the classical one. In contrast for AF $J_{1}$, the Néel phase extends into these spiral orders but to a lesser extent. In total, the domain of existence of the FM order is significantly reduced by quantum fluctuations, whereas for the Néel phase it is enhanced, compared to the classical phase diagram. Finally, two phase boundaries remain unaffected by quantum fluctuations, namely, the 2D to 3D spirals, and the one between 1D and 2D spirals.

Now, we discuss the impact of quantum fluctuations within the incommensurate spiral orders which primarily amounts to a shift of the value of the pitch vectors. Indeed, our PFFRG analysis shows that for FM $J_{1}$ the shift is such that the 1D spiral pitch vector is shifted towards that of the Néel state. At a fixed $J_{2}$, this effect is stronger for small $J_{3}$ and appears to decrease with increasing $J_{3}$ [Fig. 11(a)]. At a fixed $J_{3}$, this shift of $\mathbf{q}$ towards the Néel pitch vector increases with increasing FM $J_{2}$ [Fig. 11(c)]. Note that quantum fluctuations seem to act counterintuitively here since with increasing strength of the FM $J_{2}$ coupling one would expect the FM state to become increasingly favorable. Interestingly, in the case of AF $J_{1}$ the 1D spiral pitch vectors shift in the opposite direction, i.e., they migrate towards the FM state [see Fig. 11(a)], except for large FM $J_{2}$ where essentially no shift is observed [see Fig. 11(c)]. For the 2D spirals, we observe that when $J_{1}$ is FM, the shift of the pitch vector is towards that of the planar AF order, and the magnitude of the shift decreases with increasing $J_{3}$ [see Fig. 11(b)]. At fixed $J_{3}$, the magnitude of the shift remains essentially constant with varying $J_{2}$ [see Fig. 11(c)]. This finding is consistent with the fact that the phase boundary of the planar AF state shifts to a smaller value of $J_{2}$. However, for AF $J_{1}$ we find that the 2D spiral pitch vector is shifted towards that of the Néel state, and upon varying $J_{2}$ and $J_{3}$ the magnitude of the shift remains essentially constant [see Fig. 11(b) and Fig. 11(c)]. We do not show the changes of the $\bf q$ vector for the ${(q,q,q)}$ spiral since it covers only a tiny sliver in the quantum phase diagram.

III.3 Néel and Curie temperatures

\floatsetup

[table]capposition=bottom

\floatsetup

[table]capposition=top

The magnetic ordering temperature, i.e., the Curie temperature ($T_{C}$) for a FM and the Néel temperature ($T_{N}$) for an AF ordered state is one of the fundamental thermodynamic quantities which serves as a measure of the degree of frustration. The (numerically exact) quantum Monte Carlo (QMC) method can be employed to calculate this $T_{C}$ and $T_{N}$ in the non-frustrated region of parameter space. However, in the frustrated regime, one must resort to approximate numerical approaches to obtain estimates of $T_{C}$ and $T_{N}$. Here, we employ one- and two-loop PFFRG to estimate the ordering temperatures for non-frustrated and frustrated coupling parameters of the $J_{1}$–$J_{2}$ BCC Heisenberg model. As observed in Ref. Rück and Reuther (2018), the one-loop PFFRG is less converged when estimating critical ordering temperatures as compared to determining phase boundaries. We, therefore, carry out our calculations in both one-loop and two-loop formulation. Furthermore, we compare these estimates to those obtained by high temperature expansion (HTE), and Green’s function methods in previous studies which also serves and a benchmark test for the performance of the PFFRG. Additionally, we carry out QMC calculations for the nearest-neighbor FM and AF couplings only, and obtain estimates of $T_{C}$ and $T_{N}$ by a finite-size scaling analysis of the renormalization-group invariant quantities Binder ratio and the ratio $\xi/L$ of the second-moment correlation length $\xi$ over the lattice size $L$; see Appendix A of Ref. Parisen Toldin et al. (2015) for a discussion on the definition of $\xi$. More details on the QMC simulations and the analysis are reported in Appendix B.

We start discussing the $S=1/2$ nearest neighbor only model with FM and AF $J_{1}$ interaction. Both systems are unfrustrated such that they are amenable to a QMC calculation. Interestingly, our QMC results show that the Néel and the Curie temperatures are unequal [see Table 4], with $T_{N}$ being greater than $T_{C}$ by about $9\%$, in agreement with the findings from HTE and Green’s function methods Rushbrooke and Wood (1963); Rushbrooke et al. (1964); Oitmaa and Zheng (2004a); Juhász Junger et al. (2009); Richter et al. (2015) (see Table 4). Indeed, it is known to be a general feature of a finite spin-$S$ Heisenberg models on bipartite lattices with nonfrustrating interactions that the Néel and Curie temperatures are unequal with $T_{N}>T_{C}$ Oitmaa and Zheng (2004b). This difference between $T_{C}$ and $T_{N}$ of $\sim 9\%$ is also reflected in our two-loop PFFRG data but is slightly overestimated on the one-loop level where the difference is found to be $\sim 12\%$. Concerning absolute values of the ordering temperatures, our one-loop and two-loop results both slightly overestimate $T_{N}$ and $T_{C}$. By extending the PFFRG from one-loop to two-loop the accuracy of the results becomes significantly better, particularly, the errors of the one-loop critical ordering temperatures are approximately halved in the two-loop results. One may therefore expect that even higher loop orders might give very accurate estimates. We leave such an analysis for future studies. Finally, in Table 5, we compare our results against those for the simple cubic lattice, and also compare the $S=1/2$ ordering temperatures against the ones for the classical model to obtain the reduction due to quantum fluctuations. As expected, the ordering temperatures for the BCC lattice are larger compared to those for the SC lattice due to its higher coordination number. For the same reason, the reduction in the critical temperature for $S=1/2$ with respect to the classical value is lesser for the BCC lattice in comparison to the SC lattice.

In the presence of a frustrating $J_{2}$ interaction there is a significant reduction in both $T_{C}$ and $T_{N}$, which are found to decrease monotonically with increasing $J_{2}$ [see Fig. 12], and on approaching the transition point the ordering temperatures have a sharp drop. Our PFFRG data shows that the inequality $T_{N}>T_{C}$ remains valid up till the transition point into the stripe AF order in agreement with HTE data Oitmaa and Zheng (2004a, b); Müller, Patrick et al. (2015); Richter et al. (2015), but in contrast to results from Green’s function approach Müller, Patrick et al. (2015); Richter et al. (2015) [see Table 4]. Again, we see that the ordering temperatures from PFFRG are slightly larger than those obtained by HTE where two-loop PFFRG mostly gives better estimates. Overall, these results imply that PFFRG (particularly the two-loop formulation) correctly captures the relative behavior of ordering temperatures. The absolute values, however, might still be subject to errors of a few percent which are possibly reduced within higher-loop schemes.

IV Summary and Outlook

We have shown that frustrating the ferromagnetic and Néel antiferromagnetic orders of Heisenberg spins on a three-dimensional bipartite body-centered-cubic lattice by competing interactions up to third neighbors leads to the appearance of a rich variety of helimagnetic and collinear spin structures at the classical level. In the extreme quantum limit of $S=1/2$, our PFFRG analysis shows that the most salient feature of quantum fluctuations is the realization of an extended region of parameter space displaying quantum paramagnetic behavior. The classical phase boundaries are also found to be strongly renormalized by quantum effects, and helimagnetic pitch vectors undergo significant shifts. In total, we find that quantum effects are stronger in the case of a ferromagnetic nearest-neighbor coupling compared to an antiferromagnetic one. We have also estimated the Curie and Néel temperatures from PFFRG, and compared our results to those from quantum Monte Carlo for unfrustrated case of nearest-neighbor only antiferromagnetic and ferromagnetic models, and with available high temperature expansion data in the frustrated regime of the $J_{1}$–$J_{2}$ model. We obtain good agreement with Quantum Monte Carlo for the pure nearest-neighbor Heisenberg ferromagnet and Néel antiferromagnet and reproduce qualitative trends for frustrating couplings. However, we observe that in general the PFFRG overestimates the ordering temperatures which is partially cured by employing two-loop PFFRG.

As a future study, it will be interesting to investigate the finite-temperature classical phase diagram of the $J_{1}$–$J_{2}$–$J_{3}$ model, including its critical properties and nature of phase transitions, which has traditionally largely focussed on the Ising model Banavar et al. (1979); Velgakis and Ferer (1983); Azaria et al. (1989); Murtazaev et al. (2015, 2017), however, recent attempts have been made at the Heisenberg model for a given parameter value Murtazaev et al. (2018); Ramazanov and Murtazaev (2017). The role of disorder in determining the stability of the realized phases is another important issue worth investigating. It has been pointed out in Ref. Attig and Trebst (2017) that if one restricts the second nearest-neighbor coupling $J_{2}^{*}$ to be defined by bond-distance instead of geometrical distance (as in the current paper), then the classical Heisenberg $J_{1}$–$J_{2}^{*}$ antiferromagnet on the BCC lattice hosts spin spiral surfaces analogous to the $J_{1}$–$J_{2}$ model on the diamond lattice Bergman et al. (2007). It will be interesting to investigate the selection effects on the spiral surface due to quantum fluctuations as a function of the frustration ratio $J_{2}^{*}/J_{1}$ and spin-$S$, and in particular, examine the possibility of realizing a spiral spin liquid. Our finding of extended domains characterized by an absence of long-range dipolar magnetic order in the $S=1/2$ model lays the avenue for future numerical investigations aiming to identify the nature of the nonmagnetic phase which could potentially be host to a plethora of exotic nonmagnetic phases such as quantum spin liquids, valence-bond-crystals, and lattice-nematics or feature quadrupolar ordered phases, i.e., spin-nematic orders Andreev and Grishchuk (1984). Indeed, in the $S=1/2$ $J_{1}$–$J_{2}$–$J_{3}$ square lattice Heisenberg model, these orders were found to be stabilized Iqbal et al. (2016a). The question of the microscopic identification of the nature of the nonmagnetic phase can be addressed within the PFFRG framework itself by combining it with a self-consistent Fock-like mean-field scheme to calculate low-energy effective theories for emergent spinon excitations in $S=1/2$ systems as has been recently achieved on the square and kagome lattices Hering et al. (2018). Within this scheme, the effective spin interactions obtained from PFFRG, i.e., the two-particle vertices, act as an input for the Fock equation yielding a self-consistent approach to calculate the spinon band structures beyond a mean field treatment. However, the precise forms of such free spinon Ansätze are given by a projective symmetry group classification Wen (2002), and it will be useful to carry out a classification of the symmetry allowed mean-field quantum spin liquid and nematic states on the BCC lattice. These Ansätze would also then serve as the basis for Gutzwiller projected variational wave-function studies employing Monte Carlo methods Iqbal et al. (2011); Hu et al. (2013); Iqbal et al. (2013, 2016b).

Acknowledgments. T.M., P.G., and Y.I. thank R. Ganesh and Arnab Sen for helpful discussions. J. Richter thanks O. Götze for providing the CCM scriptfiles for the AF BCC model. We gratefully acknowledge the Gauss Centre for Supercomputing e.V. for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ). The work in Würzburg was supported by the DFG through DFG-SFB 1170 tocotronics (project B04) and the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter DFG-EXC 2147/1 ct.qmat (project-id 39085490). F.P.T. thanks the German Research Foundation (DFG) through Grant No. AS120/13-1 of the FOR 1807. Y.I. acknowledges the kind hospitality of the Helmholtz-Zentrum für Materialien und Energie, Berlin, Germany for the period May–July 2018, where part of this work was accomplished. P.G., Y.I., and T.M. acknowledge the interactions and kind hospitality at the International Centre for Theoretical Sciences (ICTS), Bengaluru, India from 29th November till 7 December, 2018 during “The $2$nd Asia Pacific Workshop on Quantum Magnetism”.

Appendix A Brief illustration of the coupled-cluster method (CCM)

We illustrate here only some features of the CCM relevant for the results shown in Fig. 9. At that we follow the lines given in Ref. Farnell et al. (2016), where the CCM was applied to the $J_{1}$–$J_{2}$ BCC model with AF $J_{1}$. For more general information on the methodology of the CCM, see, e.g., Refs.  Bishop (1991); Zeng et al. (1998); Bishop (1998); Farnell and Bishop (2004). We first mention that the CCM yields results directly in the thermodynamic limit $N\to\infty$. First we choose a normalized reference or model state $|\Phi\rangle$ that is here the classical $(\pi,\pi,\pi)$-state, see Fig. 2. Second we perform a rotation of the local axes of each of the spins such that all spins in the model state align along the negative $z$ axis. In this new set of local spin coordinates we define a complete set of mutually commuting multispin creation operators $C_{I}^{+}\equiv(C^{-}_{I})^{\dagger}$ related to this model state: $|{\Phi}\rangle=|\downarrow\downarrow\downarrow\cdots\rangle;\mbox{ }C_{I}^{+}=\hat{S}_{n}^{+}\,,\,\hat{S}_{n}^{+}\hat{S}_{m}^{+}\,,\,\hat{S}_{n}^{+}\hat{S}_{m}^{+}{s}_{k}^{+}\,,\,\ldots\;,$ $\hat{S}^{+}_{n}\equiv\hat{S}^{x}_{n}+i\hat{S}^{y}_{n}$, where the spin operators $\hat{S}^{x}_{n}$ and $\hat{S}^{y}_{n}$ are defined in the local rotated coordinate frames, and the indices $n,m,k,\ldots$ denote arbitrary lattice sites. The CCM parameterizations of the ket and bra GS eigenvectors $|\Psi\rangle$ and $\langle\tilde{\Psi}|$ of the spin system are given by $|\Psi\rangle=e^{S}|\Phi\rangle\;,\mbox{ }S=\sum_{I\neq 0}a_{I}C_{I}^{+}\;;\;$ $\langle\tilde{\Psi}|=\langle\Phi|\tilde{S}e^{-S}\;,\mbox{ }\tilde{S}=1+\sum_{I\neq 0}\tilde{a}_{I}C_{I}^{-}$. The coefficients $a_{I}$ and $\tilde{a}_{I}$ contain the CCM correlation operators, $S$ and $\tilde{S}$. They are determined by the ket-state and bra-state equations $\langle\Phi|C_{I}^{-}e^{-S}He^{S}|\Phi\rangle=0\;;\;\langle\Phi|{\tilde{S}}e^{-S}[H,C_{I}^{+}]e^{S}|\Phi\rangle=0\;;\;\forall I\neq 0.$ Each of these equations labeled by a configuration index $I$, corresponds to a certain configuration of lattice sites $n,m,k,\dots\;$. Using the Schrödinger equation, $H|\Psi\rangle=E_{0}|\Psi\rangle$, we can write the ground-state energy per site as $E_{0}=\langle\Phi|e^{-S}He^{S}|\Phi\rangle$. The magnetic order parameter (sublattice magnetization) is given by $M=-\frac{1}{N}\sum_{i=1}^{N}\langle\tilde{\Psi}|\hat{S}_{i}^{z}|\Psi\rangle$, where $\hat{S}_{i}^{z}$ is expressed in the transformed coordinate system, and $N(\rightarrow\infty)$ is the number of lattice sites. In order to truncate the expansions of $S$ and $\tilde{S}$ we use the well established LSUB$n$ approximation scheme, cf., e.g., Refs. Zeng et al. (1998); Farnell and Bishop (2004); Darradi et al. (2005, 2008); Farnell et al. (2009); Götze et al. (2011); Li et al. (2015); Bishop and Li (2017); Farnell et al. (2018). In the LSUB$n$ scheme all multi-spin correlations over distinct locales on the lattice defined by $n$ or fewer contiguous sites are retained. Using an efficient parallelized CCM code Farnell and Schulenburg we solve the CCM equations up to LSUB8. The maximum number of ket-state equations considered here is 128267. For the considered $(\pi,\pi,\pi)$ state of $J_{1}$–$J_{2}$ BCC model with FM $J_{1}$ we find that the LSUB$n$ data rapidly converge to the $n\to\infty$ limit. Thus, the difference betweeen the LSUB6 and LSUB8 CCM ground state energies $E_{0}$ (order parameter $M$) is less than 0.1% (1%). Therefore, the LSUB8 data used for the CCM curves shown in Fig. 9 practically may stand for the converged $n\to\infty$ data.

Appendix B Quantum Monte Carlo analysis

We have investigated the critical behavior of the model at a finite temperature for vanishing coupling constants $J_{2}=J_{3}=0$, considering a FM and AF $J_{1}$ interaction. Here, we fix $|J_{1}|=1$. We have simulated the model by means of the looper code Todo and Kato (2001); *looperweb of the ALPS library Albuquerque et al. (2007); Bauer et al. (2011); *alpsweb, for lattice sizes $L=8$, $12$, $16$, $24$, for a total of $2L^{3}$ lattice sites, and in an interval around the critical temperature. To compute the critical temperature we have performed a finite-size scaling Privman (1990); Campostrini et al. (2014) analysis of two renormalization-group invariant quantities. For the FM case we study the Binder ratio $U_{4}$ and the ratio $\xi/L$ of the second-moment correlation length $\xi$ over the lattice size $L$. The Binder ratio is defined as $U_{4}\equiv\langle M^{2}\rangle/\langle M^{4}\rangle$, where $M$ is the total magnetization of the system. The finite-size correlation length $\xi$ is defined in terms of the local magnetization $\hat{S}^{z}_{i}$. In the AF case, the order parameter is the staggered magnetization. Accordingly, we have analyzed the staggered Binder ratio $U_{4}^{\rm AF}\equiv\langle M_{s}^{2}\rangle/\langle M_{s}^{4}\rangle$, where $M_{s}\equiv\sum_{i}\epsilon_{i}\hat{S}^{z}_{i}$ is the staggered magnetization, and $\epsilon_{i}=1$ ($\epsilon_{i}=-1$) when the lattice site $i$ belongs to the $A$ ($B$) sublattice. As for the correlation length ratio, we have analyzed the quantity $\xi^{\rm AF}/L$, where the antiferromagnetic second-moment correlation length $\xi^{\rm AF}$ is defined in terms of the local staggered magnetization $\epsilon_{i}\hat{S}^{z}_{i}$. A discussion on the definition of a finite-size second-moment correlation length, in terms of a local order parameter, can be found in Appendix A of Ref. Parisen Toldin et al. (2015). In Fig. 13 we show the QMC estimates for the renormalization-group invariant observables considered.

Following Refs. Hasenbusch et al. (2008); Parisen Toldin et al. (2009), to analyze a renormalization-group invariant quantity $R$, we expand the corresponding scaling function and its leading scaling correction in a Taylor series around the critical temperature. To illustrate the procedure, we first consider the FM case and fit the Binder ratio $R=U_{4}$ to

[TABLE]

where $T_{c}$ is the critical temperature and $R^{*}$ is the universal value of $R$ at the critical point. Since the phase transition belongs to the classical three-dimensional Heisenberg universality class, here and in the following we fix the exponent $\nu$ to the corresponding value for such universality class $\nu=0.7112(5)$ Campostrini et al. (2002). In Eq. (10) we neglect scaling corrections. To monitor their influence, we have systematically disregarded the smallest lattice sizes. In Table 6 we report fit results as a function of $m_{\rm max}$ and the minimum lattice size taken into account $L_{\rm min}$. For a given value of $L_{\rm min}$ we observe a significant drop in the value of $\chi^{2}/\text{DOF}$ (DOF denotes the degrees of freedom) when increasing $m_{\rm max}$ from $m_{\rm max}=1$ to $m_{\rm max}=2$, whereas fits for $m_{\rm max}=3$ (not reported here) give a negligible improvement of $\chi^{2}/\text{DOF}$. This indicates that a second-order Taylor expansion of $R=U_{4}$ suitably describes the data, whereas a linear approximation is not sufficient. For fixed $m_{\rm max}=2$, on increasing $L_{\rm min}$ the value of $\chi^{2}/\text{DOF}$ decreases and we obtain a good value for $L_{\rm min}=12$, $16$. However, we also observe a systematic drift of the fitted values, which is larger than the statistical error bars. This clearly indicates that scaling corrections are relevant. To test their influence on the final results, we include them into the analysis, replacing Eq. (10) with

[TABLE]

where $k_{\rm max}$ is the Taylor expansion order of the correction to scaling term. Within the range and precision of our QMC data, fits of $R=U_{4}$ with $m_{\rm max}=2$ allow to include corrections-to-scaling with $k_{\rm max}=0$, i.e., to the leading order only, providing a suitable approximation of the scaling function and allowing to extract consistent results. As a further check, we have repeated the fits to Eq. (11) fixing the value of the correction-to-scaling exponent $\omega=0.8$, as expected for the three-dimensional $O(3)$ universality class Campostrini et al. (2002). A similar analysis has been done with the renormalization-group invariant ratio $R=\xi/L$. In this case we find that a linear approximation $m_{\rm max}=1$ is sufficient to fit the data. In Table 6 we report the results of our fit. Along these lines we have also analyzed the AF cases. As for the FM case, we have found that, within the range of our data and their precision, suitable fits of the Binder ratio $U_{4}^{\rm AF}$ require $m_{\rm max}=2$, whereas for $\xi^{\rm AF}/L$ a linear approximation $m_{\rm max}=1$ is sufficient. Corresponding fits are reported in Table 6. By judging conservatively the fit results, we extract the estimates of the critical temperatures for the FM and AF models reported in Table 4.

Fits in Table 6 allow also to determine the universal values $R^{*}$ of the renormalization-group invariant quantities at criticality. As expected by universality, they are the same for the FM and AF model. Inspecting the fit results of Table 6 we estimate

[TABLE]