Rank-2 $U(1)$ spin liquid on the breathing pyrochlore lattice

TL;DR

This paper proposes a microscopic model for a rank-2 $U(1)$ spin liquid on the breathing pyrochlore lattice, linking higher-rank gauge theories to real materials and providing experimental predictions.

Contribution

It introduces a spin liquid model based on rank-2 $U(1)$ gauge theory in a realistic magnetic system, specifically Yb-based breathing pyrochlores.

Findings

Identification of Yb-based breathing pyrochlores as candidate materials.

Explicit experimental signatures predicted for the rank-2 $U(1)$ spin liquid.

Connection established between higher-rank gauge theories and material realizations.

Abstract

Higher-rank generalisations of electrodynamics have recently attracted considerable attention because of their ability to host "fracton" excitations, with connections to both quantum stabilizer codes and holography. However, the search for higher-rank gauge theories in experiment has been greatly hindered by the lack of materially-relevant microscopic models. Here we show how a spin liquid described by rank-2 $U(1)$ gauge theory can arise in a magnet on the breathing pyrochlore lattice. We identify Yb-based breathing pyrochlores as candidate systems, and make explicit predictions for how the rank-2 $U(1)$ spin liquid would manifest itself in experiment.

Tables

Table 1. Table 1: Order parameters 𝐦 𝖷 subscript 𝐦 𝖷 {\bf m}_{\mathsf{X}} , describing how the point-group symmetry of a single tetrahedron within the pyrochlore lattice is broken by magnetic order. Order parameters transform according to irreducible representations of the point-group 𝖳 d subscript 𝖳 𝑑 {\sf T}_{d} , and are expressed in terms of linear combinations of spin-components 𝐒 i = ( S i x , S i y , S i z ) subscript 𝐒 𝑖 subscript superscript 𝑆 𝑥 𝑖 subscript superscript 𝑆 𝑦 𝑖 subscript superscript 𝑆 𝑧 𝑖 {\bf S}_{i}=(S^{x}_{i},S^{y}_{i},S^{z}_{i}) , in the global frame of the crystal axes — cf. ℋ ℋ \mathcal{H}  [Eq. ( S15 ))]. Labelling of spins within the tetrahedron follows the convention of Ross  et al.   Ross  et al. ( 2011 ) . The notation Ψ i subscript Ψ 𝑖 \Psi_{i} for ordered phases is taken from Ref.  Poole  et al. ( 2007 ) .

order	definition in terms	associated
parameter	of spin components	ordered phases
$m_{{𝖠}_{\mathrm{𝟤}}}$	$\frac{1}{2\sqrt{3}}\left(S_0^x+S_0^y+S_0^z+S_1^x-S_1^y-S_1^z-S_2^x+S_2^y-S_2^z-S_3^x-S_3^y+S_3^z\right)$	“all in-all out”
${𝐦}_{𝖤}$	$\left(\begin{array}{c}\frac{1}{2\sqrt{6}}\left(-2S_0^x+S_0^y+S_0^z-2S_1^x-S_1^y-S_1^z+2S_2^x+S_2^y-S_2^z+2S_3^x-S_3^y+S_3^z\right)\\ \frac{1}{2\sqrt{2}}\left(-S_0^y+S_0^z+S_1^y-S_1^z-S_2^y-S_2^z+S_3^y+S_3^z\right)\end{array}\right)$	$\begin{array}{c}{\Gamma }_5,\text{including}\\ {\Psi }_2\text{ and }{\Psi }_3\end{array}$
${𝐦}_{{𝖳}_{\mathrm{𝟣}+}}$	$\left(\begin{array}{c}\frac{1}{2}\left(S_0^x+S_1^x+S_2^x+S_3^x\right)\\ \frac{1}{2}\left(S_0^y+S_1^y+S_2^y+S_3^y\right)\\ \frac{1}{2}\left(S_0^z+S_1^z+S_2^z+S_3^z\right)\end{array}\right)$	collinear FM
${𝐦}_{{𝖳}_{\mathrm{𝟣},-}}$	$\left(\begin{array}{c}\frac{-1}{2\sqrt{2}}\left(S_0^y+S_0^z-S_1^y-S_1^z-S_2^y+S_2^z+S_3^y-S_3^z\right)\\ \frac{-1}{2\sqrt{2}}\left(S_0^x+S_0^z-S_1^x+S_1^z-S_2^x-S_2^z+S_3^x-S_3^z\right)\\ \frac{-1}{2\sqrt{2}}\left(S_0^x+S_0^y-S_1^x+S_1^y+S_2^x-S_2^y-S_3^x-S_3^y\right)\end{array}\right)$	non-collinear FM
${𝐦}_{{𝖳}_{\mathrm{𝟤}}}$	$\left(\begin{array}{c}\frac{1}{2\sqrt{2}}\left(-S_0^y+S_0^z+S_1^y-S_1^z+S_2^y+S_2^z-S_3^y-S_3^z\right)\\ \frac{1}{2\sqrt{2}}\left(S_0^x-S_0^z-S_1^x-S_1^z-S_2^x+S_2^z+S_3^x+S_3^z\right)\\ \frac{1}{2\sqrt{2}}\left(-S_0^x+S_0^y+S_1^x+S_1^y-S_2^x-S_2^y+S_3^x-S_3^y\right)\end{array}\right)$	Palmer–Chalker (${\Psi }_4$)

Table 2. Table 2: Coefficients a 𝖷 subscript 𝑎 𝖷 a_{\mathsf{X}} of the irrep invariants | 𝐦 𝖷 | 2 superscript subscript 𝐦 𝖷 2 |{\bf m}_{\mathsf{X}}|^{2} appearing in ℋ ℋ \mathcal{H}  [Eq.( S22 )]. Coefficients are expressed as a function of J 𝐽 J and D 𝐷 D . Here the subscripts for the A- and B-tetrahedra are suppressed.

coefficient of ${\left|{𝐦}_{𝖷}\right|}^2$ in Eq.(S22)	definition in terms of $J$ and $D$
$a_{{𝖠}_{\mathrm{𝟤}}}$	$-J-4D/\sqrt{2}$
$a_{𝖤}$	$-J+2D/\sqrt{2}$
$a_{{𝖳}_{\mathrm{𝟤}}}$	$-J-2D/\sqrt{2}$
$a_{{𝖳}_{+}}$	$3J$
$a_{{𝖳}_{-}}$	$-J+2D/\sqrt{2}$

Equations

$${\partial }_i{E}_i=0.$$

$$E_{ji}=E_{ij}\mathrm{Tr}\mathbf{E}=0,$$

$${\partial }_i{E}_{ij}={\rho }_j.$$

$${\partial }_i{E}_{ij}=0.$$

$$A_{ij}\to A_{ij}+{\partial }_i{\lambda }_j+{\partial }_j{\lambda }_i+\gamma {\delta }_{ij},$$

$$\begin{array}{rl}{\mathcal{H}}_{\mathsf{BP}}=& \sum _{⟨ij⟩\in \text{A}}\left[J_A{\mathbf{S}}_i\cdot {\mathbf{S}}_j+D_A{\hat{\mathbf{d}}}_{ij}\cdot ({\mathbf{S}}_i\times {\mathbf{S}}_j)\right]\\ & +\sum _{⟨ij⟩\in \text{B}}\left[J_B{\mathbf{S}}_i\cdot {\mathbf{S}}_j+D_B{\hat{\mathbf{d}}}_{ij}\cdot ({\mathbf{S}}_i\times {\mathbf{S}}_j)\right].\end{array}$$

$$\mathcal{H}=\frac{1}{2}\sum _{\mathsf{tet}\in \mathsf{A},\mathsf{X}}{a}_{\text{A},\mathsf{X}}{m}_{\mathsf{X}}^2+\frac{1}{2}\sum _{\mathsf{tet}\in \mathsf{B},\mathsf{X}}{a}_{\text{B},\mathsf{X}}{m}_{\mathsf{X}}^2,$$

$$J_A=J_B,D_A=D_B=0,$$

$$0<a_{{\mathsf{A}}_2}=a_{\mathsf{E}}=a_{{\mathsf{T}}_2}=a_{{\mathsf{T}}_{1-}}<a_{{\mathsf{T}}_{1+}}.$$

$${\mathbf{E}}^{\mathsf{HAF}}={\mathbf{E}}_{sym.}^{\mathsf{HAF}}+{\mathbf{E}}_{antisym.}^{\mathsf{HAF}}+{\mathbf{E}}_{\mathsf{trace}}^{\mathsf{HAF}}$$

$${\mathbf{E}}_{sym.}^{\mathsf{HAF}}=\left[\begin{array}{ccc}\frac{2}{\sqrt{3}}{m}_{\mathsf{E}}^1& m_{{\mathsf{T}}_{1-}}^z& m_{{\mathsf{T}}_{1-}}^y\\ m_{{\mathsf{T}}_{1-}}^z& -\frac{1}{\sqrt{3}}{m}_{\mathsf{E}}^1-m_{\mathsf{E}}^2& m_{{\mathsf{T}}_{1-}}^x\\ m_{{\mathsf{T}}_{1-}}^y& m_{{\mathsf{T}}_{1-}}^x& -\frac{1}{\sqrt{3}}{m}_{\mathsf{E}}^1+m_{\mathsf{E}}^2\end{array}\right],$$

$$(E_{antisym.}^{\mathsf{HAF}}{)}_{ij}=-{ϵ}_{ijk}{m}_{{\mathsf{T}}_2}^k,(E_{\mathsf{trace}}^{\mathsf{HAF}}{)}_{ij}=-{\delta }_{ij}\sqrt{\frac{2}{3}}{m}_{{\mathsf{A}}_2}.$$

$$\frac{2}{\sqrt{3}}\left[\begin{array}{c}{\partial }_x{m}_{\mathsf{E}}^1\\ -\frac{1}{2}{\partial }_y{m}_{\mathsf{E}}^1-\frac{\sqrt{3}}{2}{\partial }_y{m}_{\mathsf{E}}^2\\ -\frac{1}{2}{\partial }_y{m}_{\mathsf{E}}^1+\frac{\sqrt{3}}{2}{\partial }_y{m}_{\mathsf{E}}^2\end{array}\right]-\left[\begin{array}{c}{\partial }_y{m}_{{\mathsf{T}}_{1-}}^z+{\partial }_z{m}_{{\mathsf{T}}_{1-}}^y\\ {\partial }_z{m}_{{\mathsf{T}}_{1-}}^x+{\partial }_x{m}_{{\mathsf{T}}_{1-}}^z\\ {\partial }_x{m}_{{\mathsf{T}}_{1-}}^y+{\partial }_y{m}_{{\mathsf{T}}_{1-}}^x\end{array}\right]$$

$$-\sqrt{\frac{2}{3}}\nabla m_{{\mathsf{A}}_2}+\nabla \times {\mathbf{m}}_{{\mathsf{T}}_2}=0.$$

$$J_A,J_B>0,D_A<0,D_B=0.$$

$$a_{\mathsf{E},\text{A}}=a_{{\mathsf{T}}_{1-},\text{A}}<a_{{\mathsf{A}}_2,\text{A}},a_{{\mathsf{T}}_2,\text{A}},a_{{\mathsf{T}}_{1+},\text{A}},$$

$$⟨E_i(\mathbf{q})E_j(-\mathbf{q})⟩\propto {\delta }_{ij}-\frac{q_i{q}_j}{q^2},$$

$$\begin{split}&\langle E_{ij}({\bf q})E_{kl}(-{\bf q})\rangle\propto\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})+\frac{q_{i}q_{j}q_{k}q_{l}}{q^{4}}\\ &\quad-\frac{1}{2}\bigg{(}\delta_{ik}\frac{q_{i}q_{l}}{q^{2}}+\delta_{jk}\frac{q_{i}q_{l}}{q^{2}}+\delta_{il}\frac{q_{j}q_{k}}{q^{2}}+\delta_{jl}\frac{q_{i}q_{k}}{q^{2}}\bigg{)}\\ &\quad-\frac{1}{2}\left(\delta_{ij}-\frac{q_{i}q_{j}}{q^{2}}\right)\left(\delta_{kl}-\frac{q_{k}q_{l}}{q^{2}}\right)\;.\end{split}$$

$$J_A=J_B=1,D_A=-0.01,D_B=0.$$

$$J_A=0.57\text{meV},J_B=0.028\text{meV},$$

$$D_A=-0.16\text{meV},D_B=-0.007\text{meV},$$

$$E_{ij}=E_{ji},E_{ii}=0.$$

$${\partial }_i{E}_{ij}=0$$

$$A_{ij}=A_{ji}.$$

$$|\Psi (\mathbf{A})⟩.$$

$$-i({\lambda }_j{\partial }_i{E}_{ij}+{\lambda }_j{\partial }_i{E}_{ij})|\Psi (\mathbf{A})⟩=0.$$

$$i({\partial }_i{\lambda }_j+{\partial }_i{\lambda }_j)E_{ij}|\Psi (\mathbf{A})⟩=0.$$

$$i({\partial }_i{\lambda }_j+{\partial }_i{\lambda }_j)E_{ij}|\Psi (\mathbf{A})⟩=|\Psi (\mathbf{A}+\nabla \otimes \mathbf{\lambda }+(\nabla \otimes \mathbf{\lambda }{)}^T)⟩-|\Psi (\mathbf{A})⟩=0.$$

$$\begin{array}{r}\mathbf{A}+\nabla \otimes \mathbf{\lambda }+(\nabla \otimes \mathbf{\lambda }{)}^T,\text{ i.e., }{A}_{ij}\to A_{ij}+{\partial }_i{\lambda }_j+{\partial }_i{\lambda }_j.\end{array}$$

$$-i\gamma {\delta }_{ij}{E}_{ij}|\Psi (\mathbf{A})⟩=0.$$

$$A_{ij}\to A_{ij}+\gamma {\delta }_{ij}.$$

Full text

*Introduction. * It is of great intellectual interest and practical utility to discover novel effective laws of nature emerging from many–body systems. Traditionally, this enterprise has been entwined with the concept of broken symmetry Anderson (1972). However, a powerful alternative has proved to be the local constraints which arise from competing or “frustrated”, interactions. In the context of frustrated magnets, these can lead to the emergence of a local gauge symmetry, and thereby to quantum spin liquids, which defy all usual concepts of magnetic order, and instead exhibit fractionalised excitations and long–range entanglement Anderson (1973); Balents (2010); Savary and Balents (2017); Zhou et al. (2017). A well–studied example is quantum spin ice, a realisation of a $U(1)$ gauge theory on the pyrochlore lattice, whose emergent excitations exactly mimic conventional electrodynamics: photons, electric charges and magnetic monopoles. As such, it has attracted intense theoretical Hermele et al. (2004); Banerjee et al. (2008); Benton et al. (2012); Savary and Balents (2012); Shannon et al. (2012); Hao et al. (2014); Gingras and McClarty (2014); Kato and Onoda (2015); Chen (2017); Huang et al. (2018) and experimental Zhou et al. (2008); Ross et al. (2011); Fennell et al. (2012); Kimura et al. (2013); Sibille et al. (2015); Wen et al. (2017); Thompson et al. (2017); Sibille et al. (2018); Gao et al. (2019) investigation.

Recent work has highlighted the possibility of more exotic forms of emergent electrodynamics Xu (2006); Pretko (2017a); Rasmussen et al. (2016); Pretko (2017b), where electric and magnetic fields have the form of rank–2 (or higher–rank) tensors. These theories have modified conservation laws and gauge symmetries, resulting in some remarkable properties. Some are argued to mimic gravity Xu (2006); Benton et al. (2016); Pretko (2017c), while others are dual to elasticity theory Pretko and Radzihovsky (2018); Gromov (2019). In both cases, the charged excitations, dubbed “fractons”, have constrained mobility, and characterize a new class of topological order Chamon (2005); Shannon et al. (2004); Haah (2011); Vijay et al. (2015, 2016); Bulmash and Barkeshli (2018); Ma et al. (2018); Nandkishore and Hermele (2019); Slagle and Kim (2017); Halász et al. (2017a). Fracton models are also linked to quantum stabilizer codes Schmitz et al. (2018); Kubica and Yoshida (1836) and holography Yan (2019). None the less, these desirable properties come at a price: the local constraint required has a tensor character. As a consequence, prototypical models of fractons require rather complicated interactions Chamon (2005); Xu (2006); Haah (2011); Vijay et al. (2015, 2016), with just a handful of proposals motivated by experiment Slagle and Kim (2017); Halász et al. (2017b); You and von Oppen (2019). In the case of gapless higher–rank gauge theories, only a few concrete models exist Xu and Hořava (2010); Xu and Fisher (2007); Rasmussen et al. (2016), and even less is known about how to achieve such a phase in a real material. For this reason, realizing an emergent higher–rank electrodynamics in experiment presents a significant challenge.

In this Letter we show how a canonical rank–2 $U(1)$ [R2–U1] spin liquid can arise in a realistic model of a frustrated magnet. The model we consider is the Heisenberg antiferromagnet (HAF) on a breathing–pyrochlore (BP) lattice, perturbed by weak Dzyaloshinskii–Moriya (DM) interactions [Fig. 1a]. Working in the classical limit, relevant to a spin liquid at finite temperatures, we establish that fluctuations can be described using a tensor field satisfying the constraints required for a R2–U1 gauge theory. We use classical Monte Carlo (MC) simulation to confirm this scenario, and to explore how a R2–U1 spin liquid could be identified in experiment. We find that 4–fold pinch points (4FPP), characteristic of the R2–U1 state Prem et al. (2018), become visible in polarised neutron scattering. We discuss the application of these ideas to real materials, identifying Yb–based breathing pyrochlores as potential candidates for an R2–U1 spin liquid state. These results complement earlier work exploring gapped, fracton topological order, in models with bilinear interactions Slagle and Kim (2017); Halász et al. (2017b); You and von Oppen (2019), providing an example of an R2–U1 state, in an experimentally–motivated context.

*Review of R2–U1 theory. * Conventional, $U(1)$, electrodynamics is built around a vector field $\bf E$, subject to a Gauss law $\partial_{i}E_{i}=\rho\;$ so that, in the absence of charges,

[TABLE]

The key which unlocked the effective electrodynamics of spin ice was the realisation that, at low temperatures, in a classical limit, spins satisfied a local constraint of precisely the form of Eq. (1) Harris et al. (1997); Bramwell and Gingras (2001).

Here we consider instead an R2–U1 electrodynamics, in its self–dual, vector–charged, traceless form Pretko (2017a); Rasmussen et al. (2016); Pretko (2017b), and seek to show that, in an equivalent classical limit, spins satisfy the appropriate generalisation of Eq. (1). The R2–U1 theory is built around a rank–2 tensor electric field $\bf E$ that is symmetric and traceless,

[TABLE]

subject to a generalised Gauss’ law for a vector charge

[TABLE]

In the low–energy sector, the theory is charge free, i.e.

[TABLE]

These constraints determine the symmetry of the R2–U1 gauge field

[TABLE]

which in turn implies the form of the associated magnetic field, $B_{ij}$ Pretko (2017a, b). However the key observable properties of an R2–U1 spin liquid follow from the correlations of its electric field $E_{ij}$ Prem et al. (2018), and our goal will therefore be to show how the spins in a frustrated magnet can be described by a tensor field $E_{ij}$, satisfying the constraints Eqs. (2, 4).

*The model. * To this end, we consider a HAF, perturbed by weak DM interactions, on a “breathing” pyrochlore (BP) lattice, for which A– and B–sublattice tetrahedra have a different size

[TABLE]

Definitions of the bond–dependent vectors $\hat{\bf d}_{ij}$ Kotov et al. (2005); Poole et al. (2007); Elhajal et al. (2005); Canals et al. (2008) are given in the Supplemental Material [cf. Fig. 1a]. This model finds experimental motivation in Yb–based breathing pyrochlores, discussed below.

*Transcription to symmetry–based coordinates. * Our next step is to seek a continuum representation of Eq. (6). To accomplish this, we consider the classical limit where individual components of spin commute, and introduce a set of coarse–grained fields ${\bf m}_{\mathsf{X}}$ which transform as irreducible representations of the lattice symmetry Benton (2014); Yan et al. (2017); Benton et al. (2016). In this basis 111See Supplemental Materials for a more detailed derivation., the Hamiltonian becomes

[TABLE]

where $\mathsf{X}$ runs over irreps of the group $T_{d}$, i.e. $\{\mathsf{A_{2}},\mathsf{E},\mathsf{T_{2}},\mathsf{T_{1+}},\mathsf{T_{1-}}\}$, with the fields $m_{\mathsf{X}}$ and the coefficients $a_{\sf X}$ defined in Table I and Table II of the Supplementary Material.

Before considering the effect of DM interactions, it is helpful to explore how this approach works in the case of a known spin liquid, the HAF on a pyrochlore lattice Anderson (1956); Reimers et al. (1991); Moessner and Chalker (1998a, b); Henley (2005, 2010). Setting

[TABLE]

we find

[TABLE]

It follows that the fields $m_{\mathsf{A_{2}}},m_{\mathsf{E}},m_{\mathsf{T_{2}}},m_{\mathsf{T_{1-}}}$ are all free to fluctuate in the ground state. We can conveniently collect all of these fields in the rank–2 tensor

[TABLE]

where

[TABLE]

[TABLE]

The requirement of the continuity of the fields $m_{\mathsf{X}}$ Benton (2014) imposes the conditions

[TABLE]

We obtain exactly the same constraint if we substitute ${\bf E}^{\sf HAF}$ in Eq. (4), implying that HAF automatically satisfies one of the two constraints defining the R2–U1 spin liquid 222This system of equations can also be viewed as three independent copies of a $U(1)$ gauge theory Henley (2005)..

To convert the HAF into an R2–U1 spin liquid, we need to make the theory symmetric and traceless, and so satisfy Eq. (2). This means eliminating fluctuations of ${\bf E}^{\sf HAF}_{\sf antisym.}$ and ${\bf E}^{\sf HAF}_{\sf trace}$ from the ground state, something which can be accomplished by opening gaps to the fields $m_{\mathsf{T_{2}}}$ and $m_{\mathsf{A_{2}}}$. For the BP model, Eq. (6), this is achieved by any parameter set for which

[TABLE]

In this case, the coefficients $a_{\mathsf{X},\text{A}}$ satisfy the condition

[TABLE]

which implies that only the fields $\mathbf{m}_{\mathsf{E}}$ and $\mathbf{m}_{\mathsf{T_{1-}}}$ enter into the ground state of Eq. (7). Meanwhile, on the B–sublattice, we recover the condition Eq. (9), previously found for the HAF, which imposes the constraint Eq. (13), with the caveat that the fields $m_{\mathsf{A_{2}}}$ and $\mathbf{m}_{\mathsf{T_{2}}}$ can now be set identically equal to zero. When expressed in terms of the remaining tensor field ${\bf E}^{\sf HAF}_{\sf sym}$, this is exactly Eq. (4). It follows that, in this classical limit, an R2–U1 gauge theory, satisfying both Eq. (2) and Eq. (4) emerges as the effective description at the low–energy sector of the BP model, Eq. (6).

It is worth noting that, as in the regular pyrochlore lattice Canals et al. (2008); Chern (2010), DM interaction is only a singular perturbation in the context of the classical ground–state manifold. At finite temperature, classical spin liquids owe their stability to entropy, and a finite value of $D_{A}$ will be needed to stabilise an R2–U1 spin liquid. For exactly the same reason, introducing a finite value of $D_{B}$ does not immediately invalidate the mechanism driving the R2–U1 spin liquid, but will reduce the range of temperatures over which it is observed. We will see that both of these expectations are fulfilled by classical Monte Carlo simulations of Eq. (6), described below

Characteristic signatures of R2–U1 state. We now turn to the question of how the R2–U1 spin liquid can be identified, in both simulation and in experiment. The zero–divergence condition in spin ice, Eq. (1), manifests itself in a pinch–point singularity Henley (2005)

[TABLE]

which is observed in neutron scattering experiments Fennell et al. (2009). In the same way, the constraints associated with an R2–U1 gauge theory, Eq. (2) and Eq. (4), lead to a characteristic singularity in correlations of the tensor field $E_{ij}$ Prem et al. (2018)

[TABLE]

The three–dimensional structure of the correlation $\langle E_{\sf xy}({\bf q})E_{\sf xy}({\bf-q})\rangle$ is illustrated in Fig. 1. In the $[0kl]$ plane, correlations exhibit a conventional 2–fold pinch point, comparable to that found in spin ice [Fig. 1b]. However in the perpendicular $[hk0]$ plane, we observe a 4–fold pinch point (4FPP) [Fig. S1], which unambiguously distinguishes R2–U1 electrodynamics from lower–rank theories Prem et al. (2018).

Comparison with simulation. We can use the existence of this 4FPP as a test for the R2–U1 spin liquid in simulation. We have carried out classical Monte Carlo (MC) simulations of Eq. (6), for the parameter–set

[TABLE]

where the constraints Eq. (2) and Eq. (4) are expected to hold. The resulting correlations of $E_{ij}$, at a temperature $T=2.5\times 10^{-3}\ J_{A}$, are shown in Fig. 1d. For ${\bf q}\to 0$, these are identical to the predictions of Eq. (17), confirming that the model realizes an R2–U1 spin liquid.

*Phase diagram. * The results of simulations for a range of values of $D_{A}$ are collected in Fig. 2. At finite temperature, a finite value of $D_{A}<0$ is required to achieve a crossover from the $U(1)\times U(1)\times U(1)$ spin liquid of the pyrochlore HAF, with 2–fold pinch points, into an R2–U1 spin liquid, with 4FPP. An analytic theory of this crossover, which is controlled by the dimensionless parameter $\eta\sim|D_{A}|/k_{B}T$, is provided in Section VI of the Supplemental Material. Meanwhile, at low temperatures, sufficiently negative values of $D_{A}$ drive a first–order phase transition into an ordered state which involves the characteristic wavevector ${\bf q}={\bf W}$ (i.e. corners of the Brillouin zone) 333A more complex, multiple–${\bf q}$ ground state is hard to rule out categorically, because of the difficulty of thermalising simulations at the lowest temperatures.. In contrast, a finite value of $D_{A}>0$ leads to a continuous phase transition into a state with ${\bf q}=0$, all–out (AIAO) order.

Predictions for neutron scattering. Neutron scattering experiments do not measure correlations of $E_{ij}$ directly, but rather the spin structure factor $S^{\alpha\beta}({\bf q})=\langle S^{\alpha}({\bf q})S^{\beta}(-{\bf q})\rangle$. On general grounds Prem et al. (2018), $S^{\alpha\beta}({\bf q})$ is expected to bear witness to the singularity in Eq. (17). But exactly how 4FPPs would manifest themselves in experiment remains an open question. In Fig. 3 we present simulation results for $S^{\alpha\beta}({\bf q})$ for parameters equivalent to Fig. 3a. We find that the 4FPP is not visible in the structure factor measured by unpolarised neutron scattering [see Supplemental Material]. However the 4FPP can be resolved using polarised neutrons. In this case, it manifests itself in the spin–flip (SF) channel for neutrons polarised perpendicular to the scattering plane Fennell et al. (2009), [Fig. S1].

*Application to materials. * Breathing–pyrochlore magnets were first studied as a tractable limit of the pyrochlore HAF Harris et al. (1991); Canals and Lacroix (1998, 2000); Tsunetsugu (2001), but have since been realised in materials based on both transition–metal Okamoto et al. (2013); Tanaka et al. (2014); Nilsen et al. (2015); Wawrzyńczak et al. (2017); Okamoto et al. (2018) and rare–earth ions Kimura et al. (2014); Haku et al. (2016). Interesting parallels are also found in lacunar spinels Widmann et al. (2016); Jeong et al. (2017). To date, most theoretical work has concentrated on $SU(2)$–invariant models Harris et al. (1991); Canals and Lacroix (1998, 2000); Tsunetsugu (2001); Benton and Shannon (2015); Li et al. (2016); Essafi et al. (2017). However, in the presence of spin–orbit coupling, the symmetry of the lattice permits anisotropic exchange Rau et al. (2016); Haku et al. (2016); Savary et al. (2016); Rau and Gingras (2018). And, with respect to higher–rank gauge theories, a promising line of enquiry are Yb–based materials, where the required form of interactions appear to predominate.

One concrete example is Ba3Yb2Zn5O11 Kimura et al. (2014); Rau et al. (2016); Haku et al. (2016); Rau and Gingras (2018), where A–tetrahedra are estimated to have the coupling parameters $J_{A}\approx 0.57\ \text{meV}$, $D_{A}\approx-0.16\ \text{meV}$, with other interactions negligible. This is exactly the form of interactions needed for an R2–U1 spin liquid, a feature which is expected to be robust Rau and Gingras (2018), since it holds for a wide range of Slater–Koster overlap ratios Slater and Koster (1954). Meanwhile, exchange interactions on the larger B–tetrahedra of Ba3Yb2Zn5O11, while less well understood, appear to be orders of magnitude smaller Rau et al. (2016); Haku et al. (2016). Thus, while it seems plausible that Ba3Yb2Zn5O11 could realise a R2–U1 spin liquid, this may occur at temperatures too low to measure.

The encouraging example of Ba3Yb2Zn5O11 motivates us to consider the possibility of a magnet with similar structure, but smaller B–tetrahedra, such that the interactions on the B–sublattice become non–negligible. For concreteness, we consider a parameter set:

[TABLE]

where we assume that the interactions on the B–sublattice are of the same form as on the A–sublattice, but substantially weaker, $J_{A}/J_{B}=D_{A}/D_{B}\approx 20$. To demonstrate that the R2–U1 physics persists in the presence of finite $D_{B}$ we have used MC simulation to calculate the spin structure factor. Once again, the 4FPP associated with the R2–U1 spin liquid remains clearly visible for a range of temperatures [Fig. 3b]. The same will hold for a more general choice of interactions, as long as the anisotropic part of the exchange on the B–sublattice is sufficiently weak; for $D_{B}\sim D_{A}$, fluctuations are restricted to the local easy plane, and the R2–U1 physics will be lost.

Quantum effects. The theory of an R2–U1 spin liquid presented above is classical, so it is important to ask what might change once quantum effects are taken into account. A useful point of comparison is quantum spin ice (QSI), where quantum fluctuations leads to tunnelling between different spin configurations satisfying the “ice rules” constraint Eq. (1). This tunnelling, which occurs on loops of spins, introduces a fluctuating magnetic field ${\bf B}$, and the result, at $T=0$, is a QSL described by a the deconfined phase of a $U(1)$ quantum lattice gauge theory Hermele et al. (2004); Banerjee et al. (2008); Benton et al. (2012); Savary and Balents (2012); Shannon et al. (2012); Hao et al. (2014); Gingras and McClarty (2014); Kato and Onoda (2015); Chen (2017); Huang et al. (2018). However it is important to note that the temperature scale associated with this QSL is three orders of magnitude smaller than the range of temperatures over which Eq. (1) holds Huang et al. (2018). Moreover, since the $U(1)$ QSL is gapless, any finite temperature immediately restores classical correlations at long length scales Benton et al. (2012). As a consequence, the spin structure factor $S({\bf q})$ continues to be dominated by pinch–point singularities of the form Eq. (16), down to the lowest temperatures studied Kato and Onoda (2015).

The quantum limit of R2–U1 gauge theories has already been studied as a continuum field theory, and is qualitatively very similar to QSI Pretko (2017a, b); Prem et al. (2018). The lowest lying excitations are gapless emergent photons which modify, but do not eliminate, the singular features observed in scattering Benton et al. (2012); Prem et al. (2018). The microscopic study of quantum effects in Eq. (6) lies outside the scope of this Letter. However we anticipate that coherent gauge fluctuations will be confined to an even lower temperature scale than in QSI, by the fact that the magnetic field $B_{ij}$ is an extended object, involving third–order derivatives of $A_{ij}$ Pretko (2017a, b). For this reason the classical theory developed here should prove sufficient to interpret experiments searching for an R2–U1 in a BP material.

Summary and perspectives. In this Letter, we have used a combination of analytic field theory and classical Monte Carlo simulation to show how a rank–2 $U(1)$ [R2–U1] spin liquid, a state described by a higher–rank generalisation of electrodynamics, can arise in a pyrochlore magnet with breathing anisotropy and Dzyaloshinskii–Moriya interactions, Eq. (6) [cf. Fig. S1]. These results provide a concrete starting point for the experimental search for higher–rank gauge theories, and clarify the type of neutron scattering experiment which would be needed to resolve the 4–fold pinch points (4FPP) of a R2–U1 spin liquid [cf. Fig. 3].

This work opens a number of interesting perspectives. On the experimental side, we identify Yb based breathing–pyrochlore materials as potential candidates for a R2–U1 spin liquid state. On the theoretical side, determining the quantum ground state of Eq. (6), should ultimately prove tractable, since breathing anisotropy provides a natural control parameter for both perturbative Canals and Lacroix (1998, 2000) and variational approaches Benton et al. (2018). And, while the model studied here does not correspond to a fracton stabilizer code upon Higgsing Bulmash and Barkeshli (2018); Ma et al. (2018), the parital–confinement mechanism used to eliminate selected components of the tensorial electric field is very versatile, and easily adapted to generate other versions of R2–U1 theory Yan and Benton (shed).

Acknowledgements. The authors acknowledge helpful conversations with Jeffrey Rau and Daniel Khomskii. This work was supported by the Theory of Quantum Matter Unit, Okinawa Institute of Science and Technology Graduate University (OIST). H.Y. is also supported by Japan Society for the Promotion of Science (JSPS) Research Fellowship for Young Scientists. L. J. acknowledges financial support from the French “Agence Nationale de la Recherche” under Grant No. ANR-18-CE30-0011-01, and hospitality from Gakushuin University under Grants-in-Aid for Scientific Research on innovative areas “Topological Materials Science” (No.JP15H05852) from JSPS. The research was also supported in part by the National Science Foundation under Grant No. NSF PHY-1748958, and the KITP program ”Topological Quantum Matter: From Concepts to Realizations”.

I Rank–2 $U(1)$ gauge theory: electrostatics

Here, following Ref. Xu (2006); Pretko (2017a), we derive the relationship between electric, magnetic and gauge fields within the rank–2 $U(1)$ [R2–U1] electrodynamics considered in the main text. This section focuses on the classical electrostatics which is realized in our work. The next section will focus on the quantum dynamics of the theory, which is beyond the scope of this work but nevertheless essential for future developments.

Our starting point is an electric field described by a symmetric, traceless rank–2 tensor,

[TABLE]

Here we do not distinguish super and subscript since we are dealing with spacial indices.

The low energy sector of the electric field has vanishing vector charge, and is traceless,

[TABLE]

Here we keep all indices as subscripts but still use the Einstein summation rule. The proper rank-2 tensor with the proper Gauss law as a classical spin liquid system is achieved in this work.

II Rank–2 $U(1)$ gauge theory: dynamics

A quantum spin liquid requires quantum dynamics in addition to the emergent Gauss law. Broadly speaking, the dynamics play the role of ${\bf B}^{2}$ term in electrodynamics. They are to tunnel different classical spin liquid states between each other, leading to a long-range entangled quantum ground state and gapless photon excitations. This section explains how to derive such terms and also their implication on mobility of electric charges (fractons).

As in conventional electrodynamics, the conjugate of ${\bf E}$ is the rank-two gauge field $\bf A$, which also has to be symmetric to match the degrees of freedom,

[TABLE]

These two conditions determine the form of gauge transformation. Consider a wave-function

[TABLE]

We take a low energy configuration of $\bf E$ obeying Eq. (S2) and construct a symmetrized operator that is identical to zero to act upon the wave-function

[TABLE]

By integration by parts and assuming vanishing boundary terms, we have

[TABLE]

Since $E_{ij}$ conjugates with $A_{ij}$, it generates a transformation of $\bf A$. Thus

[TABLE]

That is, the low energy sector wave-function is invariant under gauge transformation

[TABLE]

Similarly, the traceless condition

[TABLE]

leads to another gauge symmetry

[TABLE]

Finally, the magnetic field is obtained by finding the simplest gauge-invariant quantity. In this case, it has to have three derivatives acting on the gauge field,

[TABLE]

Finally, the Gauss law, the traceless and symmetric conditions of the electric field can be used to derive:

[TABLE]

In this case, a vector charge excitation is fully fractonic, i.e., it cannot move in any direction of the system.

Further details of the phenomenology of R2–U1 phases can be found in Refs. Pretko (2017a, b).

III Derivation of Effective Field Theory

We show how a rank–2 tensor electric field, satisfying the constraint required for R2–U1 electrodynamics [Eqs. (1,2)], can be derived from a breathing pyrochlore lattice model [Eq. (6)]. The pattern of this derivation closely follows Refs. Benton et al. (2016); Benton (2014); Yan et al. (2017).

Our starting point is the breathing pyrochlore lattice with a spin on each of its sites, and nearest neighbor interactions between the spins. “Breathing” means the lattice is bi-partitioned into A- and B-tetrahedra [Fig. (1)], and each type of tetrahedron has its own interactions.

The model that hosts a rank-2 spin liquid has breathing Heisenberg antiferromagnetic interactions on both the A- and B-tetrahedra, and negative Dzyaloshinskii-Moriya (DM) interactions on A-tetrahedra only. The Hamiltonian for the model is

[TABLE]

where $\langle ij\rangle\in\text{A(B)}$ denotes nearest neightbour bonds belonging to the A(B)-tetrahedra. The sites $0,\ 1,\ 2,\ 3$ are at positions relative to the center of an A-tetrahedron

[TABLE]

where $a$ is the length of the unit cell. Vectors $\hat{\bf d}_{ij}$ are bond dependent, defined in accordance with Ref Kotov et al. (2005); Canals et al. (2008); Rau et al. (2016):

[TABLE]

Equivalently, this model can be written in a standard matrix-exchange form for a breathing-pyrochlore lattice model as

[TABLE]

where $\mathcal{J}_{\text{A,ij}}$ is a three-by-three matrix that couples spins on sub-lattice sites $i,j$ whose bond belongs to A-tetrahedra, and $\mathcal{J}_{\text{B}}$ is the coupling matrix for B-tetrahedra. In the case of $D_{B}=0$ that we are interested in, $\mathcal{J}_{\text{B}}$ is identical for any pair of $i,j$,

[TABLE]

Matrices $\mathcal{J}_{\text{A,ij}}$ are bond dependent and related to each other by the lattice symmetry. Their values are

[TABLE]

The spin degrees of freedom on each tetrahedron can be rewritten in terms of fields forming the irreducible representations of the lattice symmetry,

[TABLE]

whose definition can be found in Table 1. They are linear combinations of the spin degrees of freedom, allowing for a fully quadratic Hamiltonian:

[TABLE]

where $\mathsf{X}$ runs over irreps of the group $T_{d}$, i.e. $\{\mathsf{A_{2}},\mathsf{E},\mathsf{T_{2}},\mathsf{T_{1+}},\mathsf{T_{1-}}\}$ as listed in Eq. (S21), and the subscript A,B denotes on which type of tetrahedra they are defined. The coefficients $a_{\mathsf{X}}$ are listed in Table. 2.

For the couplings in Eq. (S15), we have on A-tetrahedra

[TABLE]

and on B-tetrahedra

[TABLE]

For $J_{A},J_{B}>0$ and $D_{A}<0$, these parameters are in order

[TABLE]

which plays the central role of dictating the low energy physics.

The irreducible representation fields are subject to constraints arising from fixed spin length

[TABLE]

for both A- and B-tetrahedra. As a consequence, the low energy sector allows the $m^{2}_{\mathsf{X}}$ corresponding to the smallest $a_{\mathsf{X}}$ to fluctuate, while all other fields have to vanish. This principle applied to our model leads to

• •

On A-tetrahedra, the fields $\mathbf{m}_{\mathsf{E}}$ and $\mathbf{m}_{\mathsf{T_{1-}}}$ can fluctuate;

• •

On A-tetrahedra, the fields ${\bf m}_{\mathsf{T_{1+}}}={\bf m}_{\mathsf{T_{2}}}=\mathbf{0},\ {m}_{\mathsf{A_{2}}}=0$;

• •

On B-tetrahedra, the fields ${m}_{\mathsf{A_{2}}},\;\mathbf{m}_{\mathsf{E}},\;\mathbf{m}_{\mathsf{T_{2}}},\;\mathbf{m}_{\mathsf{T_{1-}}}$ can fluctuate;

• •

On B-tetrahedra,

[TABLE]

Since every spin is shared by an A- and a B-tetrahedron, the fluctuating fields $\mathbf{m}_{\mathsf{E}}$ and $\mathbf{m}_{\mathsf{T_{1-}}}$ on A-tetrahedra must obey additional constraints to respect the the low-energy sector condition on B-tetrahedron imposed by Eq. (S32). Assuming that the fields are varying slowly in space such that the continuous limit can be taken, the constraint Eq. (S32) can be expressed in terms of fields living on A-tetrahedron as

[TABLE]

From this constraint we can build the symmetric, traceless, rank-two magnetic field $E_{ij}$ as

[TABLE]

such that Eq. (S33) becomes

[TABLE]

with symmetric and traceless conditions

[TABLE]

by the definition of $E_{ij}$. Hence a rank-2, traceless, vector charged magnetic field emerges at the low-energy sector from the microscopic model of Eq. (S15).

Equation (S36) constrains the form of correlations functions of $\langle E_{ij}({\bf q})E_{kl}(-{\bf q})\rangle$, in the same spirit as how the two-in-two-out condition constrains the spin-spin correlation of spin ice. It is, however, in a more complicated form. The explicit results for the traceful scalar-charged and vector-charged versions of R2–U1 are discussed in detail in Ref. Prem et al. (2018). The vector-charge field correlation is

[TABLE]

In close analogy, we derive the correlation function of our traceless vector-charged model by deducting the trace,

[TABLE]

which encodes a singularity at ${\bf q}\rightarrow 0$. Different choices of the components $E_{ij}$ and $E_{kl}$ show different patterns. A few representative ones can be found in Figs. S1,S2a.

Fig. S2b,S2c have the four-fold pinch-point (4FPP) singularity, which differentiates the rank-2 gauge theories uniquely from the conventional $U(1)$ gauge theory. It is the key signature to be looked for in experiments.

IV Alternative forms of R2–U1 spin liquid

In the main text we have shown how a classical spin liquid described by a symmetric, traceless R2–U1 gauge theory, descends from the $U(1)\times U(1)\times U(1)$ spin liquid found in the classical Heisenberg Antiferromagnet (HAF) on a pyrochlore lattice. The approach taken is very versatile, and by tuning the Hamiltonian, one can also obtain other forms of R2–U1 spin liquid.

Notice that the diagonal and off-diagonal components of ${\bf E}^{\sf HAF}_{\sf sym.}$ come from different irreps ${\bf m}_{\mathsf{E}}$ and ${\bf m}_{\mathsf{T_{1-}}}$. In the most general Hamiltonian (Eq.S22), these two irreps have their individually tunable coefficients $a_{\mathsf{E}}$ and $a_{\mathsf{T_{1-}}}$. So the symmetric part of ${\bf E}^{\sf HAF}$ can be decomposed into three components

[TABLE]

and each component can be individually tuned to be active or suppressed by choosing the proper Hamiltonian. The vector-charged Gauss law is unaffected. This allows us to build a variety of rank-2 U(1) gauge theories, including a “hollow” version with vanishing diagonal terms Ma et al. (2018) .

V Predictions for Neutron Scattering

The 4FPP is a unique pattern that differentiates the R2–U1 from vector $U(1)$ gauge theory, which only has the conventional two-fold pinch points. The 4FPPs are most unambiguously presented in the correlation function of the irrep fields as discussed in the main text. These correlation functions are, however, not directly accessible in experiment.

In magnetism, neutron scattering is widely used to measure the spin-spin correlation of the form

[TABLE]

where $\alpha,\beta=x,y,z$ are spin-component indices and $i,j=0,1,2,3$ are sub-lattice site indices.

Furthermore, with neutrons polarized in direction of unit vector $\hat{\bf v}$ perpendicular to the scattering plane, one can measure the spin-flip (SF) channel neutron scattering defined by

[TABLE]

where

[TABLE]

One can also measure the non-spin-flip (NSF) channel defined by

[TABLE]

Here we show the spin structure factor of the $[h0k]$ and $[hhk]$ plane, with zoomed-in view of the 4FPPs.

VI Temperature evolution of 4–fold pinch point into

a conventional 2–fold pinch point

The cross–over from the HAF phase to the R2–U1 phase, with decreasing temperature, is manifested in the structure factor as a cross–over between a conventional, 2–fold, pinch point and the 4–fold pinch point (4FPP) characteristic of an R2–U1 gauge theory Prem et al. (2018). To provide further quantitative details of this cross–over, here we present an analysis of correlations based on a coarse–grained field theory.

As shown above, the $U(1)\times U(1)\times U(1)$ spin liquid of the Heisenberg model can be described in terms of a non-symmetric matrix ${\bf E}$, with a Gauss’ law applied to each column. The traceless R2-U(1) spin liquid is described in the same way, but with the constraint that ${\bf E}$ be symmetric and traceless.

The following effective theory captures both cases:

[TABLE]

where the integral in Eq. (S44) is taken independently over all components of the matrix ${\bf E}$ and $\beta=1/T$ is the inverse temperature.

The limit $\delta\to\infty,\eta\to 0$, captures the Heisenberg model spin liquid, with Gauss’ law enforced on every column of ${\bf E}$ and no correlations between columns. The limit $\delta\to\infty,\eta\to\infty$ captures the R2-U(1) spin liquid, with ${\bf E}$ forced to be symmetric and traceless and still obeying Gauss’ law for each column.

In terms of the parameters of the microscopic model:

[TABLE]

since $\delta$ is the coefficient enforcing the Gauss’ law (which is generated by $J$) and $\eta$ is the coefficient enforcing the symmetric and traceless conditions (which are generated by $D_{A}$).

We can study the crossover from the Heisenberg to R2-U(1) spin liquids by first taking $\delta\to\infty$ and observing the behavior as a function of $\eta$. The crossover can be illustrated by calculating the correlation function

[TABLE]

which should have a 2-fold pinch point in the Heisenberg limit and a 4-fold pinch point in the R2-U(1) limit.

Calculating the correlation function from Eqs. (S44)-(S45) and taking the limit $\delta\to\infty$. gives us

[TABLE]

The $\eta\to 0$ limit of Eq. (S47) gives the Heisenberg limit of the correlation function (a 2-fold pinch point)

[TABLE]

Whereas, the $\eta\to\infty$ limit gives the R2-U(1) correlation function (a 4-fold pinch point):

[TABLE]

Since $\eta\sim\beta|D_{A}|$, we can see from Eq. (S47) that for a fixed ${\bf q}$ at any, fixed, finite temperature the evolution as a function of $D_{A}$, will be smooth, although at small temperatures, a small change in $D_{A}$ will give a large change in $\eta$ and hence the correlation function. Only at $T=0$ does the correlation function behave in a singular fashion as a function of $D_{A}$, but this is not surprising, since the Heisenberg limit has a highly degenerate ground state.

The progression of the correlation function as $\eta$ is increased is given in Fig (S5). This progression can either be seen as decreasing the temperature at fixed, small, $|D_{A}|$ or as increasing $|D_{A}|$ at fixed small temperature.

VII Monte Carlo simulations

Monte Carlo simulations are performed on systems of classical O(3) spins with $16L^{3}$ sites, where $L^{3}$ is the number of cubic unit cells. The spin length is $|S|=1/2$. To decorrelate the system, we use jointly the heatbath method, over-relaxation and parallel tempering. Thermalization is made in two steps: first a slow annealing from high temperature to the temperature of measurement $T$ during $t_{e}$ Monte Carlo steps (MCs) followed by $t_{e}$ MCS at temperature $T$. After thermalization, measurements are done every 10 MCs during $t_{m}=10~{}t_{e}$ MCs. All structure factors have been computed from simulations with $L=30$ and $t_{m}=5\times 10^{5}$ MCs. The phase diagram of Fig. 5 has been computed from simulations with $L=8$ and $t_{m}=10^{7}$ MCs.