Competing orders in pyrochlore magnets from a $\mathbb{Z}_2$ spin liquid perspective

TL;DR

This paper classifies $ ext{Z}_2$ spin-liquid states on the pyrochlore lattice, linking them to magnetic orders and predicting experimental signatures, thus advancing understanding of competing quantum phases in pyrochlore magnets.

Contribution

It provides a complete classification of spin-orbit-coupled $ ext{Z}_2$ spin liquids on the pyrochlore lattice using PSG, and constructs phase diagrams connecting spin liquids to magnetic orders.

Findings

Identifies intertwined magnetic and hidden orders in pyrochlore materials.

Classifies critical theories into $z=1$ and $z=2$ types with distinct experimental signatures.

Maps out phase diagrams linking spin liquids and magnetic orders.

Abstract

The pyrochlore materials have long been predicted to harbor a quantum spin liquid, an intrinsic long-range-entangled state supporting fractionalized excitations. Existing pyrochlore experiments, on the other hand, have discovered several weakly ordered states and a tendency of close competition amongst them. Motivated by these facts, we give a complete classification of spin-orbit-coupled $\mathbb{Z}_2$ spin-liquid states on the pyrochlore lattice by using the projective symmetry group (PSG) approach for bosonic spinons. For each spin liquid, we construct a mean-field Hamiltonian that can be used to describe phase transitions out of the spin liquid via spinon condensation. Studying these phase transitions, we establish phase diagrams for our mean-field Hamiltonians that link magnetic orders to specific spin liquids. In general, we find that seemingly unrelated magnetic orders are intertwined with each other and that the conventional spin orders seen in the experiments are accompanied by more exotic hidden orders. Our critical theories are categorized into $z=1$ and $z=2$ types, based on their spinon dispersion and Hamiltonian diagonalizability, and are shown to give distinct signatures in the heat capacity and the spin structure factor. This study provides a clear map of pyrochlore phases for future experiments and variational Monte Carlo studies in pyrochlore materials.

Tables

Table 1. Table 1: Most important characteristics of the 15 15 15 critical “paraphases” corresponding to phase transitions between the six 0 0 -flux ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} spin-liquid phases [labeled as 0 0 - ( n C ¯ 6 ​ S ​ n S ​ T 1 ​ n C ¯ 6 ) subscript 𝑛 subscript ¯ 𝐶 6 𝑆 subscript 𝑛 𝑆 subscript 𝑇 1 subscript 𝑛 subscript ¯ 𝐶 6 (n_{\overline{C}_{6}S}\,n_{ST_{1}}\,n_{\overline{C}_{6}}) ] and a rich variety of magnetically ordered phases. Each paraphase is labeled by the condensation momenta (see Table IV for notation) where the spinons become gapless and condense at the phase transition. For each critical theory, experimental signatures are given in terms of the dynamical critical exponent, the low-temperature behavior of the heat capacity, and the low-energy features of the dynamic spin structure factor (see Fig.  4 ). For each magnetically ordered phase obtained by spinon condensation, it is specified whether various orders are generically present ( + + ) or absent ( 0 0 ), including conventional spin orders, such as all-in-all-out (AIAO) order, XY antiferromagnetic (AFM) order [ Ψ 2 , 3 subscript Ψ 2 3 \Psi_{2,3} ], ferromagnetic (FM) order, and Palmer-Chalker (PC) order [ Ψ 4 subscript Ψ 4 \Psi_{4} ], as well as inversion-breaking “hidden” orders.

Spin-	Critical “paraphases”				Magnetically ordered phases
liquid	Condensation	Dynamical	Heat capacity:	Dynamic spin structure factor	Spin orders				Hidden
phases	momenta	exponent	$C_V\propto T^x$		AIAO	AFM	FM	PC	orders
$0$-$\left(001\right)$	$\Gamma$	$z=2$	$x=\frac{3}{2}$	Gapless at $\Gamma$	$0$	$0$	$+$	$0$	$+$
				Weak in the low-energy limit
	$L$	$z=2$	$x=\frac{3}{2}$	Gapless at $\Gamma$ and $X$	$+$	$+$	$+$	$+$	$+$
				Weak in the low-energy limit
	$\Lambda$	$z=2$	$x=1$	Gapless along $\Gamma \to X$ and $K\to \Gamma \to L\to U$	Unclear at NN level				$+$
				Low-energy weight at all momenta
$0$-$\left(010\right)$	$\Gamma$	$z=2$	$x=\frac{3}{2}$	Gapless at $\Gamma$	$0$	$0$	$+$	$+$	$0$
				Weak in the low-energy limit
	$\Lambda$	$z=1$	$x=2$	Gapless along $\Gamma \to X$ and $K\to \Gamma \to L\to U$	Unclear at NN level				$0$
$0$-$\left(100\right)$	$\Gamma$	$z=2$	$x=\frac{3}{2}$	Gapless at $\Gamma$	$+$	$+$	$0$	$0$	$0$
				Weak in the low-energy limit
	$\Lambda$	$z=1$	$x=2$	Gapless along $\Gamma \to X$ and $K\to \Gamma \to L\to U$	Unclear at NN level				$0$
$0$-$\left(101\right)$	$\Gamma$	$z=2$	$x=\frac{3}{2}$	Gapless at $\Gamma$	$+$	$+$	$0$	$0$	$+$
				Weak in the low-energy limit
	$W$	$z=1$	$x=3$	Gapless at $\Gamma$, $X$, and $\frac{2}{3}K$	$+$	$+$	$+$	$+$	$+$
				Singular in the low-energy limit
	$X$	$z=1$	$x=3$	Gapless at $\Gamma$ and $X$	$+$	$+$	$0$	$+$	$+$
				Singular in the low-energy limit
$0$-$\left(110\right)$	$\Gamma$	$z=1$	$x=3$	Gapless at $\Gamma$	$0$	$+$	$0$	$0$	$0$
				Characteristic lower edge of the spectrum
	$\Lambda$	$z=1$	$x=2$	Gapless along $\Gamma \to X$ and $K\to \Gamma \to L\to U$	Unclear at NN level				$0$
$0$-$\left(111\right)$	$\Gamma$	$z=2$	$x=\frac{3}{2}$	Gapless at $\Gamma$	$+$	$+$	$0$	$0$	$+$
				Weak in the low-energy limit
	$W$	$z=1$	$x=3$	Gapless at $\Gamma$, $X$, and $\frac{2}{3}K$	$+$	$+$	$+$	$+$	$+$
				Singular in the low-energy limit
	$X$	$z=1$	$x=3$	Gapless at $\Gamma$ and $X$	$+$	$+$	$+$	$0$	$+$
				Singular in the low-energy limit

Table 2. Table 2: Independent mean-field parameters and constraints for the sixteen PSG classes. The parameters not mentioned in this table are enforced to be zero by the constraints. The mean-field Hamiltonians for some PSG classes appear to be U(1) on the NN level, as indicated by the comment “U(1) at NN”, but recover their ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} character upon including NNN terms. Note that a nonzero onsite chemical potential μ = α 𝜇 𝛼 \mu=\alpha is allowed in all PSG classes.

Class	Independent nonzero parameters			Constraints			Note
$n_1\pi$-$\left(n_{{\overline{C}}_{6}S}{n}_{S{T}_1}{n}_{{\overline{C}}_6}\right)$	Onsite	NN	NNN	Onsite	NN	NNN
0-$\left(000\right)$	$\mu$	$a$,$c$	$A,B,D,B^{\prime }$		$c=-d$	$B=C$, $B^{\prime }=-C^{\prime }$	U(1) at NN
0-$\left(001\right)$	$\mu$	$a$,$c$,$b^{\prime }$	$A,B,D,B^{\prime }$		$c=-d$,	$B=C$, $B^{\prime }=-C^{\prime }$
0-$\left(010\right)$	$\mu$, $\nu$	$a$,$c$,$b^{\prime }$	$A,B,D,A^{\prime },B^{\prime },D^{\prime }$	${\beta }^{\prime }={\delta }^{\prime }={\gamma }^{\prime }\equiv \nu$	$c=-d$,	$B=C$, $B^{\prime }=C^{\prime }$
0-$\left(011\right)$	$\mu$,	$a$,$c$	$A,B,D,A^{\prime },B^{\prime },D^{\prime }$		$c=-d$,	$B=C$, $B^{\prime }=C^{\prime }$	U(1) at NN
0-$\left(100\right)$	$\mu$,	$b$,$c^{\prime }$	$A,B,D,B^{\prime }$		$c^{\prime }=d^{\prime }$,	$B=C$, $B^{\prime }=-C^{\prime }$
0-$\left(101\right)$	$\mu$,	$b$,$a^{\prime }$,$c^{\prime }$	$A,B,D,B^{\prime }$		$c^{\prime }=-d^{\prime }$,	$B=C$, $B^{\prime }=-C^{\prime }$
0-$\left(110\right)$	$\mu$, $\nu$	$b$,$a^{\prime }$,$c^{\prime }$	$A,B,D,A^{\prime },B^{\prime },D^{\prime }$	${\beta }^{\prime }={\delta }^{\prime }={\gamma }^{\prime }\equiv \nu$	$c^{\prime }=-d^{\prime }$,	$B=C$, $B^{\prime }=C^{\prime }$
0-$\left(111\right)$	$\mu$,	$b$,$c^{\prime }$	$A,B,D,A^{\prime },B^{\prime },D^{\prime }$		$c^{\prime }=d^{\prime }$,	$B=C$, $B^{\prime }=C^{\prime }$
$\pi$-$\left(000\right)$	$\mu$, $\nu$	$a$,$c$,$b^{\prime }$	$B,B^{\prime }$	${\beta }^{\prime }={\delta }^{\prime }={\gamma }^{\prime }=\nu$	$c=-d$	$B=-C$, $B^{\prime }=-C^{\prime }$
$\pi$-$\left(001\right)$	$\mu$	$a$,$c$	$B,B^{\prime }$		$c=-d$	$B=-C$, $B^{\prime }=-C^{\prime }$	U(1) at NN
$\pi$-$\left(010\right)$	$\mu$	$a$,$c$	$B,A^{\prime },B^{\prime },D^{\prime }$		$c=-d$	$B=-C$, $B^{\prime }=C^{\prime }$
$\pi$-$\left(011\right)$	$\mu$,	$a$,$c$,$b^{\prime }$	$B,A^{\prime },B^{\prime },D^{\prime }$		$c=-d$	$B=-C$, $B^{\prime }=C^{\prime }$
$\pi$-$\left(100\right)$	$\mu$, $\nu$	$b$,$a^{\prime }$,$c^{\prime }$	$B,B^{\prime }$	${\beta }^{\prime }={\delta }^{\prime }={\gamma }^{\prime }=\nu$	$c^{\prime }=-d^{\prime }$	$B=-C$, $B^{\prime }=-C^{\prime }$	U(1) at NN
$\pi$-$\left(101\right)$	$\mu$,	$b$,$c^{\prime }$	$B,B^{\prime }$		$c^{\prime }=d^{\prime }$	$B=-C$, $B^{\prime }=-C^{\prime }$
$\pi$-$\left(110\right)$	$\mu$	$b$,$c^{\prime }$	$B,A^{\prime },B^{\prime },D^{\prime }$		$c^{\prime }=d^{\prime }$	$B=-C$, $B^{\prime }=C^{\prime }$
$\pi$-$\left(111\right)$	$\mu$,	$b$,$a^{\prime }$,$c^{\prime }$	$B,A^{\prime },B^{\prime },D^{\prime }$		$c^{\prime }=-d^{\prime }$	$B=-C$, $B^{\prime }=C^{\prime }$

Table 3. Table 3: Paraphase boundaries of the NN ansätze.

Class	Adjacency	Paraphase boundary
0-$\left(001\right)$	$\Lambda$ vs $L$	$2ac+c^2-b^{\prime 2}=0$
	$\Gamma$ vs $L$	$2a^2+ac-b^{\prime 2}=0$ for $c>a$ and
		$ac-4c^2+b^{\prime 2}=0$ for
	$\Lambda$ vs $\Gamma$	$2a+c=b^{\prime }=0$
0-$\left(010\right)$	$\Gamma$ vs $\Lambda$	$\sqrt{4{\left(a-c\right)}^2+3{\left(\nu -b^{\prime }\right)}^2}$
		$=-2a-4c+\sqrt{3}\left|\nu +b^{\prime }\right|$
0-$\left(100\right)$	$\Gamma$ vs $\Lambda$	$\sqrt{2}b=\pm c^{\prime }$
0-$\left(101\right)$	$\Gamma$ vs $W$	$-2b^2+3c^{\prime 2}+2c^{\prime }{a}^{\prime }+a^{\prime 2}=0$, $b>0$
	$X$ vs $W$	$-2b^2+3c^{\prime 2}+2c^{\prime }{a}^{\prime }+a^{\prime 2}=0$,
0-$\left(110\right)$	$\Gamma$ vs $\Lambda$	$\sqrt{{\left(\nu +a^{\prime }-2c^{\prime }\right)}^2+2{\left(\nu -a^{\prime }\right)}^2}$
		$=4b+\sqrt{3}\left|\nu +a^{\prime }+2c^{\prime }\right|$
0-$\left(111\right)$	$\Gamma$ vs $W$	$\sqrt{2}b=\pm c^{\prime }$
	$X$ vs $W$	$b=\pm \sqrt{2}{c}^{\prime }$

Table 4. Table 4: Possible sets of condensation momenta.

Label	Description
$\Gamma$	$(0,0,0)$
$L$	$\pi ({\delta }_1,{\delta }_2,{\delta }_3)$, where ${\delta }_1,{\delta }_2,{\delta }_3\in \{1,-1\}$
$\Lambda$	$k({\delta }_1,{\delta }_2,{\delta }_3)$, where $k\in [-\pi ,\pi ]$ and ${\delta }_1,{\delta }_2,{\delta }_3\in \{1,-1\}$
$X$	$X^1=2\pi (1,0,0)$, $X^2=2\pi (0,1,0)$, $X^3=2\pi (0,0,1)$
$W$	$\pi (2,\pm 1,0)$ and all permutations of the $3$ components

Table 5. Table 5: Parametrization of the NN mean-field ansätze using generalized spherical coordinates ( ψ , θ , ϕ ) 𝜓 𝜃 italic-ϕ (\psi,\theta,\phi) for the phase diagram in Fig.  1 .

Class	Independent nonzero parameters up to NN terms
	and parameterized by
0-(001)	$(a,c,b^{\prime })=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta )$
0-(010)	$(\nu ,a,c,b^{\prime })=$
	$(\cos \psi ,\sin \psi \sin \theta \cos \varphi ,\sin \psi \sin \theta \sin \varphi ,\sin \psi \cos \theta )$
0-(100)	$(b,c^{\prime })=(\cos \varphi ,\sin \varphi )$
0-(101)	$(b,a^{\prime },c^{\prime })=(\cos \theta ,\sin \theta \cos \varphi ,\sin \theta \sin \varphi )$
0-(110)	$(\nu ,b,a^{\prime },c^{\prime })=$
	$(\cos \psi ,\sin \psi \cos \theta ,\sin \psi \sin \theta \cos \varphi ,\sin \psi \sin \theta \sin \varphi )$
0-(111)	$(b,c^{\prime })=(\cos \varphi ,\sin \varphi )$

Table 6. Table 6: Power-law exponents of the low-energy spinon density of states and the corresponding low-temperature heat capacity for critical theories of dynamical exponents z = 1 , 2 𝑧 1 2 z=1,2 where spinons condense at points or along lines.

Condensation	Dynamical	Density of states:	Heat capacity:
manifold	exponent	$g\left(ϵ\right)\propto {ϵ}^{\alpha }$	$C_V\propto T^x$
Point(s)	$z=2$	$\alpha =\frac{1}{2}$	$x=\frac{3}{2}$
	$z=1$	$\alpha =2$	$x=3$
Line(s)	$z=2$	$\alpha =0$	$x=1$
	$z=1$	$\alpha =1$	$x=2$

Table 7. Table 7: One-to-one correspondence between the potential set of momenta at which the spinons condense at the critical point and the set of momenta at which the corresponding dynamic structure factor is gapless along the high-symmetry path in Fig.  4 ; these set of momenta can be points A A \mathrm{A} or sections A → B → A B \mathrm{A}\rightarrow\mathrm{B} between two points A A \mathrm{A} and B B \mathrm{B} .

Spinon condensation	Gapless points or regions in
momenta	dynamic structure factor
$\Gamma$	$\Gamma$
$L$	$\Gamma$, $X$
$X$	$\Gamma$, $X$
$W$	$\Gamma$, $X$, $\frac{2}{3}K$
$\Lambda$	$\Gamma \to X$, $K\to \Gamma \to L\to U$

Table 8. Table 8: Analysis of the zero-momentum ( K → = 0 → 𝐾 0 \vec{K}=0 ) order parameters for the paraphases characterized by single-point and multi-point condensation (i.e., excluding line condensation). For each paraphase, the complex ( ϕ italic-ϕ \phi ) or real ( χ 𝜒 \chi ) condensation fields are specified; the zero-momentum order parameters are then bilinears of these fields with total momentum K → = 0 → 𝐾 0 \vec{K}=0 and transform under various irreducible representations of the point group O h subscript 𝑂 h O_{\text{h}} . For some paraphases, distinct order parameters are intertwined such that they must appear together at condensation; for each paraphase characterized by single-point condensation, the order parameters are arranged into classes C 𝐶 C (marked by curly brackets) such that there must be at least one nonzero order parameter from each class C 𝐶 C . The number of quadratic scalars in terms of the order parameters is also specified; if there is only one such scalar, all order parameters are in different classes and hence are maximally intertwined.

Paraphase	Condensation	Real dimension of	Decomposition into	Number of scalars	General intertwining between
	fields	order-parameter	irreducible representations	quadratic in order	irreducible representations
		space at $\overrightarrow{K}=0$	(i.e., distinct orders)	parameters	$C_1\oplus \cdots \oplus C_q$
$0$-$\left(001\right)\Gamma$	${\varphi }_{1,2}$	10	$A_{1g}\oplus T_{1g}\oplus 2T_{2u}$	$1$	$\left\{A_{1g}\right\}\oplus \left\{T_{1g}\right\}\oplus \left\{T_{2u}\right\}$
$0$-$\left(001\right)L$	${\varphi }_{1-8}$	40	$A_{1g}\oplus A_{2g}\oplus E_g\oplus 2T_{1g}\oplus 2T_{2g}$
	($2$ per ${\overrightarrow{k}}_c$)		$\oplus 2A_{1u}\oplus 2E_u\oplus 2T_{1u}\oplus 4T_{2u}$
$0$-$\left(010\right)\Gamma$	${\varphi }_{1,2}$	10	$A_{1g}\oplus T_{1g}\oplus 2T_{2g}$	$1$	$\left\{A_{1g}\right\}\oplus \left\{T_{1g}\right\}\oplus \left\{T_{2g}\right\}$
$0$-$\left(100\right)\Gamma$	${\varphi }_{1,2}$	10	$A_{1g}\oplus 3A_{2g}\oplus 3E_g$	$6$	$\left\{A_{1g}\right\}\oplus \left\{A_{2g}\right\}\oplus \left\{E_g\right\}$
$0$-$\left(101\right)\Gamma$	${\varphi }_{1,2}$	10	$A_{1g}\oplus A_{2g}\oplus E_g\oplus 2A_{1u}\oplus 2E_u$	$4$	$\left\{A_{1g}\right\}\oplus \{A_{2g},E_g\}\oplus \{A_{1u},E_u\}$
$0$-$\left(101\right)W$	${\chi }_{1-24}$	48	$2A_{1g}\oplus A_{2g}\oplus 3E_g\oplus T_{1g}$
	($8$ per $\pm {\overrightarrow{k}}_c$)		$\oplus 2T_{2g}\oplus 2A_{1u}\oplus A_{2u}$
			$\oplus 3E_u\oplus 3T_{1u}\oplus 4T_{2u}$
$0$-$\left(101\right)X$	${\chi }_{1-12}$	30	$A_{1g}\oplus 2A_{2g}\oplus 3E_g\oplus T_{2g}\oplus 2A_{1u}$
	($4$ per ${\overrightarrow{k}}_c$)		$\oplus A_{2u}\oplus 3E_u\oplus T_{1u}\oplus 2T_{2u}$
$0$-$\left(110\right)\Gamma$	${\chi }_{1,2}$	3	$A_{1g}\oplus E_g$	$1$	$\left\{A_{1g}\right\}\oplus \left\{E_g\right\}$
$0$-$\left(111\right)\Gamma$	${\varphi }_{1,2}$	10	$A_{1g}\oplus A_{2g}\oplus E_g\oplus 2A_{2u}\oplus 2E_u$	$4$	$\left\{A_{1g}\right\}\oplus \{A_{2g},E_g\}\oplus \{A_{2u},E_u\}$
$0$-$\left(111\right)W$	${\chi }_{1-24}$	48	$A_{1g}\oplus 2A_{2g}\oplus 3E_g\oplus 2T_{1g}$
	($8$ per $\pm {\overrightarrow{k}}_c$)		$\oplus T_{2g}\oplus A_{1u}\oplus 2A_{2u}$
			$\oplus 3E_u\oplus 4T_{1u}\oplus 3T_{2u}$
$0$-$\left(111\right)X$	${\chi }_{1-12}$	30	$2A_{1g}\oplus A_{2g}\oplus 3E_g\oplus T_{1g}\oplus A_{1u}$
	($4$ per ${\overrightarrow{k}}_c$)		$\oplus 2A_{2u}\oplus 3E_u\oplus 2T_{1u}\oplus T_{2u}$

Table 9. Table 9: Irreducible representations of the point group O h subscript 𝑂 h O_{\mathrm{h}} and the corresponding symmetry-breaking orders. For some representations, simple examples of order parameters are provided in terms of the spins S → i subscript → 𝑆 𝑖 \vec{S}_{i} at sites i 𝑖 i , where r → i subscript → 𝑟 𝑖 \vec{r}_{i} is the vector from site i 𝑖 i to the center of the nearest “up” tetrahedron, n → i , j subscript → 𝑛 𝑖 𝑗 \vec{n}_{i,j} is the vector from site i 𝑖 i to site j 𝑗 j , and λ i , j = ± 1 subscript 𝜆 𝑖 𝑗 plus-or-minus 1 \lambda_{i,j}=\pm 1 for bonds ⟨ i , j ⟩ 𝑖 𝑗 \langle i,j\rangle in “up” and “down” tetrahedra, respectively. Note that the scalar representation A 1 ​ g subscript 𝐴 1 g A_{1\mathrm{g}} does not break any symmetries and hence does not correspond to any order.

Irrep	Dim.	Standard name of	Simple example of order
		corresponding order	parameter in terms of spins
$A_{1g}$	$1$	(N/A)	$1$
$A_{2g}$	$1$	All-in-all-out	${\sum }_i{\overrightarrow{r}}_i\cdot {\overrightarrow{S}}_i$
$E_g$	$2$	XY antiferromagnet
$T_{1g}$	$3$	Ferromagnet	${\sum }_i{\overrightarrow{S}}_i$
$T_{2g}$	$3$	Palmer-Chalker	${\sum }_i{\overrightarrow{r}}_i\times {\overrightarrow{S}}_i$
$A_{1u}$	$1$		${\sum }_{⟨i,j⟩}{\lambda }_{i,j}\left({\overrightarrow{S}}_i\cdot {\overrightarrow{S}}_j\right)$
$A_{2u}$	$1$		${\sum }_{⟨i,j⟩}{\overrightarrow{n}}_{i,j}\cdot \left({\overrightarrow{S}}_i\times {\overrightarrow{S}}_j\right)$
$E_u$	$2$
$T_{1u}$	$3$		${\sum }_{⟨i,j⟩}\left({\overrightarrow{r}}_i\times {\overrightarrow{n}}_{i,j}\right)\times \left({\overrightarrow{S}}_i\times {\overrightarrow{S}}_j\right)$
$T_{2u}$	$3$		${\sum }_{⟨i,j⟩}{\overrightarrow{n}}_{i,j}\times \left({\overrightarrow{S}}_i\times {\overrightarrow{S}}_j\right)$

Table 10. Table 10: Correspondence between orders, irreps and basis of irreps

Irrep	Basis	Orders
$A_2$	${𝐒}_1$	all in-all out
$E$	${𝐒}_2,{𝐒}_3$	${\Psi }_2$ and ${\Psi }_3$
$T_{1,A}$	${𝐒}_4,{𝐒}_5,{𝐒}_6$	collinear FM
$T_{1,B}$	${𝐒}_7,{𝐒}_8,{𝐒}_9$	non-collinear FM
$T_2$	${𝐒}_{10},{𝐒}_{11},{𝐒}_{12}$	Palmer-Chalker

Table 11. Table 11: Critical chemical potential μ 𝜇 \mu for the 15 paraphases.

Paraphase	Critical ${\mu }_c$
0-(001)$\Gamma$	${\mu }_c=\max \{-6a,2a-8c\}$
0-(001)$L$	Largest root of ${\mu }^3+2\left(a+2c\right){\mu }^2$
	$-4\left(2a^2-4ac+2c^2+3b^{\prime 2}\right)\mu -24b^{\prime 2}\left(a+2c\right)=0$
0-(001)$\Lambda$	${\mu }_c=2a+4c$
0-(010)$\Gamma$	${\mu }_c=2\left(a+2c+\sqrt{4{\left(a-c\right)}^2+3{\left(\nu -b^{\prime }\right)}^2}\right)$
0-(010)$\Lambda$	${\mu }_c=2\left(-a-2c+\sqrt{3}\left|\nu +b^{\prime }\right|\right)$
0-(100)$\Gamma$	${\mu }_c=-6b$
0-(100)$\Lambda$	${\mu }_c=2b+4\sqrt{2}\left|c^{\prime }\right|$
0-(101)$\Gamma$	${\mu }_c=-6b$
0-(101)$W$	Largest root of ${\mu }^4-8\left(b^2+2a^{\prime 2}+4c^{\prime 2}\right){\mu }^2$
	$+64\left(2a^{\prime }+c^{\prime }\right)b{c}^{\prime }\mu -64a^{\prime }{c}^{\prime }\left(b^2+3c^{\prime 2}\right)$
	$+16{\left(b^2-3c^{\prime 2}\right)}^2-32b^2{a}^{\prime 2}+48a^{\prime 4}=0$
0-(101)$X$	${\mu }_c=2b+2\sqrt{2}\left|a^{\prime }-c^{\prime }\right|$
0-(110)$\Gamma$	${\mu }_c=6b+2\sqrt{3}\left|\nu +a^{\prime }+2c^{\prime }\right|$
0-(110)$\Lambda$	${\mu }_c=-2b+2\sqrt{{\left(\nu +a^{\prime }-2c^{\prime }\right)}^2+2{\left(\nu -a^{\prime }\right)}^2}$
0-(111)$\Gamma$	${\mu }_c=6b$
0-(111)$W$	${\mu }_c=\max \left\{\pm \sqrt{2}w+\sqrt{2}\sqrt{2b^2\mp 4\sqrt{2}bw+7w^2}\right\}$
0-(111)$X$	${\mu }_c=-2b+2\sqrt{6}\left|c^{\prime }\right|$

Equations

$${\hat{e}}_1$$

$${\hat{e}}_2$$

$${\hat{e}}_3$$

$$(r_1,r_2,r_3{)}_{\mu }$$

$$=\frac{1}{2}(r_2+r_3,r_3+r_1,r_1+r_2)+\frac{1}{2}{\hat{e}}_{\mu }.$$

$$T_1:$$

$$T_i{T}_{i+1}{T}_i^{-1}{T}_{i+1}^{-1}$$

$${\overline{C}}_6^6$$

$$S^2{T}_3^{-1}$$

$${\overline{C}}_6{T}_i{\overline{C}}_6^{-1}{T}_{i+1}$$

$$S{T}_i{S}^{-1}{T}_3^{-1}{T}_i$$

$$S{T}_3{S}^{-1}{T}_3^{-1}$$

$$({\overline{C}}_{6}S{)}^4$$

$$({\overline{C}}_6^{3}S{)}^2$$

$${\mathcal{T}}^2$$

$$\mathcal{T}\mathcal{O}{\mathcal{T}}^{-1}{\mathcal{O}}^{-1}$$

$${\hat{S}}_{{\vec{r}}_{\mu }}^{\alpha }=\frac{1}{2}{b}_{{\vec{r}}_{\mu }}^{†}{\sigma }^{\alpha }{b}_{{\vec{r}}_{\mu }}^{},\alpha =x,y,z,$$

$$G:b_{{\vec{r}}_{\mu }}\to e^{i\varphi ({\vec{r}}_{\mu })}{b}_{{\vec{r}}_{\mu }}$$

$$\sum _{\sigma =\uparrow ,\downarrow }{b}_{{\vec{r}}_{\mu },\sigma }^{†}{b}_{{\vec{r}}_{\mu },\sigma }^{}=1$$

$$\mathcal{O}:b_{{\vec{r}}_{\mu }}\to U_{\mathcal{O}}^{†}{b}_{\mathcal{O}({\vec{r}}_{\mu })}.$$

$$G_{\mathcal{O}}:b_{{\vec{r}}_{\mu }}\to e^{i{\varphi }_{\mathcal{O}}({\vec{r}}_{\mu })}{b}_{{\vec{r}}_{\mu }},$$

$$\tilde{\mathcal{O}}=G_{\mathcal{O}}\circ \mathcal{O}:b_{{\vec{r}}_{\mu }}\to e^{i{\varphi }_{\mathcal{O}}[\mathcal{O}({\vec{r}}_{\mu })]}{U}_{\mathcal{O}}^{†}{b}_{O({\vec{r}}_{\mu })},$$

$$\mathcal{T}:b_{{\vec{r}}_{\mu }}\to \mathcal{K}{U}_{\mathcal{T}}^{†}{b}_{{\vec{r}}_{\mu }},$$

$$G_{\mathcal{T}}:b_{{\vec{r}}_{\mu }}\to e^{i{\varphi }_{\mathcal{T}}({\vec{r}}_{\mu })}{b}_{{\vec{r}}_{\mu }},$$

$$\tilde{\mathcal{T}}=G_{\mathcal{T}}\circ \mathcal{T}:b_{{\vec{r}}_{\mu }}\to e^{i{\varphi }_{\mathcal{T}}({\vec{r}}_{\mu })}\mathcal{K}{U}_{\mathcal{T}}^{†}{b}_{{\vec{r}}_{\mu }}.$$

$${\mathcal{O}}_1\circ {\mathcal{O}}_2\circ \cdots =1,$$

$${\tilde{\mathcal{O}}}_1\circ {\tilde{\mathcal{O}}}_2\circ \cdots =(G_{{\mathcal{O}}_1}\circ {\mathcal{O}}_1)\circ (G_{{\mathcal{O}}_2}\circ {\mathcal{O}}_2)\circ \cdots =\mathcal{G},$$

$$G_{{\mathcal{O}}_1}$$

$${\mathcal{O}}_i\circ G_{{\mathcal{O}}_j}\circ {\mathcal{O}}_i^{-1}:b_{{\vec{r}}_{\mu }}\to e^{i{\varphi }_{{\mathcal{O}}_j}[{\mathcal{O}}_i^{-1}({\vec{r}}_{\mu })]}{b}_{{\vec{r}}_{\mu }},$$

$${\varphi }_{{\mathcal{O}}_1}({\vec{r}}_{\mu })$$

$$\mathcal{T}\circ G_{\mathcal{O}}\circ {\mathcal{T}}^{-1}:$$

$$b_{{\vec{r}}_{\mu }}\to \mathcal{K}{U}_{\mathcal{T}}^{†}{e}^{i{\varphi }_{\mathcal{O}}({\vec{r}}_{\mu })}{U}_{\mathcal{T}}{\mathcal{K}}^{†}{b}_{{\vec{r}}_{\mu }}=e^{-i{\varphi }_{\mathcal{O}}({\vec{r}}_{\mu })}{b}_{{\vec{r}}_{\mu }},$$

$${\varphi }_{\mathcal{T}}({\vec{r}}_{\mu })-{\varphi }_{\mathcal{T}}[{\mathcal{O}}^{-1}({\vec{r}}_{\mu })]-2{\varphi }_{\mathcal{O}}({\vec{r}}_{\mu })=n\pi \mathrm{mod}2\pi ,$$

$$(G\circ G_{{\mathcal{O}}_1}\circ {\mathcal{O}}_1\circ G^{-1})\circ (G\circ G_{{\mathcal{O}}_2}\circ {\mathcal{O}}_2\circ G^{-1})\circ \cdots =\mathcal{G},$$

$${\varphi }_{{\mathcal{O}}_i}({\vec{r}}_{\mu })\to {\varphi }_{{\mathcal{O}}_i}({\vec{r}}_{\mu })+\varphi ({\vec{r}}_{\mu })-\varphi [{\mathcal{O}}_i^{-1}({\vec{r}}_{\mu })],$$

$${\varphi }_{T_1}({\vec{r}}_{\mu })$$

$${\varphi }_{T_2}({\vec{r}}_{\mu })$$

Full text

I Introduction

Quantum spin liquids (QSLs) Savary and Balents (2016) are zero-temperature phases of interacting spin systems which possess intrinsic long-range entanglement and support nonlocal excitations carrying fractionalized quantum numbers. Typically, they respect all symmetries of the underlying lattice, i.e., they exhibit a lack of conventional symmetry-breaking order. The theoretical understanding of QSLs is largely in terms of emergent gauge theory, which provides a convenient mathematical framework to describe long-range entanglement, along with the nonlocal nature of the fractionalized excitations.

In frustrated magnetic systems Balents (2010), QSL ground states may control the physics even at (small) finite temperatures, as long as energy dominates over entropy. For two-dimensional spin liquids, this statement is purely asymptotic; at any nonzero temperature $T>0$, the putative QSL is adiabatically connected to a high-temperature paramagnet. However, some three-dimensional spin liquids, particularly the so-called $\mathbb{Z}_{2}$ states with Ising-like emergent gauge fields, are more robust, and can persist in the form of a distinct low-temperature phase up to a nonzero critical temperature.

While QSLs are extremely interesting from a conceptual perspective, it is far from obvious to realize them in experimental materials, or even realistic spin Hamiltonians. Traditionally, most studies considered spin-rotation-invariant Heisenberg systems on geometrically frustrated two-dimensional lattices. However, it has recently been recognized that magnetic systems with strong spin-orbit coupling provide a promising alternative avenue to QSLs Witczak-Krempa et al. (2014); Gingras and McClarty (2014); Rau et al. (2016a); Iaconis et al. (2018). In general, these systems have a large number of magnetically anisotropic terms, leading to exchange frustration as well as an extended parameter space, and are thus expected to harbor QSL ground states on a wide range of two- and three-dimensional lattices.

The most widely studied such three-dimensional structure is the pyrochlore lattice, consisting of periodically arranged corner-sharing tetrahedra. Experimentally, two large families of materials, the pyrochlore spinels and the rare-earth pyrochlores, provide vast real-world possibilities Gardner et al. (2010) to test theoretical predictions on the pyrochlore lattice. In the 2000s, it was predicted that certain antiferromagnetic pyrochlore models could support a U(1) QSL phase Hermele et al. (2004) [the “U(1)” means that the gauge field belongs to the Lie algebra of the U(1) group and that the emergent charges are characterized by integers related to the generating charge of U(1)], which is a simulacrum of electromagnetic gauge theory in high-energy physics. In 2011/2012, theoretical applications of this idea to realistic models emerged, suggesting the presence of a U(1) spin liquid in the so-called “quantum spin ice” pyrochlore materials Ross et al. (2011); Savary and Balents (2012). So far, these predictions remain to be confirmed in experiments, even though there are some promising recent developments Hirschberger et al. (2015, 2019); Gaudet et al. (2019); Sibille et al. (2016, 2015); Scheie et al. (2017); Thompson et al. (2017); Sibille et al. (2018); Tokiwa et al. (2018).

Another thread recurring in the experimental study of rare-earth pyrochlores is the close competition amongst several weakly ordered states Hallas et al. (2017). Several hints at this competition are present in the family of Yb pyrochlores, Yb2B2O7, which have a systematic structural evolution across the series B = Ge, Ti, Pt, Sn. While the germanate orders antiferromagnetically, the remaining members of the family have ferromagnetic ground states, suggesting the close proximity of at least these two phases. In each material, the specific heat is peaked at a temperature of 2-4K, while the maximum ordering temperature is 0.6K in the germanate and half or less than that in the rest of the family. These findings indicate the onset of strong spin correlations well above the ordering temperature, but an inability of the system to decide upon its ground state. The weak ferromagnetic ground state in Yb2Ti2O7 is also famously mercurial, changing its character substantially with sample variations Bowman et al. (2019). Theoretically, a classical analysis indeed finds close competition amongst several distinct phases Yan et al. (2017), but a quantum picture of this phase competition is not yet available.

In this work, we combine the two threads of phase competition and QSL physics by utilizing the connection of symmetry to emergent gauge structure. This connection is mathematically described by the projective symmetry group (PSG), proposed by Wen in 2002 Wen (2002), which encapsulates the fact that, in a QSL, the group operations of the physical symmetry group are interleaved with those of the emergent gauge group. The embedding of the physical symmetries into the PSG can then lead to a unification of distinct symmetry-breaking orders that are unrelated in classical physics. Such a unified description of seemingly unrelated magnetic orders is the main motivation behind the present study.

The PSG also offers a straightforward method to classify QSLs in the presence of symmetry. Concretely, the PSG specifies a distinct set of transformation rules for the emergent matter and gauge fields in each QSL phase, corresponding to a given PSG class. Employing the PSG method, an entire zoo of QSLs has been found on the square Reuther et al. (2014), triangular Lu (2016), kagome Lu et al. (2011), honeycomb You et al. (2012), star Choy and Kim (2009), and hyperkagome Huang et al. (2017) lattices, to give a few notable examples. Generally, these QSLs can be connected to magnetically ordered states by considering the condensation patterns that emerge when the energy of a bosonic QSL excitation is brought to zero Bergman et al. (2006); Chen (2016); Li and Chen (2017).

In this paper, we employ the PSG method to obtain a full classification of QSLs with $\mathbb{Z}_{2}$ gauge structure on the pyrochlore lattice using Schwinger bosons Sachdev (1992); Wang and Vishwanath (2006); Wang (2010); Yang and Wang (2016). While standard parton constructions also allow U(1) and SU(2) gauge structures, we consider the $\mathbb{Z}_{2}$ gauge structure for two reasons. First, it is the simplest one: quasiparticles in a $\mathbb{Z}_{2}$ QSL are weakly interacting because the gauge field itself is gapped. Second, it is also the richest one: a single U(1) PSG class can be further split into several $\mathbb{Z}_{2}$ PSG classes if the gauge symmetry is lowered from U(1) to $\mathbb{Z}_{2}$. We use Schwinger bosons rather than Abrikosov fermions Sanyal et al. (2019) to immediately obtain a bosonic excitation, the elementary Schwinger boson itself, that can condense at the phase transition out of the QSL.

As a result of our PSG analysis, we find $16$ different $\mathbb{Z}_{2}$ QSLs on the pyrochlore lattice. We use a standard mean-field description to study the [math]-flux QSLs, in which translation symmetry acts linearly (i.e., as in classical physics) on the Schwinger bosons. The PSG method also allows us to describe phase transitions from these QSLs to magnetically ordered phases. Condensing the Schwinger bosons, we identify $15$ different ordering patterns, and call them “paraphases”, since each of them actually unifies several distinct symmetry-breaking orders. We find that, generically, these orders are intertwined, necessarily appearing together at the phase transition out of the QSL, and that conventional spin orders are in many cases accompanied by inversion-breaking “hidden” orders.

The phase transitions corresponding to these $15$ paraphases fall into two dynamical classes of $z=1$ and $z=2$ quantum criticality, exhibiting critical modes with linear and quadratic dispersions, respectively. We uncover the mathematical structure discriminating between these two classes, related to Hamiltonian diagonalizability, and derive their effective field theories, along with their most important experimental signatures. In particular, we use mean-field theory to compute static and dynamic spin structure factors for each of the $15$ paraphases. Finally, by comparing the magnetic orders associated with each paraphase to those observed in experiments, we identify a set of likely QSL phases that might be relevant to real-world pyrochlore materials.

The rest of the paper is organized as follows. First, in Sec. II, we summarize our main results on the different QSL phases and the corresponding phase transitions out of them (“paraphases”). In Sec. III, we employ the PSG method, deriving the PSG classes, and constructing a mean-field theory for each PSG class. In Sec. IV, we analyze the mean-field theories of our QSL phases, describing phase transitions out of them, and establishing the two dynamical classes with critical exponents $z=1,2$. In Sec. V, we move on to the experimental signatures of our phase transitions, describing the heat capacity and the spin structure factors, and also introducing the concept of intertwined and hidden orders. Finally, in Sec. VI, we discuss our results and connect them to existing experimental data. Detailed derivations and lengthy formulas are given in the Appendices for reference.

II Main Results

From our PSG classification scheme, we find that there are $16$ different $\mathbb{Z}_{2}$ PSG classes of Schwinger bosons, corresponding to $16$ inequivalent $\mathbb{Z}_{2}$ QSL phases, on the pyrochlore lattice. Out of these $16$ different QSLs, there are eight [math]-flux QSLs and eight $\pi$-flux QSLs. For each QSL, we construct a general quadratic mean-field Hamiltonian for the Schwinger bosons containing all onsite, nearest-neighbor (NN), and next-nearest-neighbor (NNN) terms allowed by symmetry. However, for simplicity, we focus on the [math]-flux QSLs and restrict the mean-field Hamiltonian to onsite and NN terms. At such a NN level, two out of eight [math]-flux Hamiltonians have an enlarged U(1) gauge symmetry, and we thus concentrate on the remaining six [math]-flux Hamiltonians with $\mathbb{Z}_{2}$ gauge symmetry.

In each of the six corresponding $\mathbb{Z}_{2}$ QSL phases, the Schwinger bosons can be identified as elementary spinon excitations carrying fractionalized quantum numbers. If the chemical potential is tuned to its critical value, there is a phase transition driven by the condensation of these bosonic spinons. Depending on the particular patterns of spinon condensation, we describe $15$ different critical “paraphases” out of the six QSL phases. The most important characteristics of these paraphases, labeled by their parent QSL phases and the condensation momenta of the spinons, are tabulated in Table 1.

For each paraphase, the spinon spectrum is gapless at the critical point by construction. The effective field theory of the critical point is characterized by the low-energy spinon dispersion, $\omega\sim k^{z}$, in terms of the dynamical critical exponent, which is either $z=1$ or $z=2$. These two dynamical classes give rise to distinct sets of experimental signatures. For a start, the power-law exponent $x$ of the low-temperature heat capacity, $C_{V}\sim T^{x}$, is determined by the dynamical exponent $z$ and the dimensionality of the condensation manifold, i.e., if the spinons condense at points or along lines in the Brillouin zone (BZ). Also, the dynamical exponent gives rise to universal features in the static and dynamic spin structure factors, which appear on top of more detailed characteristics specific to given paraphases. In particular, when approaching zero energy, the spectral weight in the dynamic structure factor vanishes for $z=2$ but diverges for $z=1$; the divergence in the $z=1$ case is also observable as a nonanalytic behavior in the static structure factor.

To establish a connection between spinon condensation and the resulting magnetic orders, restricted to zero momentum for simplicity, we investigate the transformation rules of the possible order parameters under the point group $O_{\text{h}}$ of the pyrochlore lattice. For each paraphase, we determine which magnetic orders generically appear, concentrating in particular on the conventional spin orders seen in the experiments: the all-in-all-out, antiferromagnetic, ferromagnetic, and Palmer-Chalker orders. In doing so, we learn two important general lessons on magnetic orders obtained by spinon condensation. First, several distinct orders may be intertwined, i.e., they necessarily accompany each other, even though they are completely unrelated on the classical level. Second, the conventional spin orders may emerge together with more exotic inversion-breaking “hidden” orders.

III Projective symmetry group

III.1 Lattice symmetries

We first introduce the symmetries of the pyrochlore lattice; the convention and notation we establish here is used throughout the rest of the paper. The pyrochlore lattice consists of four FCC-type sublattices, which we label by $\mu=0,1,2,3$. To index the sites of the lattice, we use two coordinate systems: the global cartesian coordinates (GCCs) and the sublattice-indexed pyrochlore coordinates (SIPCs). The GCCs are the standard frame coordinates for the FCC cube of edge length $a=1$. The SIPCs are spanned by the lattice vectors $\hat{e}_{1}$, $\hat{e}_{2}$, and $\hat{e}_{3}$, which are expressed in GCCs as

[TABLE]

We define $\hat{\epsilon}_{i}=\frac{1}{2}\hat{e}_{i}$ ($i=1,2,3$) to be the displacement vectors of the $\mu=1,2,3$ sublattices from the $\mu=0$ sublattice, where we understand $\hat{\epsilon}_{0}=\hat{e}_{0}=0$. The relation between the SIPCs and the GCCs is then

[TABLE]

The space group of the pyrochlore lattice is the cubic space group $Fd\overline{3}m$ (No. 227), minimally generated by the translations $T_{1}$, $T_{2}$, and $T_{3}$ along the lattice vectors $\hat{e}_{1}$, $\hat{e}_{2}$, and $\hat{e}_{3}$, a sixfold rotoreflection $\overline{C}_{6}$ around the $[111]$ axis (i.e., around $\hat{e}_{1}+\hat{e}_{2}+\hat{e}_{3}$), and a non-symmorphic screw operation $S$, which is the composition of a twofold rotation around $\hat{e}_{3}$ and a translation by $\hat{\epsilon}_{3}$. These space-group generators transform the SIPCs according to

[TABLE]

Note that we can write the rotoreflection as $\overline{C}_{6}=I{C}_{3}$, where $I$ is an inversion with respect to the origin and ${C}_{3}$ is a threefold rotation around the $[111]$ axis. The generators $\{I,C_{3}\}$ are therefore equivalent to the generator $\overline{C}_{6}$; we choose a single generator $\overline{C}_{6}$ to reduce the number of generators and group relations.

The point group of the pyrochlore lattice, formally defined as the quotient group of the space group and the group of pure translations, is the cubic group $O_{\text{h}}$. This group is minimally generated by $\overline{C}_{6}$ and $S^{\prime}$, where $S^{\prime}$ is a twofold rotation around $\hat{e}_{3}$, distinguished from the space-group generator $S$ by the lack of a subsequent translation along $\hat{\epsilon}_{3}$. A detailed description of the point-group structure is given in Appendix A.

In addition to the pyrochlore space-group symmetries, time-reversal symmetry is also present in the pyrochlore materials. The corresponding time-reversal operation $\mathcal{T}$ commutes with all space-group operations and satisfies $\mathcal{T}^{2}=-1$ when acting on a half-integer spin state. The complete list of independent group relations defining the symmetry group is then

[TABLE]

The notation in Eq. (3) is understood as $i+3\equiv i$.

III.2 Projective symmetry group

In this subsection, we classify all possible $\mathbb{Z}_{2}$ quantum spin liquids that are compatible with the symmetries of the pyrochlore lattice. We first write the spins in terms of Schwinger boson bilinears as

[TABLE]

where b_{\vec{r}_{\mu}}=\left(\begin{array}[]{c}b_{\vec{r}_{\mu},\uparrow}\\ b_{\vec{r}_{\mu},\downarrow}\end{array}\right), and $\sigma^{x,y,z}$ are the Pauli matrices (also denoted by $\sigma^{1,2,3}$, respectively). Physically, the Schwinger bosons $b_{\vec{r}_{\mu}}$ describe the deconfined spinon excitations of the quantum spin liquid and, on the mean-field level, they are governed by a quadratic Hamiltonian, commonly known as the mean-field ansatz.

It is important to emphasize that the transformation in Eq. (4) is not faithful as it enlarges the local Hilbert space at each site $\vec{r}_{\mu}$. Consequently, there is a local gauge redundancy for the Schwinger bosons. Indeed, any site-dependent U(1) phase transformation

[TABLE]

leaves the spins $\hat{S}^{\alpha}_{\vec{r}_{\mu}}$ invariant. The physical Hilbert space can in principle be retained by enforcing the constraint

[TABLE]

at each site $\vec{r}_{\mu}$ of the lattice.

Under a space-group operation $\mathcal{O}$, the spins transform as $\mathcal{O}\colon\hat{S}^{\alpha}_{\vec{r}_{\mu}}\rightarrow U_{\mathcal{O}}\hat{S}^{\alpha}_{\mathcal{O}(\vec{r}_{\mu})}U^{\dagger}_{\mathcal{O}}=\frac{1}{2}b^{\dagger}_{\mathcal{O}(\vec{r}_{\mu})}U_{\mathcal{O}}\sigma^{\alpha}U_{\mathcal{O}}^{\dagger}b_{\mathcal{O}(\vec{r}_{\mu})}$, where $U_{\mathcal{O}}$ is the SU(2) rotation matrix associated with the operation $\mathcal{O}$. We therefore naïvely expect that the spinons transform as

[TABLE]

However, due to the U(1) gauge redundancy, any operation $\mathcal{O}$ is generally accompanied by a site-dependent U(1) phase transformation

[TABLE]

and the spinons thus actually transform as

[TABLE]

where the symbol “$\circ$” indicates that the gauge-enriched operation $\widetilde{\mathcal{O}}$ is a composition of the pure symmetry operation $\mathcal{O}$ and the gauge transformation $G_{\mathcal{O}}$.

Under a time reversal $\mathcal{T}$ of the system, the spins transform as $\mathcal{T}\colon\hat{S}^{\alpha}_{\vec{r}_{\mu}}\rightarrow\mathcal{K}^{\dagger}U_{\mathcal{T}}\hat{S}^{\alpha}_{\vec{r}_{\mu}}U^{\dagger}_{\mathcal{T}}\mathcal{K}$, where $U_{\mathcal{T}}=i\sigma^{2}$, while $\mathcal{K}=\mathcal{K}^{\dagger}=\mathcal{K}^{-1}$ applies complex conjugation to everything on its right. Once again, combining the naïve transformation rule for the spinons,

[TABLE]

and the accompanying U(1) phase transformation,

[TABLE]

the spinons are found to transform as

[TABLE]

Note that $[\mathcal{K},U_{\mathcal{T}}]=0$ because $U_{\mathcal{T}}$ is real.

For a quantum spin liquid, the gauge-enriched operations $\widetilde{\mathcal{O}}$ and $\widetilde{\mathcal{T}}$ generate the symmetry group of the mean-field ansatz, commonly known as the projective symmetry group (PSG). To enumerate all quantum spin liquids, we need to find all distinct PSG solutions, i.e., all gauge-inequivalent solutions for the gauge transformations $G_{\mathcal{O}}$ and $G_{\mathcal{T}}$ that are consistent with the symmetry group of the lattice, including space-group symmetries and time-reversal symmetry. In particular, for each group relation [see Eq. (3)] taking the general form of

[TABLE]

we consider the gauge-enriched group relation

[TABLE]

where $\mathcal{G}$ is a pure gauge transformation, thus corresponding to the identity operation for the spins. Being an element of the PSG by definition, $\mathcal{G}$ is also an element of the invariant gauge group (IGG), the group of all gauge transformations that leave the mean-field ansatz invariant. In most cases, such gauge transformations are exclusively “global” (i.e., site independent), and the IGG is thus a subgroup of U(1), typically $\mathbb{Z}_{2}$ or U(1), corresponding to $\mathbb{Z}_{2}$ and U(1) spin liquids, respectively. Since we are interested in classifying $\mathbb{Z}_{2}$ spin liquids, we consider IGG $=\mathbb{Z}_{2}$ in the following. The only two elements of the IGG are then $\mathcal{G}=e^{in\pi}$ with $n=\{0,1\}$.

For any group relation in terms of exclusively space-group operations, taking the form of Eq. (13), the gauge-enriched group relation in Eq. (14) can be rewritten as

[TABLE]

Using the general conjugation rule

[TABLE]

following directly from Eqs. (7) and (8), this group relation then becomes a pure phase equation:

[TABLE]

For group relations involving time reversal, special care must be taken due to the presence of the complex conjugation $\mathcal{K}$. Using the modified conjugation rule

[TABLE]

the last group relation in Eq. (3) translates into

[TABLE]

while the penultimate group relation $\mathcal{T}^{2}=-1$ gives rise to a trivial equation due to the cancellation between the phase factors $e^{i\phi_{\mathcal{T}}(\vec{r}_{\mu})}$ and $e^{-i\phi_{\mathcal{T}}(\vec{r}_{\mu})}$.

The PSG classification is obtained by listing all group relations and finding all solutions of the corresponding phase equations [see Eqs. (III.2) and (19)] for the $\mathbb{Z}_{2}$ parameters $n$ as well as the phases $\phi_{\mathcal{O}}(\vec{r}_{\mu})$ and $\phi_{\mathcal{T}}(\vec{r}_{\mu})$. We emphasize that distinct solutions, describing distinct spin liquids, must be gauge inequivalent. Indeed, by means of a general gauge transformation $G$ [see Eq. (5)], the gauge-enriched group relations in Eq. (14) can be rewritten as

[TABLE]

transforming the phases $\phi_{\mathcal{O}_{i}}(\vec{r}_{\mu})$ according to

[TABLE]

and thus indicating that two seemingly distinct solutions for the phases might in fact be equivalent.

The detailed solution of the PSG equations is presented in Appendix B. The PSG results for the phases are

[TABLE]

where $n_{1}$, $n_{\overline{C}_{6}S}$, $n_{ST_{1}}$, and $n_{\overline{C}_{6}}$ are four $\mathbb{Z}_{2}$ parameters, each being either [math] or $1$. Therefore, we find that there are $16$ gauge-inequivalent $\mathbb{Z}_{2}$ PSG classes, corresponding to distinct $\mathbb{Z}_{2}$ quantum spin liquids, which we label by the notation $n_{1}\pi$-$(n_{\overline{C}_{6}S}\,n_{ST_{1}}\,n_{\overline{C}_{6}})$. The four $\mathbb{Z}_{2}$ parameters have concrete interpretations:

• •

The parameter $n_{1}$ comes from the three PSG equations corresponding to $T_{i}T_{i+1}T_{i}^{-1}T_{i+1}^{-1}=1$, which are required by the PSG to share the same $\mathbb{Z}_{2}$ parameter. Physically, it quantifies the Aharonov-Bohm (AB) phase a spinon accumulates while moving on the closed edge of a plaquette, which is traversed by such a sequence of translations. In the case of $n_{1}=1$ ($n_{1}=0$), the AB phase is $\pi$ ([math]), corresponding to a $\pi$-flux ([math]-flux) spin liquid.

• •

The parameter $n_{\overline{C}_{6}}$ comes from the PSG equation corresponding to $\overline{C}_{6}^{6}=1$. Physically, it describes the AB phase a spinon accumulates after completing six subsequent sixfold rotoreflections. Together with $n_{ST_{1}}$, it determines whether or not the sixfold rotoreflection $\overline{C}_{6}$ acts projectively.

• •

The parameter $n_{ST_{1}}$ comes from the PSG equation corresponding to $ST_{1}S^{-1}T_{3}^{-1}T_{1}=1$. Physically, it describes the AB phase a spinon accumulates after completing the operation sequence $ST_{1}S^{-1}T_{3}^{-1}T_{1}$. Together with $n_{1}$ and $n_{\overline{C}_{6}S}$ it determines whether or not the screw operation $S$ acts projectively.

• •

The parameter $n_{\overline{C}_{6}S}$ comes from the PSG equation corresponding to $(\overline{C}_{6}S)^{4}=1$. Physically, it describes the AB phase a spinon accumulates after completing the operation sequence $(\overline{C}_{6}S)^{4}$.

III.3 Construction of mean-field ansätze

We are now in the position to construct the mean-field ansatz for each PSG class. The most general mean-field ansatz for bosonic spinons can be written as

[TABLE]

where $u^{h}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$ and $u^{p}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$ are $2\times 2$ matrices acting on spin space, and the labels “h” and “p” indicate hopping and pairing terms, respectively.

The PSG operators $\widetilde{\mathcal{O}}$ and $\widetilde{\mathcal{T}}$ are the symmetry operators of the Hamiltonian $H$, meaning $\widetilde{\mathcal{O}}\colon H\rightarrow H$ and $\widetilde{\mathcal{T}}\colon H\rightarrow H$. Since the spinons transforms under $\widetilde{\mathcal{O}}$ and $\widetilde{\mathcal{T}}$ according to Eqs. (9) and (12), the matrices $u^{h}$ and $u^{p}$ must transform as

[TABLE]

for space-group elements $\mathcal{O}\in\{T_{1},T_{2},T_{3},\overline{C}_{6},S\}$ and as

[TABLE]

for time reversal $\mathcal{T}$. The respective SU(2) matrices are

[TABLE]

where $\sigma^{0}=1_{2\times 2}$ is the identity matrix. Suppressing the site indices for simplicity, we parameterize the matrices $u^{h}$ and $u^{p}$ in the general forms

[TABLE]

where $a,b,c,d,a^{\prime},b^{\prime},c^{\prime},d^{\prime}$ are all complex. The additional factor $i\sigma^{2}$ appearing in $u^{p}$ ensures that $(a,b,c,d)$ and $(a^{\prime},b^{\prime},c^{\prime},d^{\prime})$ transform in the same way under the respective unitary conjugations $u^{h}\rightarrow Uu^{h}U^{{\dagger}}$ and $u^{p}\rightarrow Uu^{p}U^{T}$ for any $U\in$ SU(2). In both cases, the singlet parameters $a$ and $a^{\prime}$ transform as scalars, while the triplet parameters $\vec{b}=(b,c,d)$ and $\vec{b}^{\prime}=(b^{\prime},c^{\prime},d^{\prime})$ transform as SO(3) vectors. Indeed, any SU(2) rotation leaves the singlet parameters invariant and performs the corresponding SO(3) rotation on the triplet vectors: $\vec{b}\rightarrow\mathcal{R}\vec{b}$ and $\vec{b}^{\prime}\rightarrow\mathcal{R}\vec{b}^{\prime}$. For the generators $\overline{C}_{6}$ and $S$, these SO(3) rotations are

[TABLE]

while the translations $T_{1,2,3}$ correspond to trivial SO(3) rotations: $\mathcal{R}^{T_{1,2,3}}=1_{3\times 3}$.

To reduce the number of parameters in the mean-field ansatz, we first consider the effect of time reversal. Substituting Eq. (27) into Eq. (25), and taking $G_{\mathcal{T}}(\vec{r}_{\mu})=1$ from Eq. (22d), we obtain $(a,b,c,d)=(a^{*},b^{*},c^{*},d^{*})$ as well as $(a^{\prime},b^{\prime},c^{\prime},d^{\prime})=(a^{\prime*},b^{\prime*},c^{\prime*},d^{\prime*})$ and deduce that all $8$ parameters of $u^{h}$ and $u^{p}$ are real.

Turning to space-group symmetries and using Eq. (24), we can then establish relations between the respective parameters of $u^{h}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$ and $u^{p}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$ that correspond to different bonds $\langle\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}\rangle$ of the lattice. In fact, the entire mean-field ansatz in Eq. (23) can be constructed up to next-nearest-neighbor level by specifying the $8$ real parameters for each of the following three representative bonds:

• •

onsite “bond” $\vec{0}_{0}\rightarrow\vec{0}_{0}$:

[TABLE]

• •

nearest-neighbor (NN) bond $\vec{0}_{0}\rightarrow\vec{0}_{1}$:

[TABLE]

• •

next-nearest-neighbor (NNN) bond $\vec{0}_{1}\rightarrow\vec{0}_{2}-\hat{e}_{2}$:

[TABLE]

III.4 Nontrivial parameter constraints

When constructing the entire mean-field ansatz from the representative bonds in Eqs. (29)–(31), the significance of using Eq. (24) is twofold. On the one hand, most space-group elements map the representative bonds onto different bonds, thereby determining the matrices $u^{h}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$ and $u^{p}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$ for all symmetry-related bonds. On the other hand, some space-group elements map the representative bonds onto themselves, thereby leading to nontrivial constraints on the original $24$ parameters.

For simplicity, we first concentrate on the [math]-flux PSG classes. Since translation is trivial [see Eqs. (22a)–(22c)], we can restrict our attention to a single unit cell, within which bonds are mapped onto each other by elements of the point group. Since the point group $O_{\text{h}}$ consists of $48$ elements, and there are $4$ onsite, $12$ NN, and $24$ NNN bonds within a single unit cell, which can be viewed as three orbits in the point group, the orbit-stabilizer theorem implies that the onsite, NN, and NNN representative bonds are mapped onto themselves by $12$, $4$, and $2$ point-group elements, respectively. When a bond is mapped onto itself by such a point-group element, nontrivial constraints are obtained on the parameters by comparing the new and the old expressions for $u^{h}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$ and $u^{p}_{\vec{r}_{\mu},\vec{r}^{\prime}_{\nu}}$. These constraints are solved in Appendix D.

In Table 2, we present the nonzero parameters of the mean-field ansatz for each of the eight [math]-flux and each of the eight $\pi$-flux PSG classes up to NNN level, along with any constraints between the parameters. From these nonzero parameters, the entire mean-field ansatz can be constructed via Eq. (24). Note that some of the mean-field ansätze in Table 2 have an enlarged U(1) gauge symmetry at the NN level which only breaks down to $\mathbb{Z}_{2}$ when nonzero NNN terms are included.

IV Analysis of the mean-field ansätze

The previous section explains how the method of PSG can be used to obtain classes of 0-flux and $\pi$-flux mean-field ansätze, which describe distinct phases of $\mathbb{Z}_{2}$ quantum spin liquids on the mean-field level. In this section, we focus on the 0-flux mean-field ansätze and study their physical properties in great detail. Since our main goal is to explore the relationship between spin liquids and magnetic orders adjacent to them, we primarily concentrate on the critical field theories and the condensation patterns (i.e., the resulting magnetic orders).

In each mean-field ansatz, we neglect the NNN terms for simplicity, restricting our attention to onsite and NN terms. Since we are interested in $\mathbb{Z}_{2}$ spin liquids, and two out of eight 0-flux mean-field ansätze have U(1) gauge symmetry at the NN level, we only consider the remaining six mean-field ansätze in the rest of the paper.

IV.1 Symmetry properties

The PSG method is rooted in symmetry analysis, and it is important to understand how the PSG governs the symmetry of the mean-field Hamiltonians. By means of a Fourier transformation, a general mean-field Hamiltonian [see Eq. (23)] can be written in momentum space as

[TABLE]

where $B_{\vec{k}}=\left(b^{\phantom{{\dagger}}}_{\vec{k},0},b^{\phantom{{\dagger}}}_{\vec{k},1},b^{\phantom{{\dagger}}}_{\vec{k},2},b^{\phantom{{\dagger}}}_{\vec{k},3},b^{\dagger}_{-\vec{k},0},b^{\dagger}_{-\vec{k},1},b^{\dagger}_{-\vec{k},2},b^{\dagger}_{-\vec{k},3}\right)^{T}$ is a $16$-component vector of operators. The matrix $\mathcal{H}(\vec{k})$ has the standard Bogoliubov form

[TABLE]

where $U_{h}(\vec{k})=U^{\dagger}_{h}(\vec{k})$ and $U_{p}(\vec{k})=U_{p}^{T}(-\vec{k})$, corresponding to hopping and pairing terms, respectively.

The Hamiltonian matrix $\mathcal{H}(\vec{k})$ combines momenta $\pm\vec{k}$ and thus assigns a full set of physical degrees of freedom to only half of the BZ. This redundancy in the description leads to an effective charge-conjugation “symmetry”, corresponding to the matrix-level constraint

[TABLE]

where we define $U_{\mathcal{C}}=\sigma^{1}\otimes 1_{8\times 8}$. The anti-unitary charge-conjugation operator is then given by $U_{\mathcal{C}}\mathcal{K}$, where $\mathcal{K}$ denotes complex conjugation.

Considering physical symmetries, time reversal $\mathcal{T}$ gives rise to an analogous matrix-level constraint

[TABLE]

where we define $U_{\mathcal{T}}=1_{8\times 8}\otimes(i\sigma^{2})$. Correspondingly, the anti-unitary time-reversal operator is $U_{\mathcal{T}}\mathcal{K}$. Note that time reversal acts non-projectively in all PSG classes because we use gauge freedom to fix $\phi_{\mathcal{T}}(\vec{r}_{\mu})=0$.

In contrast, inversion $I=\overline{C}_{6}^{3}$ acts projectively on the spinons and generates the matrix-level constraint

[TABLE]

where $U_{I}(\vec{k})=\left(\sigma^{3}\right)^{n_{\overline{C}_{6}}}\otimes\left(U_{J}\cdot I^{2}(\vec{k})\right)\otimes\sigma^{0}$, in terms of the $4\times 4$ diagonal form-factor matrix

[TABLE]

and the diagonal matrix $U_{J}=\text{Diag}\left((-1)^{n_{ST_{1}}},1,1,1\right)$.

The symmetries $\mathcal{C}$, $\mathcal{T}$, and $I$ result in important general spectral features. First, the symmetry $I\circ\mathcal{T}$ guarantees that each energy level is doubly degenerate, according to Kramers theorem. Second, the symmetry $I\circ\mathcal{C}$ leads to an additional double degeneracy for any non-zero-energy level, which is connected to the redundant description in Eqs. (32) and (33). The two symmetries together thus result in a generic four-fold degeneracy at each energy level $E>0$ shared by momenta $\pm\vec{k}$. Note that the degeneracy may be smaller or larger at special time-reversal-invariant momenta ($\vec{k}=-\vec{k}$) because there are half as many physical degrees of freedom but, on the other hand, pure point-group symmetries (e.g., inversion) may lead to additional degeneracy.

The degeneracy of zero-energy levels is more subtle as it may be affected by the diagonalizability of the Hamiltonian matrix $\mathcal{H}(\vec{k})$. Since the low-energy physics is the main focus of our study, this issue will be addressed in a separate section (see Sec. IV.4).

IV.2 Condensation domains: a “phase diagram” for paraphases

The use of bosonic mean-field Hamiltonians, obtained from the spinon decomposition in Eq. (4), facilitates the study of phase transitions between spin liquids and magnetically ordered phases. Indeed, by lowering the chemical potential $\mu$, there is a critical chemical potential $\mu_{c}$ at which the bosonic spinons undergo Bose-Einstein condensation at some critical momenta $\vec{k}_{c}$ and the system thus develops magnetic order.

For the mean-field Hamiltonian in each PSG class, the critical chemical potential $\mu_{c}$ is a function of the mean-field parameters (see Table 11 in Appendix E for details). While the value of $\mu_{c}$ changes continuously with the mean-field parameters, and this variation of $\mu_{c}$ is thus locally analytic, it globally separates into domains across which the variation of $\mu_{c}$ is non-analytic. These domains of analyticity of $\mu_{c}$ are reminiscent of the domains of analyticity of the free energy, which define phases in thermodynamics. However, the analogy is not perfect as each such domain may give rise to several true phases on crossing the phase transition into magnetic order (i.e., when taking $\mu<\mu_{c}$). We therefore coin the word paraphases to describe the distinct domains of analyticity of $\mu_{c}$. Restated, each paraphase is a connected region of phase space in which the unstable manifold of condensation modes varies smoothly.

Following this logic, the six mean-field Hamiltonians are further divided into $15$ paraphases (see Fig. 1). The analytical expressions for the paraphase boundaries are given in Table 3, while the distinct critical momenta $\vec{k}_{c}=\Gamma,\mathrm{L},\Lambda,\mathrm{X},\mathrm{W}$ characterizing the various paraphases are explained in Table 4. Finally, the distinct expressions for the critical chemical potentials $\mu_{c}$ in the $15$ paraphases are specified in Appendix E.

Note that, for the PSG classes with $n_{ST_{1}}=1$, the PSG result for the screw operation $S$ depends on the spatial coordinates, and it is convenient to shift the entire BZ by the translation $\vec{k}\rightarrow\vec{k}-\pi(1,1,1)$. Such a shift of the BZ can be thought of as a gauge transformation of the spinons, which does not modify any physical quantities on the spin level. This shift is assumed throughout the paper and is already taken into account when specifying the condensation momenta in Fig. 1.

Note also that the region $\Lambda$ supports a one-dimensional manifold of condensation momenta. Since the only physical symmetries are discrete space-group and time-reversal symmetries, this ground-state continuum must be accidental, i.e., the result of restricting the mean-field Hamiltonians to NN level. Indeed, when including infinitesimal NNN parameters, we see that the condensation regions are reduced from $\Lambda$ to either $\Gamma$ or $\mathrm{L}$.

IV.3 Critical spectra

The critical spectra of the $15$ paraphases, corresponding to $\mu=\mu_{c}$ in each case, are shown in Fig. 4, along with the associated dynamical spin structure factors, obtained on the mean-field level. Generically, each of these spectra consists of four bands, which is consistent with the fourfold degeneracy of each band. While certain spectra have distinguishing features, not all paraphases can be fully distinguished by their spectra, as some spectral characteristics are shared by multiple paraphases. Among other features, several spectra show a quasi-mirror-reflection symmetry (in terms of energy) between two bands, which accounts for certain high-energy features in the dynamic spin structure factor (see Sec. V.2).

Most importantly, however, the critical paraphases can be divided into two classes, characterized by linear and quadratic dispersions at low energies. In terms of the dynamical critical exponent $z$, defined by $\omega\sim|k-k_{c}|^{z}$ and specified for each paraphase in Table 1, these two classes are labeled by $z=1$ and $z=2$, respectively. As we later show, paraphases with $z=1$ and $z=2$ correspond to different critical field theories, which determine the critical exponents of various physical observables, such as the heat capacity and the magnetic susceptibility, and thus lead to distinct experimental signatures.

IV.4 Hamiltonian diagonalizability

From a technical point of view, the distinction between $z=1$ and $z=2$ theories becomes evident when we try to diagonalize the Hamiltonian matrix in Eq. (33). In general, we seek a change of basis for the bosonic operators,

[TABLE]

such that the Hamiltonian in Eq. (32) is of the form

[TABLE]

where $\Lambda(\vec{k})=V^{\dagger}(\vec{k})\mathcal{H}(\vec{k})V(\vec{k})$ is a $16\times 16$ diagonal matrix, and $V(\vec{k})\in\text{SU}(8,8)$ satisfies

[TABLE]

ensuring that this change of basis is a canonical transformation. The matrices $\Lambda(\vec{k})$ and $V(\vec{k})$ can be found by solving the generalized eigenvalue problem

[TABLE]

where the eigenvalues $\lambda_{\vec{k},i}$ give the diagonal elements of the matrix $J\Lambda(\vec{k})$ and the eigenvectors $\bm{a}_{\vec{k},i}$ form the columns of the matrix $V(\vec{k})$. However, since $J\mathcal{H}(\vec{k})$ is not necessarily Hermitian (or even normal), it is not guaranteed that such a matrix $V(\vec{k})$ actually exists.

In particular, it may happen at the critical chemical potential $\mu=\mu_{c}$ that there are not enough independent eigenvectors for the zero eigenvalues $\lambda_{\vec{k}_{c}}=0$. Physically, this scenario means that we cannot diagonalize our critical Hamiltonian by a canonical transformation of bosonic creation and annihilation operators, and instead we must decompose our complex operators $B_{\vec{k}}$ according to

[TABLE]

where $\hat{X}_{\vec{k}}$ and $\hat{P}_{\vec{k}}$ are $16$-dimensional vectors of real operators, analogous to the position and momentum operators in first quantization. In terms of these new operators, the analog for the change of basis in Eq. (38) is

[TABLE]

where the $32\times 32$ matrix $W(\vec{k})$ satisfies

[TABLE]

Using this canonical change of basis, the Hamiltonian can then be brought to the diagonal form

[TABLE]

where $\hat{y}_{\vec{k},i}$ and $\hat{q}_{\vec{k},i}$ are the components of $\hat{Y}_{\vec{k}}$ and $\hat{Q}_{\vec{k}}$, respectively, and the new eigenvalues are related to the original ones by

[TABLE]

Importantly, however, unlike the original method of diagonalization, $B_{\vec{k}}\rightarrow\widetilde{B}_{\vec{k}}$, which may fail if $J\mathcal{H}(\vec{k})$ is a defective matrix, the alternative method of diagonalization, $(\hat{X}_{\vec{k}},\hat{P}_{\vec{k}})\rightarrow(\hat{Y}_{\vec{k}},\hat{Q}_{\vec{k}})$, always works.

For any zero mode $i$ at a critical momentum $\vec{k}_{c}$, we have $\alpha_{\vec{k}_{c},i}\beta_{\vec{k}_{c},i}=0$ from Eq. (46). The diagonalizability of the critical Hamiltonian $\mathcal{H}(\vec{k}_{c})$ is then determined by the following simple criterion:

• •

if $\alpha_{\vec{k}_{c},i}=\beta_{\vec{k}_{c},i}=0$, the Hamiltonian can be diagonalized in the original basis of creation and annihilation operators;

• •

otherwise, either $\alpha_{\vec{k}_{c},i}=0$, $\beta_{\vec{k}_{c},i}\neq 0$ or $\alpha_{\vec{k}_{c},i}\neq 0$, $\beta_{\vec{k}_{c},i}=0$; the Hamiltonian is not diagonalizable in any creation-annihilation-operator basis, meaning that the SU($8$,$8$) transformation is singular.

To understand how these two scenarios for the diagonalizability lead to theories of $z=2$ and $z=1$ types, respectively, we now switch to the language of path integrals and consider the critical low-energy actions.

IV.5 Effective low-energy theories

Our phase transitions from spin liquids to magnetic orders, driven by a change in the chemical potential $\mu$, are prototypes of quantum critical points (QCPs). At such a QCP, one can write down an effective theory in terms of the low-energy degrees of freedom. We assume a single condensing eigenmode obtained from the Hamiltonian $\mathcal{H}(\vec{k}_{c})$, denoted by $\widetilde{b}_{\vec{k}_{c}}$. Including spatial fluctuations, we promote this eigenmode to a field $\phi(\tau,\vec{x})$ and consider the imaginary-time action $S=\int\mathcal{L}\,d^{3}xd\tau$. If the Hamiltonian is diagonalizable, the critical Lagrangian becomes

[TABLE]

describing a massless field $\phi$ at the QCP. The corresponding action is invariant under the rescaling

[TABLE]

from which we can immediately deduce that the dynamical critical exponent is $z=2$.

However, the mass of $\phi$ should be generally considered as a tensor of real fields $\chi$ and $\pi$, which are the real and imaginary components of $\phi$, such that

[TABLE]

In the Hamiltonian language, these two components correspond to the “position” and “momentum” operators in Eq. (42). Consequently, if the Hamiltonian is not diagonalizable, only one of these components is massless at the QCP. Assuming without loss of generality that $\chi$ is massive and $\pi$ is massless, the critical Lagrangian becomes

[TABLE]

By integrating out the massive field $\chi$ and rescaling the massless field as $\pi\rightarrow r\pi$, we finally obtain

[TABLE]

This effective action is invariant under the rescaling

[TABLE]

from which we can immediately deduce that the dynamical critical exponent is $z=1$.

These two distinct QCPs, characterized by critical exponents $z=2$ and $z=1$, respectively, are reminiscent of the QCPs governing phase transitions from quantum paramagnets to XY antiferromagnets Zapf et al. (2014). If such a transition is induced by an external magnetic field, the QCP is described by the $z=2$ critical theory in Eq. (47), while if the transition is induced by pressure and is thus time-reversal symmetric, the QCP is described by the $z=1$ critical theory in Eq. (51).

IV.6 Spin condensation: order patterns

We are now ready to describe the spin orders obtained by condensing the spinons in each of the $15$ critical paraphases. When the chemical potential $\mu$ reaches its critical value $\mu_{c}$, certain spinons $\widetilde{b}_{\vec{k}_{c}}$ at critical momenta $\vec{k}_{c}$ condense and thereby acquire macroscopic occupation numbers $\langle\widetilde{b}_{\vec{k}_{c}}\rangle$. We can then use these $\langle\widetilde{b}_{\vec{k}_{c}}\rangle$ as order parameters and detect spin orders by looking at order parameter bilinears, which, according to the spinon decomposition in Eq. (4), recover spin expectation values.

So far, several types of orders have been successfully identified in pyrochlore materials, most of which do not break translation symmetry. These zero-momentum orders correspond to representations of the point group $O_{\text{h}}$ and can thus be analyzed by the standard representation theory of groups. We will defer such an effort to the next section. In this subsection, we select several paraphases with definite ordering signatures and explicitly calculate the spin expectation values via condensing spinons. This way, we capture a limited set of orders, which correspond to irreducible representations (irreps) of the tetrahedral group $T_{\text{d}}$ (see Appendix C), and show that all such orders can be obtained from at least one of the six $\mathbb{Z}_{2}$ spin liquids. We mainly restrict ourselves to NN terms in the spinon Hamiltonian but include NNN terms whenever necessary.

One must bear in mind that the simplified irrep analysis on these explicit spin-condensation orders may be incomplete. For example, we will find from such an analysis that pure all-in-all-out order may be obtained in the paraphase 0-$(100)\Gamma$, while a full representation-theory analysis in Sec. V.3 leads to Table 8, which indicates that all-in-all-out order is always intertwined with some hidden orders (i.e., it can never appear alone). Still, the naïve spin-condensation analysis in this subsection is a good starting point to build some insight into how the six spin liquids are physically distinct from each other.

IV.6.1 All-in-all-out order

We consider the paraphase 0-$(100)\Gamma$, but also remark that the paraphases 0-(101)$\Gamma$ and 0-(111)$\Gamma$ give similar results. At the critical chemical potential $\mu_{c}$, the zero-energy subspace is twofold degenerate. The zero-energy eigenvectors are obtained from Eq. (41) and are given by the time-reversal partners $\bm{a}$ and $U_{\mathcal{T}}\bm{a}^{*}$. After condensing these two modes, the corresponding operators $\widetilde{b}_{1,2}$ acquire macroscopic occupation numbers $\langle\widetilde{b}_{i}\rangle=r_{i}e^{i\phi_{i}}$, with $i=1,2$, implying $\langle B_{\vec{k}_{c}}\rangle=\sum_{i=1,2}\bm{a}_{i}r_{i}e^{i\phi_{i}}$ at the critical momentum $\vec{k}_{c}=\Gamma$. In terms of the $12$-component vector $\mathbf{S}=(\vec{S}_{0},\vec{S}_{1},\vec{S}_{2},\vec{S}_{3})$ of the spin components on the four sublattices, we have, up to a global coefficient,

[TABLE]

where $r=(r_{1}^{2}-r_{2}^{2})/(2r_{1}r_{2})$, and $\mathbf{S}^{r,c,s}$ are three equimodular and mutually orthogonal vectors (see Appendix F for detail). Using the basis for the irreducible representations of $T_{\text{d}}$ (see Appendix C), it can be shown that this paraphase generically supports two orders: the all-in-all-out order and the AFM order. One can obtain pure all-in-all-out order [see Fig. 2 for illustration] by setting particular values for the condensation parameters $r_{1,2}$ and $\phi_{1,2}$.

IV.6.2 XY antiferromagnetic order

The paraphase 0-(110)$\Gamma$ has a non-diagonalizable critical Hamiltonian, because $\mathcal{H}$ has four zero-energy eigenvalues, but the nullspace of $J\mathcal{H}$ is only two dimensional, spanned by the time-reversal partners $\bm{a}$ and $U_{\mathcal{T}}\bm{a}^{*}$. We therefore switch to the position-momentum representation $(\hat{x},\hat{p})$, according to Eq. (42). The critical Hamiltonian is then diagonalized by a basis change $(\hat{x},\hat{p})\rightarrow(\hat{y},\hat{q})$ and takes the low-energy form

[TABLE]

which contains two gapless modes $\hat{y}_{1}$ and $\hat{y}_{2}$. To minimize the energy, we must have $\langle\hat{q}_{i}\rangle=0$ and, due to the uncertainty principle, $\hat{y}_{i}$ must fluctuate maximally. In terms of $y_{i}=\langle\hat{y}_{i}\rangle$, we then find $\langle B_{\vec{k}_{c}}\rangle=\bm{b}y_{1}+U_{\mathcal{T}}\bm{b}y_{2}$ at $\vec{k}_{c}=\Gamma$ for some vector $\bm{b}$ determined by $\bm{a}$, and the final result for spin configuration becomes

[TABLE]

where $\cos\theta=y_{1}^{2}-y_{2}^{2}$ and $\sin\theta=y_{1}y_{2}/2$. This spin configuration, shown in Fig. 2, corresponds to the “XY” order of the irrep $E$ obtained in Eq. (39) of Ref. Yan et al. (2017) after a redefinition $\theta\rightarrow\theta-\frac{\pi}{12}$.

IV.6.3 Ferromagnetic order: collinear and non-collinear

For the paraphase 0-(001)$\Gamma$, all pairing terms vanish at the $\Gamma$ point at the NN level and, solving the hopping part at $\mu_{c}$, we find that the zero-energy subspace is spanned by the time-reversal partners $\bm{a}$ and $U_{\mathcal{T}}\bm{a}^{*}$. There are two cases depending on the expression for $\mu$ in terms of the mean field parameters. When $\mu=-6a$, all four spins point in the same direction, which is the collinear FM order. When $\mu=2a-8c$, the spin vector $\mathbf{S}$ follows Eq. (53), where $\mathbf{S}^{r,c,s}$ are three equimodular and mutually orthogonal vectors (see Appendix F for detail). A typical spin configuration of such ferrimagnetic nature is shown in Fig. 2.

IV.6.4 Palmer-Chalker order

The paraphase 0-(101)$\mathrm{X}$ has a non-diagonalizable critical Hamiltonian, because $\mathcal{H}$ has eight zero-energy eigenvalues, but the nullspace of $J\mathcal{H}$ is only four dimensional. Switching to the $(\hat{x},\hat{p})$ representation and diagonalizing the Hamiltonian via the basis change $(\hat{x},\hat{p})\rightarrow(\hat{y},\hat{q})$, we find that there are four gapless modes $\hat{y}_{1,2,3,4}$ for each of the three critical momenta $\vec{k}_{c}=\mathrm{X}^{1,2,3}$. The expression for $\langle B_{\vec{k}_{c}}\rangle$ thus contains $12$ real parameters: the expectation values of the maximally fluctuating modes $\hat{y}_{1,2,3,4}$ at each critical momentum. Although most choices of these condensation parameters give an order with an enlarged unit cell, some special cases respect translation symmetry. For instance, if condensation is restricted to $X^{1}$, the spin configuration, shown in Fig. 2, corresponds to a Palmer-Chalker order, transforming under the irrep $T_{2}$.

IV.6.5 Non-uniform spinon condensations and partial orders

The spin expectation values $\vec{S}_{0,1,2,3}$ for the paraphase 0-(010)$\Gamma$ have different amplitudes on different sublattices, invalidating the irrep analysis that presupposed classically ordered states of fixed-length spins. There is no a priori reason to rule out such a non-uniform spin-amplitude state. It does, however, correspond to a more “exotic” ordered phase in which the spin is more ordered on some sublattices than others. This type of partially ordered state has been proposed in the material Gd2Ti2O7 Champion et al. (2001); Javanparast et al. (2015a) and in various theoretical models.

IV.6.6 Spinon line orders

The line orders $\Lambda$ appearing in classes 0-(001), 0-(010), 0-(100), and 0-(110) have accidental degeneracies, higher than demanded by the lattice symmetry. This extra degeneracy is an artifact of the restriction to NN ansätze, and should reduce to discrete condensation momenta in the presence of further-neighbor terms. If we include infinitesimal NNN terms to the mean-field ansatz, using Table 2, we indeed see that line condensation along $\Lambda$ shrinks to point condensation at either $\Gamma$ or $\mathrm{L}$. However, if we increase these NNN terms, the condensation points are shifted away from these high-symmetry points.

Due to the large NNN parameter space, we were unable to exhaustively study the effect of NNN terms on the NN mean-field ansatz. However, for some purposes, the NN level ansätze may be adequate. For example, as we explore in the next section, the line minima contribute to substantial low-energy continua in the dynamical spin structure factor. This feature should persist at intermediate energies when small NNN terms are included.

IV.6.7 Multi-spinon condensation orders

Spinon condensation at multiple critical momenta, in the paraphases 0-(001)$\mathrm{L}$, 0-(101)$\mathrm{X}/\mathrm{W}$, and 0-(111)$\mathrm{X}/\mathrm{W}$, allows for richer physics and is often accompanied by an enlargement of the unit cell. As an example, we look at the paraphase 0-(001)$\mathrm{L}$: there are two independent zero-energy modes at each critical momentum $\mathrm{L}$, and the four inequivalent $\mathrm{L}$ momenta thus give rise to an eight-dimensional zero-energy subspace. The $16$-component zero-energy modes at these critical momenta have a complicated expression and do not form a representation of $T_{\text{d}}$, thereby leading to non-uniform spinon condensation, as discussed above. Indeed, if condensation is restricted to one of the four inequivalent $\mathrm{L}$ momenta, we find that three of the four sublattices have the same spin amplitude, while the fourth sublattice has a different one.

V Experimental signatures

V.1 Critical behavior of the heat capacity

For each critical paraphase, the low-temperature heat capacity is expected to follow a power law whose exponent is determined by the low-energy spinon density of states in the critical theory. Indeed, depending on the dynamical exponent $z$ and the spinon condensation manifold (i.e., if spinons condense at points or along lines), this low-energy density of states follows different power laws $g(\epsilon)\sim\epsilon^{\alpha}$, where the possible values of $\alpha$ are listed in Table 6. The thermal energy due to spinon excitations at temperature $T$ is then given by

[TABLE]

and the heat capacity takes the form

[TABLE]

We remark that line condensation is not stable against generic perturbations, corresponding to further-neighbor terms in the mean-field ansatz. Consequently, at the lowest temperatures, we expect that the line-condensation paraphases are governed by the same exponents as their point-condensation counterparts. Nevertheless, if the NN mean-field ansatz is a good first approximation, there is an intermediate temperature range in which the approximate line condensation in such paraphases becomes manifest and therefore the line-condensation exponents in Table 6 are experimentally observable.

V.2 Critical spin structure factors

In this subsection, we present the most direct signatures of our critical points between magnetic orders and their parent spin liquids by computing both the static and the dynamic spin structure factors for the $15$ paraphases. While our calculation is based on mean-field theory, it still serves as a reference point for classifying the possible spinon spectra in pyrochlore magnets.

The static structure factor (SSF) is defined as the spatial Fourier transform of the equal-time spin-spin correlation function,

[TABLE]

where $\alpha=x,y,z$. We calculate this quantity using the critical mean-field ansätze in Sec. IV. Writing the $16\times 16$ matrix $V(\vec{k})$ in Eq. (38) as

[TABLE]

where the $8\times 8$ blocks generally satisfy

[TABLE]

due to charge-conjugation “symmetry”, the SSF becomes (see Appendix G for a detailed derivation)

[TABLE]

in terms of the auxiliary $8\times 8$ matrices

[TABLE]

where $I(\vec{q})$ is the $4\times 4$ diagonal form-factor matrix defined in Eq. (37). The resulting SSFs for representative points in each of the $15$ paraphases are plotted in Fig. 3 for chemical potentials $\mu=\mu_{c}+10^{-\delta}$, where $\delta=1,2,\dots,9$, and $\mu_{c}$ is the critical value given in Appendix E. When numerically computing the SSF, we ensure convergence by taking a momentum-space grid that does not contain any condensation momenta $\vec{k}_{c}$.

For chemical potentials well above the critical value $\mu_{c}$, the SSFs of the $15$ paraphases (not shown here) can be partitioned into two classes, depending on the sum of the $\mathbb{Z}_{2}$ parameters $n_{\overline{C}_{6}S}+n_{ST_{1}}+n_{\overline{C}_{6}}$ characterizing the parent spin liquid. Plotted along the high-symmetry path in the BZ, the SSFs of the $n_{\overline{C}_{6}S}+n_{ST_{1}}+n_{\overline{C}_{6}}=\text{odd}$ paraphases and those of the $n_{\overline{C}_{6}S}+n_{ST_{1}}+n_{\overline{C}_{6}}=\text{even}$ paraphases resemble each other after an appropriate reflection in energy. This relation between the two classes qualitatively survives as the chemical potential approaches its critical value (see Fig. 3). For example, depending on the the sum $n_{\overline{C}_{6}S}+n_{ST_{1}}+n_{\overline{C}_{6}}$ being even or odd, the SSF has either a valley or a peak at the $\Gamma$ point. The distinction between the two behaviors can be traced back to Eqs. (61) and (V.2). Since the SSF is the squared trace norm of the matrix $U^{\alpha}$, which in turn is the sum of two matrices $W^{\alpha}$, there is a cross term from the product of the two matrices $W^{\alpha}$, physically corresponding to the spinon pairing channel $\langle b^{\dagger}_{\vec{k}_{1},\mu\sigma_{1}}b^{\dagger}_{\vec{k}_{3},\nu\sigma_{3}}\rangle\langle b_{\vec{k}_{1}-\vec{q},\mu\sigma_{2}}b_{\vec{k}_{3}+\vec{q},\nu\sigma_{4}}\rangle$ (see Appendix G), and we numerically find this cross term to be negative for the $n_{\overline{C}_{6}S}+n_{ST_{1}}+n_{\overline{C}_{6}}=\text{even}$ paraphases and positive for the $n_{\overline{C}_{6}S}+n_{ST_{1}}+n_{\overline{C}_{6}}=\text{odd}$ paraphases. Nevertheless, a deeper understanding of this connection to $n_{\overline{C}_{6}S}+n_{ST_{1}}+n_{\overline{C}_{6}}$ remains to be found.

Also, there are general differences between the SSFs of the paraphases governed by $z=1$ and $z=2$ critical theories, respectively. For most of the $z=1$ paraphases, as the chemical potential approaches its critical value, the SSF becomes a non-differentiable function at certain momenta $\vec{q}$. This feature is clearly observable in Fig. 3 for the paraphases 0-(010)$\Lambda$, 0-(100)$\Lambda$, 0-(101)$\mathrm{W}$, 0-(111)$\mathrm{W}$, 0-(111)$\mathrm{X}$ at the $\Gamma$ point and for the paraphases 0-(101)$\mathrm{W}$, 0-(101)$\mathrm{X}$, 0-(111)$\mathrm{W}$, 0-(111)$\mathrm{X}$ at the $\mathrm{X}$ point. However, not all $z=1$ paraphases conform to this rule; the SSFs of the paraphases 0-(110)$\Gamma$ and 0-(110)$\Lambda$ do not reveal any singular behavior along the high-symmetry path in the BZ. Instead, they resemble the SSFs of $z=2$ paraphases, which are smooth across the entire BZ.

To understand these features, we consider the dynamic structure factor (DSF), which provides further information on the dynamics of spinons. This quantity is defined as the spatial and temporal Fourier transform of the spin-spin correlation function,

[TABLE]

and, using the mean-field ansätze in Sec. IV, it takes the general form (see Appendix G for a detailed derivation)

[TABLE]

The critical ($\mu=\mu_{c}$) DSFs and the corresponding spinon spectra are plotted in Fig. 4 for representative points in each of the $15$ paraphases.

Focusing on universal low-energy features, we first observe that each DSF has characteristic points or regions where it is gapless (i.e., finite at small $\omega$). Since the DSF describes spin dynamics, and each spin is decomposed into two spinons, the DSF is gapless at momenta $\vec{q}$ that are appropriate sums of spinon condensation momenta $\vec{k}_{c}$ such that $\vec{q}=\vec{k}_{c,1}+\vec{k}_{c,2}$. Consequently, we can establish a one-to-one correspondence between the potential spinon condensation momenta described in Table 4 and the gapless points or regions of the DSF plotted in Fig. 4; see Table 7 for this correspondence.

We also notice that the DSF has different low-energy behavior in the paraphases governed by $z=1$ and $z=2$ critical theories, respectively. For the $z=1$ paraphases 0-(010)$\Lambda$, 0-(100)$\Lambda$, 0-(101)$\mathrm{W}$, 0-(111)$\mathrm{W}$, and 0-(111)$\mathrm{X}$, the weight of the low-energy DSF is concentrated around zero energy, while for the $z=2$ paraphases 0-(001)$\Gamma$, 0-(100)$\Gamma$, 0-(101)$\Gamma$, and 0-(111)$\Gamma$, the low-energy DSF gradually vanishes as the energy is decreased to zero.

These low-energy features in the DSF can be understood from a scaling analysis of the critical field theories in Eqs. (47) and (51). The DSF is the expectation value of a four-point correlation function in the condensation fields; using Wick’s theorem, this expectation value can be written as the convolution of two Green’s functions,

[TABLE]

where $G_{z}(\omega,\vec{k})$ are labeled by the dynamical critical exponents $z$ of corresponding field theories:

[TABLE]

Inserting $G_{z}(\omega,\vec{k})$ into Eq. (65), and evaluating the integrals over $\omega^{\prime}$ and $\vec{k}$, we obtain the scaling behaviors

[TABLE]

where $f_{1}$ and $f_{2}$ are some general functions. At the momentum with gapless DSF, corresponding to $\vec{q}=0$, the DSF at $\omega\rightarrow 0$ thus diverges in the $z=1$ case and vanishes in the $z=2$ case. This result qualitatively explains the low-energy DSF features described above. Furthermore, it elucidates why the SSFs of the $z=1$ paraphases have singularities at specific momenta, which precisely coincide with the gapless momenta of the corresponding DSFs. We stress again that there are $z=1$ paraphases, for example, 0-(110)$\Gamma$, which have SSFs and DSFs following $z=2$ behavior. Such a discrepancy may occur if a coefficient in the critical theory accidentally vanishes for the NN mean-field ansatz.

Finally, we remark that the DSFs of several paraphases have high-energy points exhibiting large spectral weights at the $\Gamma$ point. In fact, whenever such points exist, there is a quasi-mirror-reflection symmetry (in terms of energy) between two spinon bands, such that the two band energies satisfy $\lambda_{\vec{k},1}+\lambda_{\vec{k},2}=E$ for all momenta $\vec{k}$. Due to this “symmetry”, these two bands can contribute strongly at the $\Gamma$ point close to energy $E$, resulting in an increased spectral weight as well as a Dirac-like texture. However, we emphasize that the high-energy part of the DSF depends on specific details and is not to be taken too seriously; only the low-energy part of the DSF captures the universal physics in the given paraphase.

V.3 General order parameters: hidden and intertwined orders

The naïve spin-condensation analysis of magnetic orders in Sec. IV.6 is far from complete as it presumes that any zero-momentum order transforms under a representation of $T_{\text{d}}$ and thus ignores the possibility of hidden orders transforming under inversion-odd representations of the full pyrochlore point group $O_{\text{h}}=T_{\text{d}}\times C_{\text{i}}$, where $C_{\text{i}}$ is the $\mathbb{Z}_{2}$ group consisting of inversion and identity. One simple example of such a hidden order is the alternating expansion and contraction of tetrahedra realized in the “breathing” pyrochlores Okamoto et al. (2013); Kimura et al. (2014); Rau et al. (2016b). In this subsection, we analyze zero-momentum orders more comprehensively by identifying all possible order parameters in terms of the condensing spinon fields and constructing the most general Ginzburg-Landau (GL) theory that is compatible with the point group $O_{\text{h}}$ of the pyrochlore lattice. Such an analysis has been previously done for several problems building on the PSG framework Bergman et al. (2006); Balents et al. (2005).

When the spinons condense at the critical point, certain bosonic modes at the condensation momenta $\vec{k}_{c}$ become macroscopically occupied, and the expectation values of their bosonic operators thus become classical condensation fields. For the $z=2$ critical points, the condensation fields $\phi_{n}$ are complex, while for the $z=1$ critical points, the condensation fields $\chi_{n}$ are real. Importantly, these fields themselves are not valid order parameters as they carry a $\mathbb{Z}_{2}$ gauge charge and transform projectively under the point group. Indeed, the projective transformation rules of $\phi_{n}$ and $\chi_{n}$ under the generators of the point group can be explicitly obtained from the corresponding transformation rules of the original bosonic operators $b_{\vec{k}_{c},\mu}$ (see also Appendix F):

[TABLE]

where $C_{3}(\mu)=0,2,3,1$ and $S(\mu)=3,1,2,0$ for the respective sublattices $\mu=0,1,2,3$. The simplest possible order parameters are then the bilinears of the condensation fields, corresponding to total momentum $\vec{K}=0$, which are gauge invariant and transform as linear, generically reducible, representations of the point group. For each paraphase, the irrep decomposition of this reducible representation is given in Table 8. We now discuss the physical implications of this decomposition.

The scalar irrep $A_{1\mathrm{g}}$ corresponds to a quadratic invariant, i.e., a “mass” term in the GL theory, which drives the phase transition between the spin-liquid phase and the magnetically ordered phase. For almost all paraphases, it appears only once in the reducible representation, which indicates that all components of the condensation occur together by symmetry. The bilinear term transforming under the scalar irrep is $\sum_{n}\chi_{n}^{2}$, where we decompose any complex fields into real fields as $\phi_{n}=\chi_{2n-1}+i\chi_{2n}$. The effective GL theory governing the phase transition is then

[TABLE]

When $A_{1\mathrm{g}}$ appears more than once in the reducible representation [for the paraphases [math]-$(101)\mathrm{W}$ and [math]-$(111)\mathrm{X}$], it signals an accidental degeneracy, which should be lifted when going beyond the NN level.

The remaining irreps, denoted by standard labels, correspond to various order parameters that describe distinct scenarios of symmetry breaking (see Table 9). Irreps with the subscript “$\mathrm{g}$” are even under inversion and correspond to the conventional spin orders discussed in Ref. Yan et al. (2017). The single-spin order parameters of such spin orders are straightforward to detect with neutron scattering. In contrast, irreps with the subscript “$\mathrm{u}$” are odd under inversion and correspond to more unconventional hidden orders. The order parameters of these inversion-breaking orders always contain multiple spin operators and are thus harder to detect Harter et al. (2017). However, in our case, they are also accompanied by a spontaneous breaking of inversion symmetry, which may be observed as a “breathing” distortion of the pyrochlore lattice.

Table 8 indicates that one paraphase can give rise to several distinct order parameters. In general, the presence or absence of a given order parameter is determined by the particular form of the GL theory governing the phase transition. However, for some paraphases, we can argue that several distinct orders are intertwined in the sense that they always accompany one another, regardless of the GL parameters. This highly nontrivial result emerges because the magnetically ordered phases are obtained by condensing fractionalized excitations (spinons) that transform projectively under symmetries.

To analyze the general intertwining between distinct orders for a given paraphase, we form an orthogonal basis for the (real) order parameters $\{\Psi_{R,1},\cdots,\Psi_{R,N_{R}}\}$ that transform under each distinct irrep $R$. Note that $N_{R}$ is the product of the irrep dimension and the multiplicity of the irrep in the reducible representation. Since each symmetry acts on the vector $(\Psi_{R,1},\cdots,\Psi_{R,N_{R}})$ by an orthogonal matrix, the quadratic term $W_{R}=\sum_{j=1}^{N_{R}}\Psi_{R,j}^{2}$ must be a scalar transforming under $A_{1\mathrm{g}}$. This scalar can be interpreted as the “weight” of the given irrep; since it is a function of the condensation fields $\chi_{n}$, it may vanish for some special configurations of these fields, indicating the absence of the corresponding order. In contrast, the total weight of all irreps,

[TABLE]

is nonzero for all field configurations, indicating that at least one order must always be present.

For each paraphase, however, the irreps $R$ may be partitioned into classes $C_{1}\oplus\cdots\oplus C_{q}$ (see Table 8) such that the total weight of each class $C=\{R_{1},\cdots,R_{N_{C}}\}$, containing some nontrivial subset of all irreps, is proportional to the total weight of all irreps,

[TABLE]

and is thus nonzero for all configurations of the condensation fields. Consequently, at least one order from each class $C$ must always be present, regardless of the GL parameters. In the most extreme scenario, when each irrep forms its own class, such that $W_{R}\propto W_{0}$ for all irreps $R$, the orders are maximally intertwined, i.e., all of them must appear together. For certain paraphases, one can argue for this scenario by counting all possible quadratic scalars that can be formed from the order parameters or, equivalently, all possible fourth-order scalars that can be formed from the condensation fields. There is always at least one such scalar, $(\sum_{n}\chi_{n}^{2})^{2}$; however, if there is only one such scalar, it is clear that the weight $W_{R}$ of each irrep $R$ must be proportional to this scalar, and all orders must therefore be simultaneously present.

While we do not analyze the general intertwining between distinct orders in all paraphases, we observe from the particular examples studied (see Table 8) that the presence of intertwined orders is a common feature of magnetically ordered phases obtained by spinon condensation on the pyrochlore lattice. In particular, for parent spin liquids with $n_{\overline{C}_{6}}=1$, where inversion symmetry acts projectively on the spinons, we generically anticipate the (already intertwined) spin orders to be also accompanied by inversion-breaking hidden orders.

VI Discussion

VI.1 Summary

In this paper, we gave a complete classification of spin-orbit-coupled $\mathbb{Z}_{2}$ spin liquids on the pyrochlore lattice by using the PSG method for Schwinger bosons. We studied the mean-field Hamiltonians of the six [math]-flux spin liquids at the NN level and examined the critical field theories that describe phase transitions to ordered phases via spinon condensation. We found two crucially different classes of critical field theories, characterized by dynamical exponents $z=1$ and $z=2$, respectively, which have distinct properties ranging from Hamiltonian diagonalizability to experimental observables. Moreover, we investigated the zero-momentum orders obtained from spinon condensation, both by a naïve spin-condensation analysis and by the representation theory of the full pyrochlore point group $O_{\text{h}}$. We found that seemingly unrelated orders are generically intertwined with each other and that conventional spin orders are often accompanied by more exotic inversion-breaking “hidden” orders. Finally, we calculated several physical observables for our critical theories, including the heat capacity, as well as the static and dynamic spin structure factors, which may be compared with experimental data.

VI.2 Possible implications

Many pyrochlore materials have been experimentally confirmed to possess one of the spin orders discussed in this paper. For example, Yb2Pt2O7 has ferromagnetic order Cai et al. (2016), while Nd2Zr2O7 possesses all-in-all-out order Lhotel et al. (2015). Since all of these spin orders can appear as a result of spinon condensation from one of our $\mathbb{Z}_{2}$ spin liquids, one can contemplate the possibility that some of these materials are proximate to such a spin liquid.

As a particular example, one may consider Er2Ti2O7, which is confirmed to have a $\Psi_{2}$ antiferromagnetic ground state. The $\Psi_{2}$ ground state is selected from the $\Gamma_{5}$ irrep, containing both $\Psi_{2}$ and $\Psi_{3}$ states, as a result of order-by-disorder mechanism, possibly aided by virtual crystal-field effects Henley (1989); Savary et al. (2012); Oitmaa et al. (2013); Javanparast et al. (2015b); Rau et al. (2016c). This $\Psi_{2}$ ground state is quite stable, which suggests that, if it is obtained from an instability of a spin liquid, such an instability should uniquely prefer $E_{\rm g}$ order. Consulting Table 8, we see that the paraphase $\Gamma$ of the PSG class 0-(110) has a single nontrivial irrep $E_{\text{g}}$, which is not intertwined with any other orders. Hence, if Er2Ti2O7 is proximate to a spin liquid, a natural candidate for its parent spin liquid is the one corresponding to the PSG class 0-(110).

One motivation of this paper was to understand the puzzling experiments on Yb2B2O7, where B = Ge, Ti, Sn. These three compounds have distinct ground states: the Ge compound is antiferromagnetic Hallas et al. (2016a), while the Ti Gaudet et al. (2016) and Sn Yaouanc et al. (2013) compounds are ferromagnetic, at least when any order can be clearly identified. The Ti compound is also sensitive to disorder. Despite the disparate ground states, inelastic neutron scattering gives very similar spectra for all three materials Hallas et al. (2016b), consisting of continuum weight over the entire Brillouin zone down to the lowest energies resolvable in the measurements. This observation suggests that the relevant excitations are characteristic of some common underlying structure, which is distinct from the usual spin waves tied to the individual ordered states. The approach in this paper gives one possible explanation: the excitations may be the spinons of a parent spin-liquid state.

To identify a potential parent spin liquid, we seek a PSG class from which both antiferromagnetic and ferromagnetic orders can be obtained through the same condensation paraphase. It is clear from Table 1 that such classes exist; the classes 0-(001), 0-(101) and 0-(111) all satisfy this criterion. Therefore, the proximity to a spin liquid corresponding to either of these classes can potentially explain the observed excitation spectra. Looking at the dynamic spin structure factors in Fig. 4, we indeed see that many of the critical structure factors in these classes [e.g., 0-(001)L, 0-(101)W, and 0-(111)X] have a large scattering continuum over the entire Brillouin zone down to a very small fraction of the spin-excitation bandwidth. It would be interesting to attempt a more quantitative comparison with the experimental data, which would require, at the very least, a careful consideration of effects beyond mean-field theory.

If the scattering continua in the Yb pyrochlores are reflections of a parent spin liquid, it also suggests that hidden order may be present in these materials Hallas et al. (2016b). Indeed, from the last column of Table 1, we see that the paraphases 0-(001), 0-(101) and 0-(111) all include hidden orders breaking inversion symmetry. Searching for such inversion-breaking orders may be an incisive test of the physical picture presented in this paper; if such an order is identified, a full characterization may be assisted by the associated order parameters in Table 9. We note that hidden order may also participate in the specific-heat anomalies of the Yb pyrochlores Hallas et al. (2016b).

VI.3 Future directions

The present paper explored the physics of proximity to a broad class of quantum spin liquids on the pyrochlore lattice. Nevertheless, several assumptions in the analysis could be modified or relaxed in future work. We focused on $\mathbb{Z}_{2}$ spin liquids and used the framework of bosonic spinons; it would be interesting to consider $U(1)$ spin liquids and explore fermionic spinons as well. The fermionic approach does not, however, provide a simple mean-field way to study magnetic instabilities, which is straightforward with bosonic spinons by condensing them.

In addition, the PSG results may be further exploited even within the framework of bosonic spinons. We concentrated on the [math]-flux NN mean-field Hamiltonians for simplicity, assuming that NNN terms do not qualitatively change our results. This assumption, however, is not necessarily true; in certain cases, a NNN term one-tenth as strong as a NN term can already change the condensation momenta. Moreover, the $\pi$-flux PSG classes may exhibit interesting physics of their own. These PSG classes have a fourfold enlarged unit cell due to nontrivial translational PSG along the $\hat{e}_{2}$ and $\hat{e}_{3}$ directions, which leads to a $64\times 64$ mean-field Hamiltonian in terms of the parameters in Table 2. Multi-spinon condensation may further enlarge the magnetic unit cell. In turn, this enlargement results in a complex spinon spectrum that probably requires a more computational approach.

We also presumed that the full symmetry group of the pyrochlore lattice is preserved at the level of the spin Hamiltonian. However, there is a large family of “breathing” pyrochlore materials Okamoto et al. (2013); Kimura et al. (2014); Rau et al. (2016b) that explicitly break inversion symmetry ($Fd\overline{3}m\rightarrow F\overline{4}3m$) by expansion and contraction of alternating tetrahedra. One material from this family, Ba3Yb2Zn5O11, was reported to remain disordered down to 0.38 K Kimura et al. (2014), and a gauge mean-field theory, distinct from the spinon approach in this paper, predicts that this material may experience a non-symmetry-breaking transition between a paramagnet and a quantum spin ice Savary et al. (2016). It would be interesting to see how this material (and the phase transition predicted for it) fits into a pyrochlore PSG classification.

Finally, the PSG method can be connected to the energetics of realistic spin Hamiltonians. Indeed, our mean-field spinon states can in principle be used as variational wave functions, as can their so far unexplored fermionic counterparts. Calculating variational energies for these wave functions would require a major effort in variational Monte Carlo in three dimensions; it is well beyond the present work but is quite worthwhile to explore.

Acknowledgements.

We acknowledge Yi-Zhuang You, Yuan-Ming Lu, and Bill Jacob for useful discussions. The work of G.B.H. was supported at ORNL by Laboratory Director’s Research and Development funds and at the KITP by the Gordon and Betty Moore Foundation’s EPiQS Initiative through Grant No. GBMF4304. C.L. and L.B. were supported by the DOE Office of Science’s Basic Energy Sciences program under Award No. DE-FG02-08ER46524.

Appendix A Point-group structure

The space group $Fd\overline{3}m$ belongs to the cubic crystal system with point group $O_{\text{h}}$. The point group $O_{\text{h}}$ has order 48 and is the symmetry group of a pyrochlore primitive cell – a pair of corner sharing tetrahedra. It has a direct product structure $O_{\text{h}}\cong S_{4}\times\mathbb{Z}_{2}$, which can be understood as following.

We label the seven vertices by $\mu^{\pm}$, where $\mu=0,1,2,3$ is the sublattice index and “$+$” (“$-$”) denotes the upper (lower) tetrahedon (where $0^{+}=0^{-}$ is the shared corner), then the symmetry operations in $O_{\text{h}}$ are permutations over two sets $\{0,1,2,3\}$ and $\{+,-\}$. The generators include a 3-fold rotation $C_{3}=(123)$, a screw operation (modding out translations) $S=(03)(+-)$ and an inversion $I=(+-)$, written in terms of the cycle notation for permutations. We also define the operation $\Sigma=S\circ I=(03)$ for future convenience. (We can also define $\overline{C}_{6}=C_{3}\circ I$ to reduce the number of generators, since equivalently $C_{3}=\overline{C}_{6}^{4}$ and $I=\overline{C}^{3}_{6}$.) The inversion $I$ is the generator of the $\mathbb{Z}_{2}$ group, therefore we can write $O_{\text{h}}\cong S_{4}\cup(I\circ S_{4})$, where $I\circ S_{4}$ is the coset of $S_{4}$ left-composed by $I$. The subgroup $S_{4}$ corresponds exactly to the tetrahedron group, $T_{\text{d}}$. The 24 elements of the group $S_{4}\simeq T_{\text{d}}$ are generated by $\Sigma$ and $C_{3}$ as (where it is understood $4\equiv 0$)

[TABLE]

Eqs. (72) will be useful in determining the mean field ansätze parameter constraints for the PSG classes.

Appendix B Solving PSG equations

The space group part of the PSG equations are

[TABLE]

The corresponding phase equations are

[TABLE]

where both Eq. (74a) and Eq. (74d) stand for three equations $i=1,2,3$, and Eq. (74e) stands for two equations $i=1,2$.

First we solve Eq. (74a). Due to gauge freedom of second type, we can use a gauge transformation to achieve $\phi_{T_{1}}(r_{1},r_{2},r_{3})_{\mu}=\phi_{T_{2}}(0,r_{2},r_{3})_{\mu}=\phi_{T_{3}}(0,0,r_{3})=0$. Then Eq. (74a) gives

[TABLE]

Using Eq. (75) to solve Eq. (74d)

[TABLE]

we get $n_{1}=n_{2}=n_{3}$, and

[TABLE]

Then using Eq. (75) to solve Eq. (74e) and (74f), we get

[TABLE]

Using Eqs. (75), (77) and (78) to solve Eq. (74g) and Eq. (74h) we get

[TABLE]

and

[TABLE]

Then, from Eq. (74b) we get

[TABLE]

and Eq. (74c) gives

[TABLE]

and

[TABLE]

Eqs. (80a) and (82) further imply that

[TABLE]

This completes solving the inter-unit cell part of the space group PSG equations. We can use some of the remaining gauge freedom to simplify results. In order to use the IGG freedom we notice that Eq. (73c), (73d) and (73e) have operators that appear odd number of times. According to our analysis in the main text, we can set $n_{S}=0$, and two out of the three parameters $n_{\overline{C}T_{i}}$ to be zero, which together with Eq. (84) means that $n_{\overline{C}_{6}T_{1}}=n_{\overline{C}_{6}T_{2}}=n_{\overline{C}_{6}T_{3}}=0$. The independent $\mathbb{Z}_{2}$ parameters at this point are

[TABLE]

Then we add time reversal operation. From

[TABLE]

where $\mathcal{O}\in\{T_{1},T_{2},T_{3},\overline{C}_{6},S\}$, we get

[TABLE]

where Eq. (87a) stands for three equations $i=1,2,3$. From Eq. (87a) we get

[TABLE]

From Eq. (87b) we get $n_{\mathcal{T}T_{1}}=n_{\mathcal{T}T_{2}}=n_{\mathcal{T}T_{3}}\equiv n_{\mathcal{T}T}$, and

[TABLE]

therefore $n_{\mathcal{T}\overline{C}_{6}}+n_{\mathcal{T}T}=n_{\overline{C}_{6}}$. Finally Eq. (87c) gives $n_{\mathcal{T}T}=0$, $n_{ST_{1}}=n_{ST_{2}}=n_{\mathcal{T}S}-n_{1}$, and

[TABLE]

Lastly the equation $\mathcal{T}^{2}=-1$ gives no constraint.

Now we have solved all the inter-unit cell part of the PSG equations. The intra-unit cell part gives

[TABLE]

Then we use the gauge freedom of second type. Note under gauge transformation

[TABLE]

we have $\phi_{\mathcal{O}}(\vec{r}_{\mu})\rightarrow\phi_{\mathcal{O}}(\vec{r}_{\mu})+\phi(\vec{r}_{\mu})-\phi[\mathcal{O}^{-1}(\vec{r})_{\mu}]$, where $\phi$ is an artitrary $U(1)$ phase, we have $\phi_{\overline{C}_{6}}(\vec{0})_{0}\rightarrow\phi_{\overline{C}_{6}}(\vec{0})_{0}$, $\phi_{\overline{C}_{6}}(\vec{0})_{i}\rightarrow\phi_{\overline{C}_{6}}(\vec{0})_{i}+\phi_{i}-\phi_{i+1}$, $\phi_{S}(\vec{0})_{0}\rightarrow\phi_{S}(\vec{0})_{0}+\phi_{0}-\phi_{3}$, $\phi_{S}(\vec{0})_{1,2}\rightarrow\phi_{S}(\vec{0})_{1,2}$, $\phi_{S}(\vec{0})_{3}\rightarrow\phi_{S}(\vec{0})_{3}+\phi_{3}-\phi_{0}$, and $\phi_{\mathcal{T}}(\vec{0})_{\mu}\rightarrow\phi_{\mathcal{T}}(\vec{0})_{\mu}+2\phi_{\mu}$. Then, we can choose the value of $\phi_{\mu}$ to fix

[TABLE]

and

[TABLE]

where $p_{\mu},m_{0}$ and $m_{1,2}$ are all $\mathbb{Z}_{2}$ parameters.

Note we still have a discrete gauge freedom: we can choose a particular sublattice $\nu$ and define gauge transformation

[TABLE]

then Eq. (93) is preserved but the relative phase of $\phi_{\overline{C}_{6}}$ can be changed. By choosing $\nu=1,2,3$ we can use gauge (95) to fix $p_{1}=p_{2}=p_{3}\equiv p$. Furthermore, we can use the global $\mathbb{Z}_{2}$ freedom for $\phi_{\overline{C}_{6}}(\vec{r}_{\mu})$ and $\phi_{S}(\vec{r}_{\mu})$ to fix $p_{0}=0$ and $m_{1}=0$. Then, let $\nu=0,3$, we can use gauge (95) to fix $m_{0}=0$. By checking Eqs. (91), we have $n_{S\overline{C}_{6}}=n_{\overline{C}_{6}}+n_{1}+n_{ST_{1}}$, $p=n_{1}+n_{ST_{1}}$ and $m_{2}=n_{\overline{C}_{6}S}$. The final solution is presented in Eq. (22).

Appendix C Basis for irreps of $T_{d}$

This appendix gives the representation analysis result for the spins $\mathbf{S}$ on a single tetrahedron, which can be equally applied to pyrochlore lattices with a $\Gamma$ point order. The twelve-component spin $\mathbf{S}$ form a 12-dimensional representation of the tetrahedon group $T_{\text{d}}$. The group $T_{\text{d}}$ has irreducible representation (irrep)$A_{1},A_{2},E,T_{1},T_{2}$. $\mathbf{S}$ can be decomposed into irreps $A_{2},E,T_{1,A},T_{1,B}$ and $T_{2}$. The corresponding basis and orders are listed in Table. 10. This is simply a reproduction of TABLE III in Ref. Yan et al. (2017).

The basis are

[TABLE]

Appendix D Derivation of the mean-field Hamiltonians

In this section we present the solution of mean field constraints from PSG classes, up to NNN level.

For the on-site bond $\vec{0}_{0}\rightarrow\vec{0}_{0}$, the 12 group elements that map the bond back are

[TABLE]

which give constraints

[TABLE]

The only allowed onsite hopping term is $\alpha$ which is simply the chemical potential $\mu$. The allowed pairing terms are $\beta^{\prime}=\gamma^{\prime}=\delta^{\prime}\equiv\nu$ when and only when $(n_{1}+n_{ST_{1}},n_{\overline{C}_{6}})=(1,0)$. Note the singlet pairing term is not allowed, $\alpha^{\prime}=0$, due to hermiticity.

Consider the case for the NN bond $\vec{0}_{0}\rightarrow\vec{0}_{1}$. The four group elements that map the bond back are [see Appendix A for details]

[TABLE]

using Eqs.(72), (24), (26) and (30), we get constraints

[TABLE]

we get

• •

$n_{\overline{C}_{6}S}=0$: $b=0$, $c=-d$, $a$, $c$ free independent; $a^{\prime}=c^{\prime}=d^{\prime}=0$, we have

• –

$n_{1}+n_{ST_{1}}=n_{\overline{C}_{6}}$: no NN pairing term allowed;

• –

$n_{1}+n_{ST_{1}}=n_{\overline{C}_{6}}+1$: $b^{\prime}$ free.

• •

$n_{\overline{C}_{6}S}=1$: $a=c=d=0$, $b$ free; $b^{\prime}=0$, and

• –

$n_{1}+n_{ST_{1}}=n_{\overline{C}_{6}}$: $a^{\prime}=0$, $c^{\prime}=d^{\prime}$, $c^{\prime}$ free;

• –

$n_{1}+n_{ST_{1}}=n_{\overline{C}_{6}}+1$: $c^{\prime}=-d^{\prime}$, $a^{\prime}$, $c^{\prime}$ free independent.

For the NNN bond $\vec{0}_{1}\rightarrow\vec{0}_{2}-\hat{e}_{2}$, it can be checked that only the identity $(1)$ and the element $(12)(+-)=S\circ C_{3}\circ S\circ C^{-1}_{3}\circ S\circ C_{3}$ map the bond back, which gives

[TABLE]

Therefore for hopping

• •

$n_{1}=0$: $B=C$, $A,B,D$ free independent;

• •

$n_{1}=1$: $A=D=0$, $B=-C$, $B$ free,

and for pairing

• •

$n_{ST_{1}}=0$: $A^{\prime}=D^{\prime}=0$, $B^{\prime}=-C^{\prime}$, $B^{\prime}$ free;

• •

$n_{ST_{1}}=1$: $B^{\prime}=C^{\prime}$, $A^{\prime},B^{\prime},D^{\prime}$ free independent.

These results are listed in Table 2.

Appendix E Critical chemical potential $\mu$

The critical chemical potential $\mu_{c}$ for the 15 paraphases is listed in Table 11.

Appendix F Condensation results

The three vectors $\mathbf{S}^{r}$, $\mathbf{S}^{c}$, and $\mathbf{S}^{s}$ for the paraphase 0-$(100)\Gamma$, mentioned in Eq. (53), are

[TABLE]

The three vectors $\mathbf{S}^{r}$, $\mathbf{S}^{c}$, and $\mathbf{S}^{s}$ for the paraphase 0-$(001)\Gamma$, mentioned in Sec. IV.6.3, are

[TABLE]

In writing down the Ginzburg-Landau theory for the paraphase 0-(010)$\Gamma$, the transformation rules of $\phi_{1},\phi_{2},\overline{\phi}_{1},\overline{\phi}_{2}$ are

[TABLE]

[TABLE]

where the definition of $\delta,\Delta$ and $\zeta$ has been given in Section IV.6. We see that in this case the fields transform to their complex conjugates under $O_{\text{h}}$.

The only quartic term invariant under $O_{\text{h}}$ is

[TABLE]

In writing down the Ginzburg-Landau theory for the paraphase 0-(100)$\Gamma$, the transformation rules of $\phi_{1,2}$ under $\overline{C}_{6}$ and $S$ are recorded by the following matrices:

[TABLE]

There are six quartic terms that are invariant under $O_{\text{h}}$. Three of them can be written as $\Phi_{i}^{2}$, where

[TABLE]

In writing down the Ginzburg-Landau theory for the paraphase 0-(101)$\Gamma$, the transformation rules of $\phi_{1,2}$ under $\overline{C}_{6}$ and $S$ are recorded by the following matrices:

[TABLE]

the extra factor of $i$ for $\overline{C}_{6}$ is due to $n_{\overline{C}_{6}}=1$.

In writing down the Ginzburg-Landau theory for the paraphase 0-(110)$\Gamma$, the transformation rules of $\chi_{1,2}$ under $\overline{C}_{6}$ and $S$ are recorded by the following matrices:

[TABLE]

the only quadratic and quartic order parameter are the trivial one: $(\chi_{1}^{2}+\chi_{2}^{2})^{i}$, $i=1,2$. At sextic order, there are two terms allowed:

[TABLE]

In writing down the Ginzburg-Landau theory for the paraphase 0-(111)$\Gamma$, the transformation rules of $\phi_{1,2}$ under $\overline{C}_{6}$ and $S$ are recorded by the following matrices:

[TABLE]

Appendix G Structure factor

We present the calculation of the static and dynamic structure factors. Importantly, we must distinguish the phase factors in the above two expressions: the $e^{i\vec{q}\cdot(\vec{r}_{\mu}-\vec{r}^{\prime}_{\nu})}$ in structure factor $\mathcal{S}^{\alpha}$ is global, namely $\vec{r}_{\mu}=\vec{r}+\hat{\varepsilon}_{\mu}$, which keeps track the relative displacement between sublattices, and this comes from the definition of structor factor. The phase factor in Fourier transformation, on the other hand, must agree with the convention we choose in Fourier transforming the Hamiltonian into $\vec{k}$ space: remember that Block Hamiltonian sets the displacement between sublattices to zero, therefore we must also set the displacement between sublattices to zero in the Fourier transforms, i.e. $\vec{r}_{\mu}=\vec{r}$.

[TABLE]

From Eq. (59) we write

[TABLE]

one can show using Wick’s theorem that

[TABLE]

where the first, second, and third terms come from the channels $\left\langle b^{\dagger}_{\vec{k}_{1},\mu\sigma_{1}}b^{\dagger}_{\vec{k}_{3},\nu\sigma_{3}}\right\rangle\left\langle b_{\vec{k}_{1}-\vec{q},\mu\sigma_{2}}b_{\vec{k}_{3}+\vec{q},\nu\sigma_{4}}\right\rangle$, $\left\langle b^{\dagger}_{\vec{k}_{1},\mu\sigma_{1}}b_{\vec{k}_{3}+\vec{q},\nu\sigma_{4}}\right\rangle\left\langle b^{\dagger}_{\vec{k}_{3},\nu\sigma_{3}}b_{\vec{k}_{1}-\vec{q},\mu\sigma_{2}}\right\rangle$, and $\left\langle b^{\dagger}_{\vec{k}_{1},\mu\sigma_{1}}b_{\vec{k}_{1}-\vec{q},\mu\sigma_{2}}\right\rangle\left\langle b^{\dagger}_{\vec{k}_{3},\nu\sigma_{3}}b_{\vec{k}_{3}+\vec{q},\nu\sigma_{4}}\right\rangle$, respectively. In the expression for structure factors the third term becomes $\left\langle\hat{S}^{\alpha}_{\vec{r}_{\mu}}\right\rangle$ and vanishes due to time reversal symmetry. Therefore

[TABLE]

where we have used Eq. (60) and defined $I(\vec{q})$ as in Eq. (37). We used the notation “$\colon$” to denote that the corresponding rows (columns) are retained in the matrix: for example, the notaton $\left(V^{\dagger}_{12}(\vec{k}_{1})\right)_{\rho_{1}\tau_{1},\colon}$ denotes the $\rho_{1}\tau_{1}$-row of the matrix $V^{\dagger}_{12}(\vec{k}_{1})$ (where all columns in this row are retained), and the notation $\left(V_{11}(\vec{k}_{1}-\vec{q})\right)_{\colon,\rho_{2}\tau_{2}}$ denotes the $\rho_{2}\tau_{2}$-column of the matrix $V_{11}(\vec{k}_{1}-\vec{q})$ (where all rows in this column are retained).

Furthermore, Define $W^{\alpha}(\vec{k},\vec{q})=V^{\dagger}_{12}(\vec{k}_{1})\left(I(\vec{q})\otimes\sigma^{\alpha}\right)V_{11}(\vec{k}_{1}-\vec{q})$, we have

[TABLE]

Notice that if we pick the four terms relevant for a given $\vec{k}_{1}$ (we omit the $\alpha$ index for simplicity):

[TABLE]

therefore Eq. (115) can be written as Eqs. (61)–(V.2) in the main text.

The dynamic structure factor is defined in Eq. (63). We have

[TABLE]

in Heisenberg representation, $\widetilde{b}(t)=e^{iHt}\widetilde{b}e^{-iHt}=e^{i\lambda\widetilde{b}^{\dagger}\widetilde{b}t}\widetilde{b}e^{-i\lambda\widetilde{b}^{\dagger}\widetilde{b}t}=\widetilde{b}e^{-i\lambda t}$, we have

[TABLE]

again we neglect the first term which vanishes in dynamical structure factor due to time reversal symmetry. So that

[TABLE]

plug this equation into Eq. (117) we obatin Eq. (64).

We see that the dynamic structure factor simply “disperses” the static structure factor according to the energy level of each matrix element of $V^{\dagger}_{12}(\vec{k}_{1})\left(I(\vec{q})\otimes\sigma^{\alpha}\right)V_{11}(\vec{k}_{1}-\vec{q})$.