Schrieffer-Wolff transformation) lies entirely in the ice manifold, then using
Eq. (19), we get

$\matrixelement{\Psi}{\hat{H}}{\Psi}=\left(\bra{\Psi}\hat{U}_{S}^{\dagger}\right)\hat{H}_{\text{eff}}\left(\hat{U}_{S}\ket{\Psi}\right).$
(30)

For the ground state, $\ket{\Psi_{g}}$ of the full Hamiltonian $\hat{H}$,
$\hat{U}_{S}\ket{\Psi_{g}}$ lies entirely in the ice manifold. Thus, we pick
an ansatz wavefunction for $\hat{U}_{S}\ket{\Psi}$ that also lies entirely in
the ice manifold and compute the expectation value of $\hat{H}_{\text{eff}}$
in our ansatz state to get the energy. Before describing our ansatz states in
Sec. III.1.3, we first consider the limit $\Omega=0$ in the next section.

####  III.1.2 Classical ground state of the long-range Hamiltonian

Here, we will find the ground state selected by long-range interactions in the
limit $\Omega=0$ where there are no quantum fluctuations. The Hamiltonian is
$\hat{H}_{\text{cl}}=\hat{H}_{0}+\hat{H}_{LR}$. We find the ground state by
going to the Fourier space. Since the pyrochlore lattice is an FCC lattice
with a four-site basis, we use the notation $\hat{S}^{z}_{\mathbf{r},\mu}$ for
spins where $\mathbf{r}$ is an FCC lattice vector and $\mu\in\\{0,1,2,3\\}$
labels the sites within the basis. The spin $\hat{S}^{z}_{\mathbf{r},\mu}$ is
physically located at $\mathbf{r}+\mathbf{e}_{\mu}/2$ where $\mathbf{e}_{\mu}$
are the vectors joining a diamond A site to a neighboring diamond B site. (See
Fig. 1(b) for the precise definition.) As we are considering the classical
limit in this section, we drop hats on quantities which would otherwise be
operators. The Fourier transform of $S_{\mathbf{r},\mu}^{z}$ is

$S^{z}_{\mathbf{r},\mu}=\frac{1}{\sqrt{N_{\text{u.c.}}}}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}S^{z}_{\mathbf{k},\mu},$
(31)

where $N_{\text{u.c.}}$ is the number of FCC unit cells. Substituting this in
$H_{\text{cl}}$, we get

$H_{\text{cl}}=\sum_{\mu,\nu,\mathbf{k}}V_{\mu\nu,\mathbf{k}}S^{z}_{\mathbf{k},\mu}S^{z}_{-\mathbf{k},\nu},$
(32)

where $\mathbf{k}$ is a vector in the Brillouin zone of the FCC lattice and
$V_{\mu\nu,\mathbf{k}}$ is the Fourier transform of the van der Waals
potential:

$V_{\mu\nu,\mathbf{k}}=\frac{V}{2}\sum_{\mathbf{r}}e^{i\mathbf{k}\cdot\mathbf{r}}\left(\frac{a}{|\mathbf{r}+(\mathbf{e}_{\mu}-\mathbf{e}_{\nu})/2|}\right)^{6}.$
(33)

Diagonalizing the matrix $V_{\mu\nu,\mathbf{k}}$ for each $\mathbf{k}$ gives

$H_{\text{cl}}=\sum_{\mu,\mathbf{k}}\varepsilon_{\mathbf{k},\mu}|S^{{}^{\prime}z}_{\mathbf{k},\mu}|^{2},$
(34)

where $S^{{}^{\prime}z}_{\mathbf{k},\mu}$ is related to
$S^{z}_{\mathbf{k},\nu}$ through a multiplication by a unitary matrix
$U_{\mu\nu,\mathbf{k}}$ which diagonalizes $V_{\mu\nu,\mathbf{k}}$:
$S^{{}^{\prime}z}_{\mathbf{k},\mu}=\sum_{\nu}U_{\mu\nu,\mathbf{k}}S^{z}_{\mathbf{k},\nu}$.
Recall that $S^{z}_{\mathbf{r},\mu}$ is either $+1/2$ or $-1/2$. This imposes
the following constraint:

$\sum_{\mathbf{k},\mu}|S^{{}^{\prime}z}_{\mathbf{k},\mu}|^{2}=\sum_{\mathbf{k},\mu}|S^{z}_{\mathbf{k},\mu}|^{2}=\sum_{\mathbf{r},\mu}\left(S^{z}_{\mathbf{r},\mu}\right)^{2}=N_{\text{u.c.}}.$
(35)

From Eq. (34), the energy can be interpreted as a weighted sum of
$\varepsilon_{\mathbf{k},\mu}$ with the corresponding weights being
$|S^{{}^{\prime}z}_{\mathbf{k},\mu}|^{2}$. Because of the constraint in Eq.
(35), the energy is minimized by having the full weight on the smallest
$\varepsilon_{\mathbf{k},\mu}$ and no weight on the rest of the
$\varepsilon_{\mathbf{k},\mu}$. This holds provided that such a configuration
of $S^{{}^{\prime}z}_{\mathbf{k},\mu}$ in the momentum space corresponds to
some configuration in the real space where $S^{z}_{\mathbf{r},\mu}$ are $\pm
1/2$.

Calculating the Fourier transform of the long-range potential, Eq. (33), and
its eigenvalues $\varepsilon_{\mathbf{k},\mu}$, we find that the minimum of
$\varepsilon_{\mathbf{k},\mu}$ occurs for $\mathbf{k}=\bm{0}$ and is triply
degenerate. In particular,

$V_{\mu\nu,\mathbf{k}=\bm{0}}=\begin{pmatrix}v_{1}&v_{2}&v_{2}&v_{2}\\\
v_{2}&v_{1}&v_{2}&v_{2}\\\ v_{2}&v_{2}&v_{1}&v_{2}\\\
v_{2}&v_{2}&v_{2}&v_{1}\\\ \end{pmatrix},$ (36)

where $v_{1}=0.113V$ and $v_{2}=1.12V$. Its eigenvalues are
$\varepsilon_{\bm{0},0}=3.46V$ and
$\varepsilon_{\bm{0},1}=\varepsilon_{\bm{0},2}=\varepsilon_{\bm{0},3}=-1.004V$.
The unitary that diagonalizes the above matrix also relates
$S^{{}^{\prime}z}_{\bm{0},\mu}$ to $S^{z}_{\bm{0},\nu}$ as

$\begin{pmatrix}S^{{}^{\prime}z}_{\bm{0},0}\\\ S^{{}^{\prime}z}_{\bm{0},1}\\\
S^{{}^{\prime}z}_{\bm{0},2}\\\ S^{{}^{\prime}z}_{\bm{0},3}\\\
\end{pmatrix}=\frac{1}{2}\begin{pmatrix}1&1&1&1\\\ 1&1&-1&-1\\\ 1&-1&1&-1\\\
1&-1&-1&1\\\ \end{pmatrix}\begin{pmatrix}S^{z}_{\bm{0},0}\\\
S^{z}_{\bm{0},1}\\\ S^{z}_{\bm{0},2}\\\ S^{z}_{\bm{0},3}\\\ \end{pmatrix}.$
(37)

Since $\varepsilon_{\bm{0},1},\varepsilon_{\bm{0},2}$ and
$\varepsilon_{\bm{0},3}$ are the minimum eigenvalues, the energy is minimized
by having all the weight distributed between
$S^{{}^{\prime}z}_{\bm{0},1},S^{{}^{\prime}z}_{\bm{0},2}$ and
$S^{{}^{\prime}z}_{\bm{0},3}$ and no weight on the remaining
$S^{{}^{\prime}z}_{\mathbf{k},\mu}$, that is,
$S^{{}^{\prime}z}_{\mathbf{k}\neq\bm{0},\mu}=0$ and
$S^{{}^{\prime}z}_{\bm{0},0}=0$. There indeed exist states satisfying these
two conditions. The first condition,
$S^{{}^{\prime}z}_{\mathbf{k}\neq\bm{0},\mu}=0$, implies that the ground state
is a $\mathbf{k}=\bm{0}$ state, while the second condition,
$S^{{}^{\prime}z}_{\bm{0},0}=0$, implies that the ground state satisfies the
ice rule (so that the total spin, which is $S^{{}^{\prime}z}_{\bm{0},0}$ is
0), see Eq. (37). There are six such states, and we refer to them as the “ice
ferromagnet” or “ice FM” states. One of these is shown in Fig. 7.

Figure 7: An ice ferromagnet state. It is an ice rule obeying state (i.e.,
$n_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}}=2$
on every tetrahedron) with $\mathbf{k}=\bm{0}$. All the up-pointing tetrahedra
are copies of each other. The same is true for the down-pointing tetrahedra.
There are six (${}^{4}C_{2}$) such states, and together they make up the
ground subspace of $\hat{H}_{\text{cl}}$.

####  III.1.3 Ansatz wavefunctions for the ordered state and for the quantum
spin liquid

We now assume that, as one increases $\Omega$ starting from $\Omega=0$, the
ground state remains adiabatically connected to the ice ferromagnet derived in
the previous section till the point where it undergoes the putative phase
transition to the QSL. Therefore, our ansatz for the ordered state is

$\ket{\Psi_{\text{ord}}}=\hat{U}^{\dagger}_{S}\ket{\Psi_{\text{IFM}}},$ (38)

where $\ket{\Psi_{\text{IFM}}}$, a product state in the $\hat{S}^{z}$ basis,
is the $\mathbf{k}=\bm{0}$ ice ferromagnet defined Sec. III.1.2. This
configuration is given by
$S^{z}_{\mathbf{r},\mu}=\frac{1}{2}\varepsilon_{\mu}$ (independent of
$\mathbf{r}$), where
$\left(\varepsilon_{0},\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}\right)\equiv\left(1,1,-1,-1\right)$.
We note that there are six such choices for $\varepsilon_{\mu}$ that satisfy
the ice rule. We pick one such choice, but our calculations are not sensitive
to which one we pick. $\ket{\Psi_{\text{IFM}}}$ lives entirely in the ice
manifold. Left-multiplication by $\hat{U}^{\dagger}_{S}$ takes it back to the
original Hilbert space with ice rule violations.

Our ansatz wave function for the spin liquid state is

$\ket{\Psi_{\text{QSL}}}=\hat{U}^{\dagger}_{S}\ket{\Psi_{\text{RK}}},$ (39)

where $\ket{\Psi_{\text{RK}}}$ is a uniform superposition of all dimer
coverings [65] of the diamond lattice (with
$n_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}}=2$).
$\ket{\Psi_{\text{RK}}}$ lives in the ice manifold. Like before, we left-
multiply it by $\hat{U}^{\dagger}_{S}$ to take it back to the original Hilbert
space. The justification for our choice is the following.
$\ket{\Psi_{\text{RK}}}$ is the ground state of the dimer model at the RK
point [see Eq. (10)]. When the RK potential is zero, $\ket{\Psi_{\text{RK}}}$
has an energy expectation value of
$-4N_{\text{u.c.}}J_{\text{ring}}\overline{n}_{\text{flip}}$, where
$\overline{n}_{\text{flip}}$ is the average fraction of flippable hexagons in
the RK wavefunction. We find numerically that
$\overline{n}_{\text{flip}}=0.1757$ (also calculated in Ref. [48]). Therefore,
the energy of $\ket{\Psi_{\text{RK}}}$ is
$-0.7028J_{\text{ring}}N_{\text{u.c.}}$ which is not too far from the ground
state energy of the dimer model (12) found in Ref. [51] to be
$-0.756J_{\text{ring}}N_{\text{u.c.}}$. Even though $\ket{\Psi_{\text{RK}}}$
has slightly higher energy, it has the advantage of being simpler to sample by
classical Monte Carlo. This explains our choice.

For comparison, we will also calculate the energy of a different ordered state
$\ket{\Psi_{\text{ord}}^{\prime}}=\hat{U}_{S}^{\dagger}\ket{\Psi_{\text{IAFM}}}$
that we call an ice antiferromagnet. Here $\ket{\Psi_{\text{IAFM}}}$ is an ice
rule obeying state with ordering wave vector
$\mathbf{k}=\pi\left(\mathbf{b}_{1}+\mathbf{b}_{2}\right)$, where
$\mathbf{b}_{1}$, $\mathbf{b}_{2}$ and $\mathbf{b}_{3}$ are primitive
reciprocal lattice vectors of the FCC lattice satisfying
$\mathbf{a}_{i}\cdot\mathbf{b}_{j}=\delta_{ij}$. This state is known elsewhere
in literature as the $2\pi(001)$ state (this nomenclature uses an enlarged
cubic unit cell of the FCC lattice) [78, 79, 76].

####  III.1.4 Numerical results—energy expectation values and phase diagram

We now describe our computation of the expectation value of
$\hat{H}_{\text{eff}}$ [see Eq. (22)] in $\ket{\Psi_{\text{RK}}}$,
$\ket{\Psi_{\text{IFM}}}$ and in $\ket{\Psi_{\text{IAFM}}}$. While the
expectation value in $\ket{\Psi_{\text{IFM}}}$ and $\ket{\Psi_{\text{IAFM}}}$
can be computed straightforwardly, the expectation value in
$\ket{\Psi_{\text{RK}}}$ requires classical Monte Carlo sampling. We use a
system with $8\times 8\times 8$ unit cells (i.e., containing $2048$ pyrochlore
sites) with periodic boundary conditions in the $\mathbf{a}_{1}$,
$\mathbf{a}_{2}$, and $\mathbf{a}_{3}$ directions. We restrict our sampling to
sectors in which the total electric flux piercing through any 2D torus cross-
section (as defined in Sec. IV B of Ref. [48]) is 0. Our sampling is done
using loop moves as described in Refs. [48, 78, 79] – in each Monte Carlo run,
we perform $512\times 500,000$ loop moves. We calculate
$\bar{n}_{\text{flip}}$, $H_{\text{LR}}$, $W^{(2)}_{\text{LR}}$ and
$L^{(2)}_{\text{LR}}$ after every 512 loop moves, i.e., we take 500,000 data
points. We calculate $M^{(2)}_{\text{LR}}$, $W^{(3)}_{\text{LR}}$, and
$W^{(4)}_{\text{LR}}$ after every $512\times 10,000$ loop moves, i.e. we take
50 data points. We repeat this procedure for 9 independent runs in order to
calculate the uncertainties. Our results are summarized in Table 1.

Operator | $\ket{\Psi_{\text{RK}}}$ | $\ket{\Psi_{\text{IFM}}}$ | $\ket{\Psi_{\text{IAFM}}}$
---|---|---|---
$\hat{R}$ | $0.70288(4)N_{\text{u.c.}}$ | $0$ | $0$
$\hat{H}_{\text{LR}}$ | $2.6037(1)\times 10^{-2}N_{\text{u.c.}}$ | $-0.4002\times 10^{-2}N_{\text{u.c.}}$ | $3.8722\times 10^{-2}N_{\text{u.c.}}$
$\hat{W}^{(2)}_{\text{LR}}$ | $1.11778(1)\times 10^{-3}N_{\text{u.c.}}$ | $0.01642\times 10^{-3}N_{\text{u.c.}}$ | $1.4994\times 10^{-3}N_{\text{u.c.}}$
$\hat{L}^{(2)}_{\text{LR}}$ | $-2.7467(3)\times 10^{-4}N_{\text{u.c.}}$ | $-0.0829\times 10^{-4}N_{\text{u.c.}}$ | $-7.5662\times 10^{-4}N_{\text{u.c.}}$
$\hat{M}^{(2)}_{\text{LR}}$ | $2.96(3)\times 10^{-3}N_{\text{u.c.}}$ | $0.073\times 10^{-3}N_{\text{u.c.}}$ | $6.66\times 10^{-3}N_{\text{u.c.}}$
$\hat{W}^{(3)}_{\text{LR}}$ | $5.25(4)\times 10^{-5}N_{\text{u.c.}}$ | $-0.00665\times 10^{-5}N_{\text{u.c.}}$ | $5.81\times 10^{-5}N_{\text{u.c.}}$
$\hat{W}^{(4)}_{\text{LR}}$ | $-3.57(2)\times 10^{-6}N_{\text{u.c.}}$ | $-0.0309\times 10^{-6}N_{\text{u.c.}}$ | $-5.35\times 10^{-6}N_{\text{u.c.}}$
Table 1: The expectation values of the operators in the left column in ansatz
wavefunctions $\ket{\Psi_{\text{RK}}}$, $\ket{\Psi_{\text{IFM}}}$ and
$\ket{\Psi_{\text{IAFM}}}$ respectively. The operator $\hat{R}$ is defined as
$\hat{R}=\sum_{\hexagon}\ket{\raisebox{-2.84544pt}{ \leavevmode\hbox
to9.79pt{\vbox to10.94pt{\pgfpicture\makeatletter\hbox{\hskip
4.89597pt\lower-5.46773pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}{}{{}}{} {}{} {}{} {}{} {}{}
{}{}
{}{}{}\pgfsys@moveto{3.84448pt}{1.8529pt}\pgfsys@lineto{0.0pt}{4.26773pt}\pgfsys@lineto{-3.69595pt}{2.13388pt}\pgfsys@lineto{-3.69597pt}{-2.13387pt}\pgfsys@lineto{0.0pt}{-4.26773pt}\pgfsys@lineto{3.69595pt}{-2.13388pt}\pgfsys@lineto{3.69597pt}{2.13387pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}{}\pgfsys@moveto{3.69597pt}{-2.13387pt}\pgfsys@moveto{4.69597pt}{-2.13387pt}\pgfsys@curveto{4.69597pt}{-1.58157pt}{4.24826pt}{-1.13387pt}{3.69597pt}{-1.13387pt}\pgfsys@curveto{3.14368pt}{-1.13387pt}{2.69597pt}{-1.58157pt}{2.69597pt}{-2.13387pt}\pgfsys@curveto{2.69597pt}{-2.68616pt}{3.14368pt}{-3.13387pt}{3.69597pt}{-3.13387pt}\pgfsys@curveto{4.24826pt}{-3.13387pt}{4.69597pt}{-2.68616pt}{4.69597pt}{-2.13387pt}\pgfsys@closepath\pgfsys@moveto{3.69597pt}{-2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}{}\pgfsys@moveto{0.0pt}{4.26773pt}\pgfsys@moveto{1.0pt}{4.26773pt}\pgfsys@curveto{1.0pt}{4.82002pt}{0.55229pt}{5.26773pt}{0.0pt}{5.26773pt}\pgfsys@curveto{-0.55229pt}{5.26773pt}{-1.0pt}{4.82002pt}{-1.0pt}{4.26773pt}\pgfsys@curveto{-1.0pt}{3.71544pt}{-0.55229pt}{3.26773pt}{0.0pt}{3.26773pt}\pgfsys@curveto{0.55229pt}{3.26773pt}{1.0pt}{3.71544pt}{1.0pt}{4.26773pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{4.26773pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}{}\pgfsys@moveto{-3.69597pt}{-2.13387pt}\pgfsys@moveto{-2.69597pt}{-2.13387pt}\pgfsys@curveto{-2.69597pt}{-1.58157pt}{-3.14368pt}{-1.13387pt}{-3.69597pt}{-1.13387pt}\pgfsys@curveto{-4.24826pt}{-1.13387pt}{-4.69597pt}{-1.58157pt}{-4.69597pt}{-2.13387pt}\pgfsys@curveto{-4.69597pt}{-2.68616pt}{-4.24826pt}{-3.13387pt}{-3.69597pt}{-3.13387pt}\pgfsys@curveto{-3.14368pt}{-3.13387pt}{-2.69597pt}{-2.68616pt}{-2.69597pt}{-2.13387pt}\pgfsys@closepath\pgfsys@moveto{-3.69597pt}{-2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{3.69597pt}{2.13387pt}\pgfsys@moveto{4.69597pt}{2.13387pt}\pgfsys@curveto{4.69597pt}{2.68616pt}{4.24826pt}{3.13387pt}{3.69597pt}{3.13387pt}\pgfsys@curveto{3.14368pt}{3.13387pt}{2.69597pt}{2.68616pt}{2.69597pt}{2.13387pt}\pgfsys@curveto{2.69597pt}{1.58157pt}{3.14368pt}{1.13387pt}{3.69597pt}{1.13387pt}\pgfsys@curveto{4.24826pt}{1.13387pt}{4.69597pt}{1.58157pt}{4.69597pt}{2.13387pt}\pgfsys@closepath\pgfsys@moveto{3.69597pt}{2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{-3.69597pt}{2.13387pt}\pgfsys@moveto{-2.69597pt}{2.13387pt}\pgfsys@curveto{-2.69597pt}{2.68616pt}{-3.14368pt}{3.13387pt}{-3.69597pt}{3.13387pt}\pgfsys@curveto{-4.24826pt}{3.13387pt}{-4.69597pt}{2.68616pt}{-4.69597pt}{2.13387pt}\pgfsys@curveto{-4.69597pt}{1.58157pt}{-4.24826pt}{1.13387pt}{-3.69597pt}{1.13387pt}\pgfsys@curveto{-3.14368pt}{1.13387pt}{-2.69597pt}{1.58157pt}{-2.69597pt}{2.13387pt}\pgfsys@closepath\pgfsys@moveto{-3.69597pt}{2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{-4.26773pt}\pgfsys@moveto{1.0pt}{-4.26773pt}\pgfsys@curveto{1.0pt}{-3.71544pt}{0.55229pt}{-3.26773pt}{0.0pt}{-3.26773pt}\pgfsys@curveto{-0.55229pt}{-3.26773pt}{-1.0pt}{-3.71544pt}{-1.0pt}{-4.26773pt}\pgfsys@curveto{-1.0pt}{-4.82002pt}{-0.55229pt}{-5.26773pt}{0.0pt}{-5.26773pt}\pgfsys@curveto{0.55229pt}{-5.26773pt}{1.0pt}{-4.82002pt}{1.0pt}{-4.26773pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{-4.26773pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>\leavevmode\hbox to9.79pt{\vbox to10.94pt{\pgfpicture\makeatletter\hbox{\hskip
4.89597pt\lower-5.46773pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}{}{{}}{} {}{} {}{} {}{} {}{}
{}{}
{}{}{}\pgfsys@moveto{3.69597pt}{2.13387pt}\pgfsys@lineto{0.0pt}{4.26773pt}\pgfsys@lineto{-3.69595pt}{2.13388pt}\pgfsys@lineto{-3.69597pt}{-2.13387pt}\pgfsys@lineto{0.0pt}{-4.26773pt}\pgfsys@lineto{3.69595pt}{-2.13388pt}\pgfsys@lineto{3.69597pt}{2.13387pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}{}\pgfsys@moveto{3.69597pt}{2.13387pt}\pgfsys@moveto{4.69597pt}{2.13387pt}\pgfsys@curveto{4.69597pt}{2.68616pt}{4.24826pt}{3.13387pt}{3.69597pt}{3.13387pt}\pgfsys@curveto{3.14368pt}{3.13387pt}{2.69597pt}{2.68616pt}{2.69597pt}{2.13387pt}\pgfsys@curveto{2.69597pt}{1.58157pt}{3.14368pt}{1.13387pt}{3.69597pt}{1.13387pt}\pgfsys@curveto{4.24826pt}{1.13387pt}{4.69597pt}{1.58157pt}{4.69597pt}{2.13387pt}\pgfsys@closepath\pgfsys@moveto{3.69597pt}{2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}{}\pgfsys@moveto{-3.69597pt}{2.13387pt}\pgfsys@moveto{-2.69597pt}{2.13387pt}\pgfsys@curveto{-2.69597pt}{2.68616pt}{-3.14368pt}{3.13387pt}{-3.69597pt}{3.13387pt}\pgfsys@curveto{-4.24826pt}{3.13387pt}{-4.69597pt}{2.68616pt}{-4.69597pt}{2.13387pt}\pgfsys@curveto{-4.69597pt}{1.58157pt}{-4.24826pt}{1.13387pt}{-3.69597pt}{1.13387pt}\pgfsys@curveto{-3.14368pt}{1.13387pt}{-2.69597pt}{1.58157pt}{-2.69597pt}{2.13387pt}\pgfsys@closepath\pgfsys@moveto{-3.69597pt}{2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}{}\pgfsys@moveto{0.0pt}{-4.26773pt}\pgfsys@moveto{1.0pt}{-4.26773pt}\pgfsys@curveto{1.0pt}{-3.71544pt}{0.55229pt}{-3.26773pt}{0.0pt}{-3.26773pt}\pgfsys@curveto{-0.55229pt}{-3.26773pt}{-1.0pt}{-3.71544pt}{-1.0pt}{-4.26773pt}\pgfsys@curveto{-1.0pt}{-4.82002pt}{-0.55229pt}{-5.26773pt}{0.0pt}{-5.26773pt}\pgfsys@curveto{0.55229pt}{-5.26773pt}{1.0pt}{-4.82002pt}{1.0pt}{-4.26773pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{-4.26773pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{3.69597pt}{-2.13387pt}\pgfsys@moveto{4.69597pt}{-2.13387pt}\pgfsys@curveto{4.69597pt}{-1.58157pt}{4.24826pt}{-1.13387pt}{3.69597pt}{-1.13387pt}\pgfsys@curveto{3.14368pt}{-1.13387pt}{2.69597pt}{-1.58157pt}{2.69597pt}{-2.13387pt}\pgfsys@curveto{2.69597pt}{-2.68616pt}{3.14368pt}{-3.13387pt}{3.69597pt}{-3.13387pt}\pgfsys@curveto{4.24826pt}{-3.13387pt}{4.69597pt}{-2.68616pt}{4.69597pt}{-2.13387pt}\pgfsys@closepath\pgfsys@moveto{3.69597pt}{-2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{4.26773pt}\pgfsys@moveto{1.0pt}{4.26773pt}\pgfsys@curveto{1.0pt}{4.82002pt}{0.55229pt}{5.26773pt}{0.0pt}{5.26773pt}\pgfsys@curveto{-0.55229pt}{5.26773pt}{-1.0pt}{4.82002pt}{-1.0pt}{4.26773pt}\pgfsys@curveto{-1.0pt}{3.71544pt}{-0.55229pt}{3.26773pt}{0.0pt}{3.26773pt}\pgfsys@curveto{0.55229pt}{3.26773pt}{1.0pt}{3.71544pt}{1.0pt}{4.26773pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{4.26773pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}{}\pgfsys@moveto{-3.69597pt}{-2.13387pt}\pgfsys@moveto{-2.69597pt}{-2.13387pt}\pgfsys@curveto{-2.69597pt}{-1.58157pt}{-3.14368pt}{-1.13387pt}{-3.69597pt}{-1.13387pt}\pgfsys@curveto{-4.24826pt}{-1.13387pt}{-4.69597pt}{-1.58157pt}{-4.69597pt}{-2.13387pt}\pgfsys@curveto{-4.69597pt}{-2.68616pt}{-4.24826pt}{-3.13387pt}{-3.69597pt}{-3.13387pt}\pgfsys@curveto{-3.14368pt}{-3.13387pt}{-2.69597pt}{-2.68616pt}{-2.69597pt}{-2.13387pt}\pgfsys@closepath\pgfsys@moveto{-3.69597pt}{-2.13387pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>In the
RK wavefunction,
$\matrixelement{\Psi_{\text{RK}}}{\hat{R}}{\Psi_{\text{RK}}}=4\bar{n}_{\text{flip}}N_{\text{u.c.}}$.
To calculate expectation values in $\ket{\Psi_{\text{RK}}}$, we have used
classical Monte Carlo sampling.

With these values at hand, we calculate the expectation value of
$\hat{H}_{\text{eff}}$ using Eq. (22) in $\ket{\Psi_{\text{RK}}}$,
$\ket{\Psi_{\text{IFM}}}$, and $\ket{\Psi_{\text{IAFM}}}$, and the result is
plotted in Fig. 8. As we turn on $\Omega$, the transition point $\Omega_{C}$
can be determined within our approximation as the $\Omega$ for which the
energy of the ice ferromagnet becomes higher than that of the RK wavefunction,
as calculated using Eq. (22). We find

$\Omega_{C}=0.43927(1)V.$ (40)

Figure 8: $\langle\hat{H}_{\text{eff}}\rangle$ in $\ket{\Psi_{\text{RK}}}$,
$\ket{\Psi_{\text{IFM}}}$, and $\ket{\Psi_{\text{IAFM}}}$ calculated by
inserting the values in Table 1 in Eq. (22).

There is an important question on whether our use of perturbation theory is
justified. First, we argue that treating $\hat{H}_{\text{LR}}$ perturbatively
is justified. $\left\\{\hat{H}_{\text{LR}}\right\\}$,
$\left\\{\hat{W}^{(2)}_{\text{LR}},\hat{L}^{(2)}_{\text{LR}},\hat{M}^{(2)}_{\text{LR}}\right\\}$,
$\left\\{\hat{W}^{(3)}_{\text{LR}}\right\\}$, and
$\left\\{\hat{W}^{(4)}_{\text{LR}}\right\\}$ are sets of operators that are
first, second, third, and fourth order respectively in $\hat{H}_{\text{LR}}$.
As we can see from Table 1, the expectation values of these operators in
$\ket{\Psi}_{\text{RK}}$ drops by an order of magnitude each time one goes one
order higher in $\hat{H}_{\text{LR}}$. Next, is perturbation theory in
$\hat{H}_{\Omega}$ justified, given that our calculated $\Omega_{C}$ is
outside the $\Omega\ll V$ regime? We observe that the leading contribution to
$J_{\text{ring}}$ that we dropped,
$\frac{33833}{2592}\frac{(\Omega_{C})^{8}}{V^{7}}=0.018V$ [74], is smaller
than the one we kept, $\frac{63}{16}\frac{(\Omega_{C})^{6}}{V^{5}}=0.028V$. If
we had kept higher order contributions to $J_{\text{ring}}$, it would only
decrease the energy of the QSL relative to the ice ferromagnet and ice
antiferromagnet. Further, the energy of the QSL that we present is a
conservative estimate since we used the RK wavefunction which has higher
energy than the true ground state of Hamiltonian (9). This gives us hope that
our result obtained using perturbation theory is qualitatively correct.

Within our approximation, for $\Omega<\Omega_{C}$, the ground state is an ice
ferromagnet, an ordered state satisfying the ice rule. For $\Omega>\Omega_{C}$
but also close to $\Omega_{C}$, the ground state is in the QSL phase, i.e.,
the deconfined phase of a $U(1)$ gauge theory. From the point of view of the
QSL, the ordered ice ferromagnet state is obtained when monopole excitations
of the spin liquid proliferate, and the monopole-antimonopole string operator,
to be defined in Sec. IV.2, Eq. (53), acquires an expectation value. As a
consequence of this, the fractional “electric charges”, or spinons, get
confined [46, 47]. The monopole creation operator (see Sec. IV.2 and Ref.
[48]) is diagonal in the $\hat{S}^{z}$ basis, and acts in the sector that
obeys the ice rule. It is thus plausible that the confined phase is indeed the
ice ferromagnet. While our calculation provides microscopic intuition for this
transition, we emphasize that, to prove the existence of, locate and
characterize this transition accurately, one needs to do a more careful
quantum Monte Carlo calculation.

###  III.2 Large $\Omega$—Higgs transition

From the Hamiltonian in Eq. (6), it is clear that, in the limit $\Omega\gg V$,
the ground state is a transverse-field-polarized (TFP) state, i.e., a product
state of $(\ket{g}-\ket{r})_{i}$ at each site $i$. Thus, as $\Omega$ is
increased away from $\Omega_{C}$, the system should eventually go through a
phase transition from the putative QSL phase into the TFP phase. From the
point of view of the QSL, this is a Higgs transition because the operator
$\hat{S}^{x}$ that acquires expectation value in the TFP phase creates a pair
of “electric”-charge excitations in the spin liquid. The perturbation theory
in $\Omega/V$ that we performed in Sec. III.1 relies on the ability to go to a
basis where the Hilbert space decouples into ice rule obeying and ice rule
disobeying sectors separated by an energy gap of $V$. But the ground state in
the $\Omega\gg V$ limit (TFP) straddles both of these sectors. So we do not
expect perturbation theory in $\Omega/V$ to capture the phase transition into
the TFP phase that contains the $\Omega\to\infty$ ground state. Hence, we will
present an indirect reasoning below. In the $\Omega\ll V$ limit,
$\hat{H}_{\text{LR}}$ was important, since it was the dominant term splitting
the degeneracy in the ice manifold. On the other hand, in the vicinity of the
putative Higgs transition, $\hat{H}_{\text{LR}}$ may not be as important since
the largest term in $\hat{H}_{\text{LR}}$ has magnitude $V/27$, and as
justified above using Table 1, the effect of $\hat{H}_{\text{LR}}$ is indeed
perturbative. Therefore, we drop $\hat{H}_{\text{LR}}$ as a zeroth-order
approximation for calculating the Higgs transition point. The resulting
Hamiltonian is the transverse field Ising model on the pyrochlore lattice.
Refs. [62] and [74] studied this model and found the transition point
$\Omega_{H}$ to be at $\Omega_{H}=0.55(5)V$ and $0.6V$ respectively. This
leads us to expect that, in the window $0.44<\Omega<0.55$, the ground state
may be a QSL, leading us to sketch the phase diagram shown in Fig. 9.

Figure 9: Approximate ground state phase diagram of the Hamiltonian in Eq.
(6). The ground state for $\Omega=0$ was calculated exactly to be an ice
ferromagnet (ice FM) in Sec. III.1.2. We assume that, as $\Omega$ is
increased, no phase transition occurs to a different ordered state. The
transition point from the ice ferromagnet (confined phase) to the QSL
(deconfined phase) at $\Omega_{C}\approx 0.44V$ is obtained by comparing
energies of ansatz wavefunctions in the effective Hamiltonian obtained using
perturbation theory in $\hat{H}_{\Omega}$ and $\hat{H}_{LR}$. For the Higgs
transition to the TFP phase, we make an approximation by dropping
$\hat{H}_{LR}$, in which case $\Omega_{H}$ was calculated in Ref. [62] to be
0.55(5).

Within our approximation, $\Omega_{C}<\Omega_{H}$ and there is a window where
the QSL is the ground state. However, the introduction of
$\hat{H}_{\text{LR}}$ may result in a lowering of the energy of the TFP state
relative to the QSL. Calculating this effect and verifying that this does not
bring down $\Omega_{H}$ far enough to destroy the QSL phase requires a more
careful calculation which is beyond the scope of this work.

In the remainder of this section, we provide some intuition for the Higgs
transition by performing a gauge mean field theory (gMFT) calculation
introduced in Ref. [80].

####  III.2.1 Gauge mean field theory—Higgs transition

The main idea of this approach is to first recast the microscopic Hamiltonian
as an exact gauge theory by introducing ancillary degrees of freedom followed
by a mean-field decoupling of the interactions. This theory involves bosonic
charges hopping in the presence of a fluctuating gauge field whose mean-field
value is chosen self-consistently. If this mean-field gauge-field
configuration is such that the hopping amplitudes of the bosonic charges is 0,
then the theory is in a confined phase. If not, the theory is in the
deconfined phase as long as the bosons do not condense. If the bosonic charges
condense, then the theory is in a Higgs phase, which is adiabatically
connected to the TFP state.

Concretely, the construction is as follows. For $\mathbf{r}\in A$, where $A$
is a sublattice of the diamond lattice,

$\hat{S}^{+}_{\mathbf{r}\rightarrow\mathbf{r}+\mathbf{e}_{\mu}}=\hat{\Phi}^{\dagger}_{\mathbf{r}}\hat{\mathsf{s}}^{+}_{\mathbf{r}\rightarrow\mathbf{r}+\mathbf{e}_{\mu}}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}},$
(41)

where
$\hat{S}^{+}_{\mathbf{r}\rightarrow\mathbf{r}+\mathbf{e}_{\mu}}\equiv\hat{S}^{+}_{\mathbf{r}+\mathbf{e}_{\mu}/2}=\hat{S}^{+}_{\mathbf{r},\mu}$
(and similarly
$\hat{\mathsf{s}}^{+}_{\mathbf{r}\rightarrow\mathbf{r}+\mathbf{e}_{\mu}}\equiv\hat{\mathsf{s}}^{+}_{\mathbf{r}+\mathbf{e}_{\mu}/2}=\hat{\mathsf{s}}^{+}_{\mathbf{r},\mu}$)
lives on a bond of the diamond lattice connecting sites $\mathbf{r}$ and
$\mathbf{r}+\mathbf{e}_{\mu}$ (recall that centers of the bonds of the diamond
lattice are sites of the pyrochlore lattice). $\hat{\mathsf{s}}^{z}$ is also a
spin-1/2 operator and has eigenvalues $\pm 1/2$. Here,
$\hat{\Phi}^{\dagger}_{\mathbf{r}}$ serves as a raising operator for
$\hat{Q}_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}_{\mathbf{r}}}\equiv\eta_{\mathbf{r}}(\hat{n}_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}_{\mathbf{r}}}-2)$,
where $\eta_{\mathbf{r}}=1$ for $\mathbf{r}\in A$ and $\eta_{\mathbf{r}}=-1$
for $\mathbf{r}\in B$. For convenience, we drop the symbol  from now on.
$\hat{Q}_{\mathbf{r}}$ and $\hat{\Phi}^{\dagger}_{\mathbf{r}}$ satisfy the
commutation relation:
$\commutator{\hat{Q}_{\mathbf{r}}}{\hat{\Phi}^{\dagger}_{\mathbf{r}}}=\hat{\Phi}^{\dagger}_{\mathbf{r}}$.
Note that $\hat{\Phi}_{\mathbf{r}}$ is not a canonical boson but a rotor
satisfying

$\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\Phi}_{\mathbf{r}}=1.$ (42)

To recover the original spin Hilbert space, one imposes the constraint that
the total gauge charge at $\mathbf{r}$ is

$\hat{Q}_{\mathbf{r}}=\eta_{\mathbf{r}}\sum_{\mu}\hat{\mathsf{s}}^{z}_{\mathbf{r}+\eta\mathbf{e}_{\mu}/2}.$
(43)

Rewriting the Hamiltonian (6) in terms of the fictitious variables,
$\hat{Q}_{\mathbf{r}},\hat{\Phi}_{\mathbf{r}}$ and
$\hat{\mathsf{s}}_{\mathbf{r},\mu}$ we get

$\displaystyle\hat{H}=$ $\displaystyle\frac{V}{2}\sum_{\mathbf{r}\in
A,B}\hat{Q}^{2}_{\mathbf{r}}-\frac{\Omega}{2}\sum_{(\mathbf{r}\in
A),\mu}\left(\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\mathsf{s}}^{+}_{\mathbf{r}\rightarrow\mathbf{r}+\mathbf{e}_{\mu}}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}}+\text{
H.c.}\right)$ (44)
$\displaystyle+\frac{1}{2}\sum_{\mathbf{r},\mathbf{r}^{\prime}\in
A}\sum_{\mu,\nu}V_{\mu\nu}(\mathbf{r}-\mathbf{r}^{\prime})\hat{\mathsf{s}}^{z}_{\mathbf{r}\mu}\hat{\mathsf{s}}^{z}_{\mathbf{r}^{\prime},\nu},$

where
$V_{\mu\nu}(\mathbf{r}-\mathbf{r}^{\prime})=V\left(\frac{a}{\mathbf{r}-\mathbf{r}^{\prime}+\mathbf{e}_{\mu}/2-\mathbf{e}_{\nu}/2}\right)^{6}$
whenever $(\mathbf{r},\mu)$ and $(\mathbf{r}^{\prime},\nu)$ are distinct and
are not nearest-neighbors. $V_{\mu\nu}(\mathbf{r}-\mathbf{r}^{\prime})$ is 0
otherwise.

Following Ref. [80], we perform the zeroth-order mean-field decoupling:
$\hat{\Phi}^{\dagger}\hat{\Phi}\hat{\mathsf{s}}\rightarrow\hat{\Phi}^{\dagger}\hat{\Phi}\expectationvalue{\hat{\mathsf{s}}}+\expectationvalue{\hat{\Phi}^{\dagger}\hat{\Phi}}\hat{\mathsf{s}}-\expectationvalue{\hat{\Phi}^{\dagger}\hat{\Phi}}\expectationvalue{\hat{\mathsf{s}}}$
and
$\hat{\mathsf{s}}\hat{\mathsf{s}}\rightarrow\hat{\mathsf{s}}\expectationvalue{\hat{\mathsf{s}}}+\expectationvalue{\hat{\mathsf{s}}}\hat{\mathsf{s}}-\expectationvalue{\hat{\mathsf{s}}}\expectationvalue{\hat{\mathsf{s}}}$
(where $\hat{\mathsf{s}}$ could either be $\hat{\mathsf{s}}^{+}$,
$\hat{\mathsf{s}}^{-}$, or $\hat{\mathsf{s}}^{z}$). Upon doing so, the
Hamiltonian decouples into a Hamiltonian of bosons hopping on the diamond
lattice and a Hamiltonian of spins in an external field, which itself is set
self-consistently by the Green’s function of the bosons. Before solving the
resulting theory, one needs to enforce the constraints (42) and (43) using
Lagrange multipliers $\lambda_{\mathbf{r}}$ and $v_{\mathbf{r}}$,
respectively. Within the mean-field theory, it is assumed that these Lagrange
multipliers take a spatially homogeneous value at the saddle point. We then
find the minimum value of $\Omega^{\text{MF}}_{H}$ such that, for any
$\Omega>\Omega^{\text{MF}}_{H}$, it is possible to self-consistently choose
$\lambda$ only by macroscopically occupying a boson mode. This
$\Omega^{\text{MF}}_{H}$ marks the location of the Bose-Einstein-condensation
transition (or Higgs transition within the mean-field theory). We find
$\Omega^{\text{MF}}_{H}\approx 0.7V$. In Appendix A, we present more details
of this calculation. An artifact of this technique is that, although we
include long-range interactions in our calculation, they do not play any role
at the saddle point near the Higgs transition. Therefore, the final steps and
result of our calculation are identical to the ones carried out in [68].

In Appendix A, we also point out a major limitation of this technique in the
small-$\Omega$ limit that may not have been appreciated in previous literature
on gauge mean field theory.

###  III.3 Comments on dynamical state preparation

So far, we have focused on the nature of the ground state of Hamiltonian (6)
as a function of $\Omega/V$. However, what is often experimentally relevant is
the nature of the state prepared by a ramping of parameters in a finite amount
of time. This issue was studied in Ref. [81] for a $\mathbb{Z}_{2}$ gauge
theory in the context of the experiment in Ref. [26]. Since the main ideas of
Ref. [81] are general enough, here we will present an adaptation of the
conclusions of Ref. [81] to our setting.

The excitations of a $U(1)$ QSL are gapless “photons”, magnetic monopoles, and
“electric charges” (spinons). The transition of a QSL to an ice ferromagnet is
driven by the condensation of monopoles, while the transition to the TFP phase
is driven by the condensation of spinons. The gapless “photons” are not
directly involved in these transitions. Also, a state with “photon” modes
excited on top of a QSL state is still in the deconfined phase of the $U(1)$
gauge theory. This allows us to ignore “photons” in this section. Since the
confined phase, ice ferromagnet has an extensive number of monopoles, we use
the difference per unit cell between the energies of the QSL and ice
ferromagnet states as a proxy for the monopole energy scale. At $\Omega=0$,
this difference is
$\left(\expectationvalue{\hat{H}_{\text{LR}}}_{\text{QSL}}-\expectationvalue{\hat{H}_{\text{LR}}}_{\text{IFM}}\right)/N_{\text{u.c.}}\approx
0.03V$ (see Table 1), which is much smaller than the spinon energy scale (see
Fig. 10 for a sketch).

Figure 10: A qualitative sketch of the energy scales (per unit cell) in our
problem. For $\Omega>\Omega_{C}$, the ground state is a $U(1)$ QSL. Ice
ferromagnet is the ordered state obtained when monopoles proliferate, i.e.,
the ice ferromagnet has an extensive number of monopoles. We therefore use the
energy difference per unit cell between the QSL and the ice ferromagnet at
$\Omega=0$, obtained in Table 1, as a proxy for the monopole energy scale.
This scale $\sim 0.03V$ is much smaller than the spinon energy scale
(“electric charge”), which is $\sim$ V.

Suppose one starts with an initial state (for a small $\epsilon\sim\Omega/V$)

$\ket{\Psi_{(t=0)}}=\otimes_{i}\left(\ket{g}_{i}+\epsilon\ket{r}_{i}\right),$
(45)

which is the ground state in the limit of large negative $\delta/V$ and small
$\Omega/V$. As shown in Sec. II, the classical ground state lies in the ice
manifold when $\delta\in\left(2.46V,4.46V\right)$. Now suppose that $\delta$
is ramped up from its initial large negative value to a value in this range
such that the ramp is adiabatic with respect to the spinon gap $V$, but is
sudden with respect to the monopole scale $\sim 0.03V$, while keeping
$\Omega/V\ll 1$. Using arguments in Ref. [81], this protocol will not prepare
the ground state, which, from Fig. 9, is an ice ferromagnet. Instead, it will
(approximately) project out violations of the ice rule (due to adiabaticity
with respect to the spinon scale) from the initial state $\ket{\Psi_{(t=0)}}$.
The resulting final state is

$\ket{\Psi_{\text{final}}}\approx\hat{\mathcal{P}}\left\\{\otimes_{i}\left(\ket{g}_{i}+\epsilon\ket{r}_{i}\right)\right\\}=\ket{\Psi_{\textrm{RK}}},$
(46)

where $\hat{\mathcal{P}}$ is the projector onto the ice manifold. The
projected wavefunction is an equal-weight superposition of all coverings,
which is simply the RK wavefunction and which lies in the QSL phase [48].
There is one catch to the above argument—the spinon gap closes during the
above ramp. So it is impossible to be sudden with respect to the monopole
scale and yet be strictly adiabatic with respect to the spinon gap throughout
the ramp. For a short duration (while the ramp is going through the spinon gap
closing), adiabaticity with respect to the spinon gap will be violated. By the
Kibble-Zurek mechanism, the resulting state is composed of finite-size puddles
of QSL-like regions with a nonzero density of spinons interspersed [81, 82,
83, 84]. Thus, in summary, there are two different ways in which one can
prepare a $U(1)$ QSL-like state in experiment and study a confinement-
deconfinement transition111We note that the confinement-deconfinement
transition of $U(1)$ gauge theory in 3+1D is strictly speaking, a ground state
transition [85, 71]. Therefore, in this paper, when we use the phrase
confinement-deconfinement transition, we mean signatures of this transition in
a finite-size state prepared in finite time..

  1. 1.

$\Omega/V\ll 1$: Perform a ramp of $\delta$ starting from a large negative
value and ending in the range $(2.46V,4.46V)$ for a fixed $\Omega/V\ll 1$ such
that the ramp is adiabatic with respect to $V$ (spinon gap) but sudden with
respect to the monopole scale ($\sim 0.03V$). Even though the ground state is
not a QSL for these parameters, this procedure would create puddles of QSL-
like regions by the argument in Ref. [81]. To see a deconfinement-confinement
transition, the ramp of $\delta$ should be slowed down and, once it is
adiabatic with respect to the monopole gap, an ordered, i.e. confined state
will be prepared.

  2. 2.

Adiabatic: Perform a ramp of $\delta$ starting from a large negative value and
ending in the range $(2.46V,4.46V)$ and a ramp in $\Omega$ starting from
$\Omega/V\ll 1$ and ending in a final value $\Omega_{f}$, such that both ramps
are adiabatic with respect to the monopole scale always. The two ramps can be
performed simultaneously, or such that the ramp in $\delta$ precedes the ramp
in $\Omega$. This would approximately create the ground state of Hamiltonian
(6). As the final value $\Omega_{f}$ goes through $\Omega_{C}$ ($\Omega_{H}$),
the nature of the final state prepared this way goes through a confinement-
deconfinement (Higgs) transition.

Once a state is prepared by either of the above schemes, one needs to devise
measurements that can tell whether the state is in the confined phase or in
the deconfined phase. We address this in the following section.

##  IV Diagnosis of the quantum spin liquid

Access to wavefunction snapshots in the $\hat{S}^{z}$ basis, combined with
access to unitary evolution, allows one to use the Rydberg-atom platform to
measure non-local observables, a feature generally unavailable in traditional
condensed matter systems. In this section, we describe some measurable
correlators which can be used to observe the signatures of a quantum spin
liquid state. In this section, we assume that the detuning is chosen such that
$\rho=2$.

###  IV.1 Plaquette-plaquette correlators

The plaquette operators are off-diagonal in the $\hat{S}^{z}$ basis and map
one dimer covering to another. They are important in distinguishing a coherent
quantum superposition from a classical admixture of states. We define two
plaquette operators $\hat{X}_{P}$ and $\hat{Y}_{P}$ for a hexagonal plaquette
$P$ of the pyrochlore lattice as

$\displaystyle\hat{X}_{P}$
$\displaystyle=\left(\hat{S}_{1}^{+}\hat{S}_{2}^{-}\hat{S}_{3}^{+}\hat{S}_{4}^{-}\hat{S}_{5}^{+}\hat{S}_{6}^{-}+\text{H.c.}\right),$
(47) $\displaystyle\hat{Y}_{P}$
$\displaystyle=-i\hat{S}_{1}^{+}\hat{S}_{2}^{-}\hat{S}_{3}^{+}\hat{S}_{4}^{-}\hat{S}_{5}^{+}\hat{S}_{6}^{-}+\text{H.c.},$

where $1,2,\ldots,6$ denote the sites around a plaquette $P$. Either of the
two correlators, $\langle\hat{X}_{P}\hat{X}_{P^{\prime}}\rangle$ and
$\langle\hat{Y}_{P}\hat{Y}_{P^{\prime}}\rangle$, of the plaquette operators on
two plaquettes $P$ and $P^{\prime}$ can distinguish a QSL phase from other
phases including a classical spin ice (see Table 2).

Correlator | Confined (Ice FM) | Deconfined (QSL) | Higgs (TFP) | Classical Spin Ice
---|---|---|---|---
$\langle\hat{X}_{P}\hat{X}_{P^{\prime}}\rangle-4$ | Exp. or faster decay | $1/R^{8}$ | Nonzero const. | Exp. or faster decay
$\langle\hat{Y}_{P}\hat{Y}_{P^{\prime}}\rangle$ | Exp. or faster decay | $1/R^{4}$ | Exp. or faster decay | Exp. or faster decay
$\expectationvalue{\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})}$ | Nonzero const. | Exp. decay | Exp. decay | Exp. or faster decay
$\chi_{\mathcal{C}}^{E}$ | Nonzero const.222Distinguishing this non-zero constant from zero for $\chi_{\mathcal{C}}^{E}$ in the confined phase (Ice FM) and for $\chi_{\mathcal{C}}^{M}$ in the Higgs phase (TFP) may be practically challenging. | Exp. or faster decay | Nonzero const. | Exp. or faster decay
$\chi_{\mathcal{C}}^{M}$ | Nonzero const. | Exp. or faster decay | Nonzero const.2 | Exp. or faster decay
$\langle\hat{S}^{z}_{\mathbf{r}i}\hat{S}^{z}_{\mathbf{r^{\prime}}j}\rangle$ | Nonzero const. | $1/R^{4}$ | Exp. decay | $1/R^{4}$
Table 2: Behavior of various correlators. $\hat{X}_{P}$ and $\hat{Y}_{P}$ are
plaquette operators defined in Eq. (47).
$\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})$
is a monopole string operator defined in Eq. (53). $\chi_{\mathcal{C}}^{E}$
and $\chi_{\mathcal{C}}^{M}$ are BFFM order parameters defined in Eq. (60) and
Eq. (64), respectively. In this table, we have omitted the form factors
multiplying $1/R^{4}$ and $1/R^{8}$ that are provided in Eqs. (50) and (65).

We compare the two correlators and provide protocols to measure them. We
assume throughout that the two plaquettes $P$ and $P^{\prime}$ do not have any
sites in common.

Using the mapping between the spins and the effective $U(1)$ gauge theory from
Eq. (13), we see that the operators $\hat{X}_{P}$ and $\hat{Y}_{P}$ are equal
to (twice) the cosine and the sine of the magnetic field, respectively:

$\displaystyle\hat{X}_{P}$
$\displaystyle=2\cos(\hat{\theta}_{1}-\hat{\theta}_{2}+\hat{\theta}_{3}-\hat{\theta}_{4}+\hat{\theta}_{5}-\hat{\theta}_{6})=2\cos(\hat{b}_{\bm{\mathsf{r}},\mu}),$
(48) $\displaystyle\hat{Y}_{P}$
$\displaystyle=2\sin(\hat{\theta}_{1}-\hat{\theta}_{2}+\hat{\theta}_{3}-\hat{\theta}_{4}+\hat{\theta}_{5}-\hat{\theta}_{6})=2\sin(\hat{b}_{\bm{\mathsf{r}},\mu}),$

where $\bm{\mathsf{r}}$ belongs to the dual diamond lattice [see Fig. 1(c)],
and $\mu\in\\{0,1,2,3\\}$ labels the direction of magnetic field.
$\hat{b}_{\bm{\mathsf{r}},\mu}$ is along $\bm{\mathsf{u}}_{\mu}$, which are
vectors joining an A site of the dual diamond lattice to its neighboring B
sites. These vectors are perpendicular to the plaquettes of the pyrochlore
lattice, see Fig. 11.

Figure 11: Notation for the plaquette correlators. $P$ and $P^{\prime}$ are
two hexagonal plaquettes of the pyrochlore lattice.
$\bm{\mathsf{r}},\bm{\mathsf{r^{\prime}}}$,
$\bm{\mathsf{r}}+\bm{\mathsf{u}}_{\mu}$, and
$\bm{\mathsf{r^{\prime}}}+\bm{\mathsf{u}}_{\nu}$ are the sites of the dual
diamond lattice. $\bm{\mathsf{u}}_{\mu}$ and $\bm{\mathsf{u}}_{\nu}$ are
vectors perpendicular to $P$ and $P^{\prime}$.

The effective theory in the deconfined phase (QSL) is Maxwell
electromagnetism.

In the $3+1$D continuum Maxwell electromagnetism, the correlators of the
Cartesian components of the magnetic field $\hat{b}_{\bm{\mathsf{r}},i}$ for
$i\in\\{x,y,z\\}$ can be expressed as Gaussian integrals which evaluate to
[48]

$\langle\hat{b}_{\bm{0},i}\hat{b}_{\bm{\mathsf{R}},j}\rangle_{0}\propto\frac{1}{\mathsf{R}^{4}}\left(2\frac{\mathsf{R}_{i}\mathsf{R}_{j}}{\mathsf{R}^{2}}-\delta_{ij}\right)\equiv\mathcal{C}_{ij}^{B}(\bm{\mathsf{R}}),$
(49)

where $\langle\cdot\rangle_{0}$ is the expectation value with respect to the
Maxwell action.

The correlators of the magnetic field operators
$\hat{b}_{\bm{\mathsf{r}},\mu}$ for $\mu\in\\{0,1,2,3\\}$ on the pyrochlore
plaquettes in the deconfined phase are obtained from Eq. (49) by taking
components of the Cartesian magnetic field along the vectors
$\bm{\mathsf{u}}_{\mu}$. The result is

$\displaystyle\langle\hat{X}_{P}\hat{X}_{P^{\prime}}\rangle$
$\displaystyle-4\propto\frac{1}{\mathsf{R}^{8}}\left[\sum_{k,l}(\bm{\mathsf{u}}_{\mu})_{k}(\bm{\mathsf{u}}_{\nu})_{l}\left(2\frac{\mathsf{R}_{l}\mathsf{R}_{k}}{\mathsf{R}^{2}}-\delta_{k,l}\right)\right]^{2},$
(50) $\displaystyle\langle\hat{Y}_{P}\hat{Y}_{P^{\prime}}\rangle$
$\displaystyle\propto\frac{1}{\mathsf{R}^{4}}\left[\sum_{k,l}(\bm{\mathsf{u}}_{\mu})_{k}(\bm{\mathsf{u}}_{\nu})_{l}\left(2\frac{\mathsf{R}_{l}\mathsf{R}_{k}}{\mathsf{R}^{2}}-\delta_{k,l}\right)\right],$

where the summation is over $k,l\in\\{x,y,z\\}$,
$\bm{\mathsf{R}}=\bm{\mathsf{r}}-\bm{\mathsf{r}}^{\prime}$, and $\mathsf{R}$
is assumed to be large compared to the monopole correlation length. The
factors inside the square brackets are geometric factors, which depend on the
direction of the vectors $\bm{\mathsf{u}}_{\mu}$, $\bm{\mathsf{u}}_{\nu}$, and
$\bm{\mathsf{R}}$, but are independent of the distance $\mathsf{R}$ between
the two plaquettes. Ref. [48] also separately studied the correlators
precisely at the RK point (which sits at the phase boundary between deconfined
and confined phases) where the effective field theory differs from the regular
Maxwell theory. In the RK wavefunction, while the behavior of the plaquette
correlators differs from Eq. (50), it is still a power law with a slower decay
[48]. We note that, if the experimentally prepared state is close to an RK
wavefunction (see discussion in Sec. III.3), then this distinction will be
important.

Outside of the deconfined phase (QSL) of the compact $U(1)$ gauge theory,
Maxwell electromagnetism is no longer the effective theory, and the behavior
of the correlators is different. The ice ferromagnet state is an ordered state
with the spins pointing in the $z$-direction, thus the two plaquette
correlators are expected to decay exponentially with $\mathsf{R}$. The TFP
phase has $\langle\hat{S}^{+}\rangle\neq 0$, and hence the plaquette X
correlator approaches a nonzero constant at large $\mathsf{R}$, while the
plaquette Y correlator is 0 (or decays exponentially or faster in
$\mathsf{R}$). Lastly, in a classical spin liquid, which consists of an
incoherent mixture of exponentially many dimer coverings, the plaquette
correlators decay at least exponentially (see Table 2).

Since the plaquette correlators involve off-diagonal operators, they cannot be
read out directly from the snapshots of a Rydberg-atom array. However, we show
that they can be measured by evolving the system under a modified Hamiltonian
for a specific time duration followed by measurement of a diagonal operator
[37, 26]. We describe the protocols to measure both plaquette X and plaquette
Y correlators in the sections below.

####  IV.1.1 Plaquette X correlator

For a state $\ket{\Psi_{0}}$ completely within the ice manifold, the
expectation value of $\hat{X}_{P}\hat{X}_{P^{\prime}}$ is the same as that of
a product of $\hat{S}^{x}$ operators on the 12 sites of $P$ and $P^{\prime}$,
that is,

$\langle\Psi_{0}|\hat{X}_{P}\hat{X}_{P^{\prime}}|\Psi_{0}\rangle=\langle\Psi_{0}|\prod_{i=1}^{12}2\hat{S}^{x}_{i}|\Psi_{0}\rangle,$
(51)

where the sites $i=1,2,\ldots,6$ are on the plaquette $P$ and sites
$i=7,8,\ldots,12$ are on the plaquette $P^{\prime}$. This can be seen by
writing $2\hat{S}^{x}_{i}=\hat{S}^{+}_{i}+\hat{S}^{-}_{i}$ and noticing that
the only terms that preserve the ice rule are ring exchanges over $P$ and
$P^{\prime}$. When the remaining terms act on a state in the ice manifold,
they either take the state outside of the ice manifold or annihilate it. Thus
the expectation value of these remaining operators in $\ket{\Psi_{0}}$ is
zero. For example,
$\hat{S}_{1}^{+}\hat{S}_{2}^{+}\hat{S}_{3}^{+}\hat{S}_{4}^{-}\hat{S}_{5}^{-}\hat{S}_{6}^{-}\hat{S}_{7}^{+}\hat{S}_{8}^{-}\hat{S}_{9}^{+}\hat{S}_{10}^{-}\hat{S}_{11}^{+}\hat{S}_{12}^{-}$
acting on a state in the ice manifold would either annihilate this state or
give a state that violates the ice rule on four of the tetrahedra surrounding
$P$. The protocol to measure the plaquette X correlator is as follows.

We abruptly change the phase and the amplitude of the Rabi frequency, so that
the new Hamiltonian is
$\hat{H}_{\text{Y}}\approx\Omega_{Y}\sum_{i}\hat{S}^{y}_{i}$ with
$\Omega_{Y}\gg V$. (Achieving $\Omega_{Y}\gg V$ may require working with atom
spacings that are sufficiently large and/or with Rydberg principal quantum
numbers that are sufficiently low, but not low enough to make Rydberg lifetime
a problem.) It is assumed that this change of the Hamiltonian is done
sufficiently rapidly so that the sudden approximation is valid and the state
of the system does not change. We evolve the system under $\hat{H}_{Y}$ for a
time $t_{Y}=\pi/(2\Omega_{Y})$, which amounts to a $\pi/2$ pulse about the
$y$-axis, and then measure all atoms in the $\\{\ket{g},\ket{r}\\}$ basis.
Thanks to Eq. (51), this allows us to compute the plaquette X correlator.

Eq. (51) assumes that the state $\ket{\Psi_{0}}$ is in the ice manifold.
However, the ground state $\ket{\Psi_{g}}$ of the system is not completely in
the ice manifold. The error introduced because of assuming $\ket{\Psi_{g}}$ to
be in the ice manifold is of sixth order in $\Omega/V$, as we show in Appendix
B. The intuitive reason for this is that six factors of $\hat{H}_{\Omega}$ are
required to give a state that can have nonzero overlap with
$\prod_{i=1}^{12}2\hat{S}_{i}^{x}\ket{\Psi_{0}}$.

We note that, for $\Omega/V\gg 1$ (i.e. in the TFP phase), the experimentally
measured quantity
$\langle\Psi_{g}|\prod_{i=1}^{12}2\hat{S}_{i}^{x}|\Psi_{g}\rangle$ is not the
same as the plaquette X correlator
$\langle\Psi_{0}|\hat{X}_{P}\hat{X}_{P^{\prime}}|\Psi_{0}\rangle$, but their
behaviors are nevertheless the same, i.e., both approach a nonzero constant as
the distance between the plaquettes increases.

The plaquette X correlator behaves as
$\langle\hat{X}_{P}\hat{X}_{P^{\prime}}\rangle-4\propto 1/R^{8}$ (the angular
dependence is not shown here) in the QSL phase. Although it is a power law,
the decay is very rapid, and it might be practically difficult to distinguish
it from an exponential decay. This issue is less prominent in the plaquette Y
correlator, and we now provide a protocol to measure it.

####  IV.1.2 Plaquette Y correlator

This protocol relies on the fact that, for a state $\ket{\Psi_{0}}$ within the
ice manifold, we have an identity similar to Eq. (51):

$\langle\Psi_{0}|\hat{Y}_{P}\hat{Y}_{P^{\prime}}|\Psi_{0}\rangle=\langle\Psi_{0}|\prod_{i=1}^{6}(2\hat{S}^{x}_{2i-1})(2\hat{S}^{y}_{2i})|\Psi_{0}\rangle.$
(52)

Recall that sites $1,2,\ldots,6$ are on $P$, while sites $7,8,\ldots,12$ are
on $P^{\prime}$. This protocol is the same as the protocol for measuring the
plaquette X correlator, except that now the $\pi/2$ pulses on sites $2i$ for
$i=1,2,\ldots,6$ are about the $x$-axis on the Bloch sphere while the $\pi/2$
pulses on sites $2i-1$ for $i=1,2,\ldots,6$ are around the $y$-axis. After
applying these $\pi/2$ pulses, $\prod_{i=1}^{12}2\hat{S}_{i}^{z}$ is measured
by taking a snapshot of the array. Note that, for any set of non-overlapping
plaquettes, we can simultaneously measure $Y$ correlators between all pairs of
plaquettes in this set.

This protocol requires control over individual sites and has an error of order
$(\Omega/V)^{6}$ (the estimation of the error is similar to the one for the
plaquette X correlator done in Appendix B). The advantage of the plaquette Y
correlator over the plaquette X correlator is that the former decays as
$1/R^{4}$, which is slower than the decay of the plaquette X correlator and
can be easier to observe experimentally.

###  IV.2 Monopole-monopole correlator

In the deconfined phase, monopoles are gapped. Therefore, the expectation
value of an (equal-time) operator that creates a string with a monopole and
antimonopole at its endpoints should decay exponentially with the length of
the string. On the other hand, in the confined phase, monopoles are condensed,
and hence the expectation value should approach a nonzero constant as the
length of the string increases. In the continuum, the following operator
inserts a string that creates a monopole at $\bm{\mathsf{r}}_{1}$ and an
antimonopole at $\bm{\mathsf{r}}_{2}$ [48]:

$\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})\sim
e^{i\int\differential[3]{\bm{\mathsf{r^{\prime}}}}\mathbf{\mathcal{A}}(\bm{\mathsf{r^{\prime}}})\cdot\hat{\bm{e}}(\bm{\mathsf{r^{\prime}}})}.$
(53)

Here $\mathbf{\mathcal{A}}(\bm{\mathsf{r^{\prime}}})$ is a classical (non
single-valued) vector potential such that the flux $\phi_{\Sigma}$ of
$\mathbf{\mathcal{B}}=\gradient\times\mathbf{\mathcal{A}}$ through a closed
surface $\Sigma$ is

$\phi_{\Sigma}\equiv\oint_{\Sigma}\mathbf{\mathcal{B}}\cdot d\mathbf{S}=2\pi
qQ_{\Sigma},$ (54)

where $Q_{\Sigma}=1$ when $\Sigma$ encloses $\bm{\mathsf{r}}_{1}$ and not
$\bm{\mathsf{r}}_{2}$, $Q_{\Sigma}=-1$ when $\Sigma$ encloses
$\bm{\mathsf{r}}_{2}$ and not $\bm{\mathsf{r}}_{1}$, and $Q_{\Sigma}=0$
otherwise. $q$ is an integer and denotes the “charge” of the monopole string.
For simplicity, we will set $q=1$ in this section. We clarify that
$\mathbf{\mathcal{B}}$ and $\phi_{\Sigma}$ are classical numbers and are
different from $\hat{\mathbf{b}}$ and $\hat{\Phi}_{\Sigma}$ which are
operators. $\hat{\mathbf{b}}\equiv\gradient\times\hat{\mathbf{a}}$, for gauge-
field (operator) $\hat{\mathbf{a}}$, and $\hat{\Phi}_{\Sigma}$ is defined as

$\hat{\Phi}_{\Sigma}\equiv\oint_{\Sigma}\hat{\mathbf{b}}\cdot
d\mathbf{S}=2\pi\hat{m},$ (55)

where $\hat{m}$ takes integer eigenvalues. The form of the monopole string
operator is chosen so that it increases the flux through $\Sigma$ by $2\pi
Q_{\Sigma}$, i.e.,

$\commutator{\hat{\Phi}_{\Sigma}}{\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})}=2\pi
Q_{\Sigma}\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2}).$
(56)

Figure 12: (a) The “polyhedron” formed by four puckered hexagons of the
diamond lattice is shown in orange. The centers of these “polyhedra” form the
dual diamond lattice. (b) The center of the “polyhedron” in (a) is also the
center of a truncated tetrahedron (shown in red) of the pyrochlore lattice.

We now adapt this operator to the Rydberg setting. Consider the diamond
lattice formed by the centers of tetrahedra of the pyrochlore lattice, Fig.
1(b). Unlike the continuum, it is now important to specify that the endpoints
of the monopole string $\bm{\mathsf{r}}_{1}$ and $\bm{\mathsf{r}}_{2}$ belong
to the dual diamond lattice [see Fig. 1(c)], whose sites are centers of
“polyhedra” made of four puckered-hexagonal “plaquettes” of the diamond
lattice333In terms of the original pyrochlore lattice, the vertices of the
dual diamond lattice are centers of the truncated tetrahedra [see Fig. 12(b)]
which fill the voids between the tetrahedra., see Fig. 12(a). Let
$\bm{x}\equiv\mathbf{r}+\mathbf{e}_{\mu}/2$ be a site on the pyrochlore
lattice, where $\mathbf{r}$ is an A-site of the diamond lattice.
$\mathcal{A}_{\bm{x}}\equiv\mathcal{A}_{\mathbf{r},\mathbf{r}+\mathbf{e}_{\mu}}$
is the discrete version of $\mathcal{A}$ integrated (Fig. 1(b) shows the
vectors $\mathbf{e}_{\mu}$). along the line pointing from the center of an A
tetrahedron (centred at $\mathbf{r}$) to the B tetrahedron (centred at
$\mathbf{r}+\mathbf{e}_{\mu}$) such that the two tetrahedra touch at $\bm{x}$.

$\mathcal{A}_{\bm{x}}$ is required to satisfy the discrete version of Eq.
(54), and hence depends on $\bm{\mathsf{r}}_{1}$, $\bm{\mathsf{r}}_{2}$, the
“magnetic field” configuration $\mathbf{\mathcal{B}}$ and the gauge choice for
$\mathcal{A}_{\bm{x}}$. For the pyrochlore lattice, we have

$\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})=e^{i\sum_{\bm{x}\in\text{pyrochlore}}\mathcal{A}_{\bm{x}}\left(\hat{n}_{\bm{x}}-1/2\right)}.$
(57)

This operator is purely diagonal in the $\hat{n}_{\bm{x}}$ basis (i.e., in the
$\hat{S}^{z}$-basis). So, experimentally, one can calculate this phase for
each snapshot and average over shots.

Theoretically, one expects

$\absolutevalue{\expectationvalue{\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})}}\sim\begin{cases}e^{-\absolutevalue{\bm{\mathsf{r}}_{2}-\bm{\mathsf{r}}_{1}}/\lambda},&\text{deconfined
phase,}\\\ \text{constant,}&\text{confined phase},\end{cases}$ (58)

where $\lambda$ is a correlation length that depends on the monopole gap and
the “photon” velocity.

Figure 13: An example of the monopole string operator
$\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})$
for which we provide $\mathcal{A}_{\bm{x}}$ explicitly. In our example, the
string carries $2\pi$ flux through a tube with a width of 7 puckered hexagons
of the diamond lattice. (a) A schematic of the tube running along the
$z$-direction. The diamond lattice (whose vertices are centers of tetrahedra
of the pyrochlore lattice) can be seen as ABC stacking of layers of
“honeycomb” lattices made of chair-like puckered hexagons. The tube consists
of three types of layers shown in yellow, orange, and cyan. Each layer is made
of 7 puckered hexagons. To convey a sketch, we depict such a layer by a big
hexagon with some thickness. (b) A side view of the stack showing three of its
layers, which then need to be repeated in the $z$ direction to get the entire
string. For bonds $\bm{x}$ with arrows, the value of $\mathcal{A}_{\bm{x}}$ is
written next to the bond. For bonds $\bm{x}$ without arrows,
$\mathcal{A}_{\bm{x}}=0$. (c) Top view of three of the layers of the stack. It
can be seen from all three sub-figures (a)-(c) that the flux through any
closed surface $\Sigma$ that completely encloses an integer number of layers,
such that the bottom layer is included but not the top, is $2\pi$. However, if
$\Sigma$ partially encloses a layer, then $\Phi_{\Sigma}$ is 0. This
difficulty in defining arbitrary integer multiples of $2\pi$ flux through a
volume enclosed by a finite number of plaquettes has been observed before
[48]. Therefore in our construction, $\bm{\mathsf{r}}_{1}$ and
$\bm{\mathsf{r}}_{2}$ have to be seen as being smeared across 7 points of the
dual diamond lattice below the bottom layer and above the top layer
respectively, in order to be consistent with Eq. (54).

In Fig. 13, we provide an example of one configuration of the classical
numbers $\mathcal{A}_{\bm{x}}$ that defines a monopole string operator. Below,
we comment on the freedom in choosing $\mathcal{A}_{\bm{x}}$.

####  IV.2.1 Choice of $\mathcal{A}$

The classical numbers $\mathcal{A}_{\bm{x}}$ should of course obey the
constraint that the flux of $\curl\mathcal{A}$ through a closed surface
$\Sigma$ is $2\pi Q_{\Sigma}$, as mentioned above. However, one still has a
freedom in the choice of $\mathcal{A}$ in the following two respects:

  1. 1.

Freedom in the arrangement of the field lines of $\curl\mathcal{A}$. For
example, they can be confined to a thin tube connecting $\bm{\mathsf{r}}_{1}$
and $\bm{\mathsf{r}}_{2}$, or be spread out according to Coulomb’s law, or be
something in between. Different such arrangements, due to their different
energy costs, would differ in sub-leading corrections to the exponentially
decaying behavior, but the leading behavior would be unchanged. In Fig. 13, we
provide a choice of $\mathcal{A}$, such that the monopole string is localized
to a thin tube.

  2. 2.

For a fixed choice of field lines, we still have a gauge choice for
$\mathcal{A}$. Consider a gauge transformation
$\mathcal{A}_{\mathbf{r},\mathbf{r}+\mathbf{e}_{\mu}}\to\mathcal{A}_{\mathbf{r},\mathbf{r}+\mathbf{e}_{\mu}}+\lambda_{\mathbf{r}+\mathbf{e}_{\mu}}-\lambda_{\mathbf{r}}$,
where $\lambda_{\mathbf{r}}$ is an $\mathbf{r}$-dependent real number. It
results in

$\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})\to\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\rightarrow\bm{\mathsf{r}}_{2})e^{-i\sum_{\mathbf{r}}\lambda_{\mathbf{r}}\eta_{\mathbf{r}}(\hat{n}_{\mathchoice{\includegraphics[height=2.15277pt]{tetrahedron}}{\includegraphics[height=3.6597pt]{tetrahedron}}{\includegraphics[height=3.22916pt]{tetrahedron}}{\includegraphics[height=1.07639pt]{tetrahedron}}_{\mathbf{r}}}-2)},$
(59)

where $\eta_{\mathbf{r}}=1$ for $\mathbf{r}\in A$ and $\eta_{\mathbf{r}}=-1$
for $\mathbf{r}\in B$. In the $\Omega/V\ll 1$ limit, we have
$\hat{n}_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}_{\mathbf{r}}}=2$,
so the expectation value is invariant under the gauge transformation. Away
from this limit, a gauge transformation on
$\mathcal{A}_{\mathbf{r},\mathbf{r}+\mathbf{e}_{\mu}}$ generically results in
a physical transformation on the monopole string operator. However, as long as
the external field $h=0$ [$h$ is defined in Eq. (5)], by particle-hole
symmetry, we have
$\expectationvalue{\hat{n}_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}_{\mathbf{r}}}}=2$,
and
$\expectationvalue{(\hat{n}_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}_{\mathbf{r}}}-2)^{2}}$
is bounded. Hence we do not expect the gauge transformation on
$\mathcal{A}_{\mathbf{r},\mathbf{r}+\mathbf{e}_{\mu}}$ to qualitatively change
the behavior of Eq. (58).

###  IV.3 BFFM order parameter

It is known that the confined and deconfined phases of a gauge theory without
matter fields can be distinguished by the scaling of the Wilson loops
$W_{\mathcal{L}}=\expectationvalue{e^{i\oint_{\mathcal{L}}A_{\mu}dx^{\mu}}}$,
where $A_{\mu}$ is the gauge field and $\mathcal{L}$ is a closed loop. In the
deconfined phase, the Wilson loop follows the perimeter law,
$W_{\mathcal{L}}\propto e^{-\text{Perimeter of }\mathcal{L}}$, while in the
confined phase, it follows the area law, $W_{\mathcal{L}}\propto
e^{-\text{Area of }\mathcal{L}}$. However, in the presence of matter fields
(which are generically always present), the Wilson loop follows the perimeter
law in both phases [86, 87], and it cannot be used to distinguish them. The
Bricmont-Frölich-Fredenhagen-Marcu (BFFM) order parameter is useful in such
cases and has a different behavior in the two phases [88, 89, 90, 91, 92, 93,
37, 26]. The BFFM order parameter, denoted here by $\chi^{E}_{\mathcal{C}}$,
is defined as

$\chi_{\mathcal{C}}^{E}=\frac{\absolutevalue{\expectationvalue{e^{i\sum_{\mathcal{C}}\hat{a}_{\mathbf{r}\mathbf{r}^{\prime}}}+\text{H.c.}}}}{\sqrt{\absolutevalue{\expectationvalue{e^{i\sum_{\mathcal{L}}\hat{a}_{\mathbf{r}\mathbf{r}^{\prime}}}+\text{H.c.}}}}},$
(60)

where $\mathcal{C}$ is an open curve and $\mathcal{L}$ is the closed loop
formed by combining $\mathcal{C}$ with its mirror image about a plane that
intersects $\mathcal{C}$ only at its end points. This order parameter detects
long-range order in the “electric charge”-creation string. In the Higgs phase,
“electric charges” are condensed, and hence $\chi_{\mathcal{C}}^{E}$
approaches a nonzero constant. In the deconfined phase, the numerator in Eq.
(60) (calculated on an open curve) decays to zero faster than the denominator
(calculated on a closed loop, giving the Wilson loop), as the length of
$\mathcal{C}$ is increased. Therefore, in the deconfined phase,
$\chi_{\mathcal{C}}^{E}$ goes to 0 as the length of $\mathcal{C}$ is
increased. In the confined phase, it was argued in Ref. [89] that while both
the numerator and the denominator go to zero as the length of $\mathcal{C}$ is
increased, the limit of their ratio approaches a constant. However
distinguishing this constant from zero in finite systems for finite length of
$\mathcal{C}$ may be difficult. Below we explain how to measure
$\chi_{\mathcal{C}}^{E}$.

Using the mapping from spin operators to gauge fields, Eqs. (13) and (15), we
see that

$e^{i\sum_{\mathcal{C}}\hat{a}_{\mathbf{r}\mathbf{r}^{\prime}}}\simeq\hat{S}_{1}^{+}\hat{S}_{2}^{-}\hat{S}_{3}^{+}\cdots,$
(61)

where the product of $\hat{S}^{+}$ and $\hat{S}^{-}$ operators is over the
sites on the curve $\mathcal{C}$. The denominator in $\chi_{\mathcal{C}}^{E}$
has a similar expression in terms of spin operators.

From the point of view of measurement, it is more convenient to consider
another quantity, which has the same behavior as $\chi_{\mathcal{C}}^{E}$ in
the three phases, defined as:

$\tilde{\chi}_{\mathcal{C}}^{E}\equiv\frac{\absolutevalue{\expectationvalue{\prod_{i\in\mathcal{C}}\hat{S}_{i}^{x}}}}{\sqrt{\absolutevalue{\expectationvalue{\prod_{i\in\mathcal{L}}\hat{S}_{i}^{x}}}}},$
(62)

In the transverse-field-polarized (Higgs) phase,
$\tilde{\chi}_{\mathcal{C}}^{E}$ approaches a nonzero constant, just like
$\chi_{\mathcal{C}}^{E}$. Now, we argue that even in the QSL and confined
phases, $\tilde{\chi}_{\mathcal{C}}^{E}$ and $\chi_{\mathcal{C}}^{E}$ have the
same behavior. For a state $\ket{\Psi}$ that dominantly lies in the ice
manifold, with corrections from outside the ice manifold being of order
$\Omega/V$ (such as the ground state $\ket{\Psi_{g}}$), we have

$\langle\Psi|\hat{S}_{1}^{+}\hat{S}_{2}^{-}\hat{S}_{3}^{+}\cdots+\text{H.c.}|\Psi\rangle\\\
=\langle\Psi|(2\hat{S}_{1}^{x})(2\hat{S}_{2}^{x})(2\hat{S}_{3}^{x})\cdots|\Psi\rangle+\Theta\left((\Omega/V)^{L}\right),$
(63)

where $L$ is the number of sites on $\mathcal{C}$. The correction is of order
$(\Omega/V)^{L}$ by an argument similar to the one used to show that the error
is sixth order in the protocol to measure the plaquette X correlator (see
Appendix B). Thus, for small $\Omega/V$, $\chi_{\mathcal{C}}^{E}$ and
$\tilde{\chi}_{\mathcal{C}}^{E}$ are equal up to order $(\Omega/V)^{L}$.

The numerator and the denominator of $\tilde{\chi}_{\mathcal{C}}^{E}$ can be
measured by applying $\pi/2$ pulses about the $y$-axis and measuring, from the
snapshots, products of $\hat{S}^{z}$ along $\mathcal{C}$ and $\mathcal{L}$.
This procedure is similar to the protocol to measure the plaquette X
correlator, described in Sec. IV.1.1.

The operator $e^{i\sum_{\mathcal{C}}\hat{a}_{\mathbf{r}\mathbf{r}^{\prime}}}$
creates two opposite “electric charges” at the endpoints of $\mathcal{C}$. So
a magnetic analogue of $\chi^{E}_{\mathcal{C}}$ can also be defined, where the
numerator is the expectation value of the operator that creates a monopole and
an antimonopole at the endpoints of $\mathcal{C}$. Such an order parameter,
$\chi_{\mathcal{C}}^{M}$, detects long-range order in the monopole string
operator and is given by

$\chi_{\mathcal{C}}^{M}=\frac{\expectationvalue{\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}\left(\bm{\mathsf{r}}_{1}\xrightarrow{\mathcal{C}}\bm{\mathsf{r}}_{2}\right)}}{\sqrt{\expectationvalue{\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}\left(\bm{\mathsf{r}}_{1}\xrightarrow{\mathcal{L}}\bm{\mathsf{r}}_{1}\right)}}},$
(64)

where
$\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}(\bm{\mathsf{r}}_{1}\xrightarrow{\mathcal{C}}\bm{\mathsf{r}}_{2})$
inserts a monopole-antimonopole string along $\mathcal{C}$ and was defined in
Eq. (57). In this section, we use the notation where the path of the monopole-
antimonopole string is explicitly written in the argument of
$\hat{\mathcal{M}}^{\dagger}\hat{\mathcal{M}}$. Since this operator is
diagonal in the $\hat{S}^{z}$ basis, it can be measured straightforwardly from
the snapshots of the Rydberg-atom array.

In the confined phase, monopoles are condensed, so $\chi_{\mathcal{C}}^{M}$
should be a nonzero constant. In the deconfined phase, by the argument of Ref.
[89], the numerator of Eq. (64) decays to zero faster than the denominator as
the length of $\mathcal{C}$ increases. Therefore, in the deconfined phase,
$\chi_{\mathcal{C}}^{M}$ goes to zero as the length of $\mathcal{C}$
increases. In the Higgs phase, even though there is no long-range order in the
monopole string and both the numerator and denominator go to zero, by the
argument in Ref. [89], the ratio (i.e. $\chi_{\mathcal{C}}^{M}$) approaches a
nonzero constant as the length of $\mathcal{C}$ increases. But distinguishing
this non-zero constant from zero in finite-size numerics and experiment may be
challenging (similar to the situation for $\chi_{\mathcal{C}}^{E}$ in the
confined phase).

The behavior of the BFFM order parameters in various phases is summarized in
Table 2.

Before proceeding, we note that our protocols to measure the plaquette
correlators and the BFFM order parameter $\chi_{\mathcal{C}}^{E}$ work in the
limit $\Omega/V\ll 1$, which is outside the window in which the ground state
of Hamiltonian (6) is a QSL. However, we explained in Sec. III.3 that it is
possible to dynamically prepare finite puddles of QSL regions even in the
$\Omega/V\ll 1$ limit when the ground state is not a QSL. Our protocols can
then be applicable.

###  IV.4 Two-point $\hat{S}^{z}$ correlator

Consider two spins $\hat{S}^{z}_{\mathbf{r},\mu}$ and
$\hat{S}^{z}_{\mathbf{r}^{\prime},\nu}$ located on the sites
$\mathbf{r}+\mathbf{e}_{\mu}/2$ and $\mathbf{r}^{\prime}+\mathbf{e}_{\nu}/2$,
where $\mathbf{r}$ and $\mathbf{r}^{\prime}$ are the centers of two up-
pointing tetrahedra and $\mu,\nu\in\\{0,1,2,3\\}$ label the sites of the
tetrahedra (see Fig. 14). From the mapping of spins to gauge theory, Eqs. (13)
and (15), it can be seen that the two-point correlator of these two spins
$\langle\hat{S}^{z}_{\mathbf{r},\mu}\hat{S}^{z}_{\mathbf{r}^{\prime},\nu}\rangle$
is the same as the two-point correlator of the electric field.

The effective theory in the deconfined phase is the Maxwell electromagnetism.
In the $3+1$D continuum Maxwell electromagnetism, the correlator of the
Cartesian components of the electric field $\hat{e}_{\mathbf{r},i}$ for
$i\in\\{x,y,z\\}$ can be expressed as Gaussian integral which evaluate to [48]

$\langle\hat{e}_{\bm{0},i}\hat{e}_{\mathbf{R},j}\rangle_{0}\propto\frac{1}{R^{4}}\left(2\frac{R_{i}R_{j}}{R^{2}}-\delta_{ij}\right),$
(65)

where $\langle\cdot\rangle_{0}$ denotes expectation value with respect to the
Maxwell action.

The correlator of the electric field operators $\hat{e}_{\mathbf{r},\mu}$ for
$\mu\in\\{0,1,2,3\\}$ along the links of the diamond lattice are obtained from
Eq. (65) by taking components of the Cartesian electric field along the
vectors $\mathbf{e}_{\mu}$. The result is

Figure 14: Notation for the two-point $\hat{S}^{z}$ correlator. $\mathbf{r}$
and $\mathbf{r}^{\prime}$ are the positions of the centers of the tetrahedra.
$\mathbf{e}_{\mu}$ are the vectors joining the center of an up-pointing
tetrahedron to the centers of its neighboring down-pointing tetrahedra.

$\langle\hat{S}^{z}_{\mathbf{r},\mu}\hat{S}^{z}_{\mathbf{r^{\prime}},\nu}\rangle=\sum_{k,l\in\\{x,y,z\\}}(\mathbf{e}_{\mu})_{k}(\mathbf{e}_{\nu})_{l}\langle\hat{e}_{\mathbf{r},k}\hat{e}_{\mathbf{r}^{\prime},l}\rangle_{0},$
(66)

In the confined phase (ice ferromagnet), which is primarily diagonal in the
$\hat{S}^{z}$ basis, this correlator should approach a constant for large
$R=|\mathbf{r}-\mathbf{r^{\prime}}|$. On the other hand, in the TFP phase,
which is primarily a product state in the $\hat{S}^{x}$ basis, this correlator
decays at least exponentially with $R$ (see Table 2). Since
$\hat{S}^{z}_{\mathbf{r},\mu}\hat{S}^{z}_{\mathbf{r^{\prime}},\nu}$ is a
diagonal operator, its correlator can be measured experimentally by capturing
snapshots of the Rydberg-atom array and averaging over them.

##  V Discussion

In this work, we have presented a proposal to prepare and detect the
deconfined phase of the $U(1)$ gauge theory in 3+1 dimensions on a Rydberg
atom simulator. We first showed that laser-driven neutral atoms trapped in a
pyrochlore lattice using optical tweezer arrays naturally realise a $U(1)$
quantum spin liquid as the ground state when the laser detuning lies in a
specified window and the interactions between Rydberg atoms are restricted to
nearest-neighbor. We then studied the effect of van der Waals interactions
beyond nearest-neighbor. In the classical limit obtained by dropping the Rabi
frequency term, we showed that long-range interactions break the degeneracy to
select an ice ferromagnet as the ground state. We then studied the competition
between the long-ranged interactions that prefer an ordered state and quantum
fluctuations that prefer a QSL state, by calculating the energies in ansatz
wavefunctions using perturbation theory. We found that, for Rabi frequencies
greater than $\Omega_{C}\approx 0.44V$, the ground state is a QSL within our
approximation. When $\Omega$ is increased further, we argued that the QSL goes
into a transverse-field-polarized state via a Higgs transition. While we have
focused on the ground state, we also commented on the effect of dynamical
state preparation in deciding the nature of the prepared state. We then
provided experimental protocols for measuring the plaquette correlators,
Bricmont-Frölich-Fredenhagen-Marcu order parameters, the monopole-monopole
correlator, and the “electric field” correlator that can distinguish a QSL
phase from ordered phases.

Our ground state phase diagram is the result of an approximate calculation.
While it is possible that the true phase diagram differs from what we found,
we note that there are other knobs one can tune to get a desired phase
diagram. Dressed states created from multiple Rydberg and possibly ground
levels can be used to customize the interactions away from the isotropic
$1/r^{6}$ form we considered in this paper [94, 95, 96, 97, 98, 99]. Designing
a dressing scheme compatible with the symmetries of the pyrochlore lattice and
exploring the resulting phase diagrams is an interesting direction for future
work. We also note that our proposal requires two Rydberg excitations per
tetrahedron, meaning that it lies outside of the Rydberg-blockade regime and
is therefore sensitive to imperfections and thermal fluctuations in nearest-
neighbor spacing. It will therefore be useful to extend our proposal to the
blockade regime of one excitation per tetrahedron. While previous numerical
work on dimer models have required a nonzero RK potential (6-body term) to
achieve this, it will be worthwhile to study if one can engineer long-range
Rydberg interactions that stabilize a spin-liquid in the blockade regime.

Figure 15: A lattice made of corner-sharing tetrahedra different from the
pyrochlore lattice. The lattice consists of $ABAB\ldots$ stacking of the blue
($A$) and the orange ($B$) layers. A configuration satisfying
$n_{\mathchoice{\includegraphics[height=3.01389pt]{tetrahedron}}{\includegraphics[height=5.12358pt]{tetrahedron}}{\includegraphics[height=4.52083pt]{tetrahedron}}{\includegraphics[height=1.50694pt]{tetrahedron}}}=2$
is shown here.

One can also look for other lattices that could realize a $U(1)$ QSL ground
state. One such possibility is a lattice of corner-sharing tetrahedra where
all up-pointing tetrahedra (and separately all down-pointing tetrahedra) form
a hexagonal close-packed lattice shown in Fig. 15. If only nearest-neighbor
interactions are considered between atoms positioned on the sites of this
lattice, then, by perturbation theory in $\Omega/V$ for a particular range of
detunings, one gets ring exchange terms similar to the ones obtained in Sec.
II.1, and the system maps onto a dimer model. It is not known if this dimer
model is in the QSL phase when the RK potential is zero and long range van der
Waals interactions are included. Another open problem is to construct lattices
where a dimer model can be realized within the blockade regime without the RK
potential.

Next, we note that, formally, a distinction between the confined and
deconfined phases exists only in the thermodynamic limit. Experimentally,
there are two finiteness effects that can be important. First, a realistic
three dimensional Rydberg array will likely have a relatively small linear
dimension. Some of the correlators presented in Sec. IV require asymptotic
behavior in distance to distinguish different phases. Second, as found in Ref.
[81] and mentioned in Sec. III.3, a finite-time state preparation scheme would
generically prepare puddles of spin-liquid regions as opposed to an entire
spin liquid. It is therefore necessary to quantitatively study how the
behavior of the correlators is modified under these conditions.

We also note that, to translate field-theory observables into microscopic
variables, we relied on the perturbative limit of small $\Omega/V$. However in
the phase diagram that we found, the region where the spin liquid is a ground
state does not satisfy $\Omega/V\ll 1$. Understanding how the field-theory
operators (e.g. plaquette, monopole, and electric-field operators) get
renormalized away from the perturbative limit is important both from
fundamental and practical standpoints.

Our work is a proposal to prepare a gapless $U(1)$ spin liquid using unitary
evolution. An interesting research direction would be to come up with schemes
that also use projective measurements to expedite the state preparation along
the lines of Refs. [44, 100]. One can also explore how other exotic phases of
matter such as fractons and 3+1D topological order can potentially be realized
on a Rydberg simulator.

###### Acknowledgements.

We thank Nikita Astrakhantsev, Peter Lunts, Nathan Schine, Alexander
Schuckert, Dayal Singh, and Ashvin Vishwanath for discussions. J.S., G.N., and
V.G. were supported by NSF DMR-2037158, US-ARO Contract No.W911NF1310172, and
Simons Foundation. J.S. and A.V.G. were supported in part by AFOSR, NSF QLCI
(award No. OMA-2120757), DoE QSA, DoE ASCR Accelerated Research in Quantum
Computing program (award No. DE-SC0020312), the DoE ASCR Quantum Testbed
Pathfinder program (award No. DE-SC0019040), NSF PFCQC program, ARO MURI,
AFOSR MURI, and DARPA SAVaNT ADVENT.

## References

  * Anderson [1973] P. W. Anderson, Resonating valence bonds: A new kind of insulator?, Mater. Res. Bulletin 8, 153 (1973).
  * Savary and Balents [2017a] L. Savary and L. Balents, Quantum spin liquids: A review, Rep. Prog. Phys. 80, 016502 (2017a).
  * Senthil  _et al._ [2004] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. Fisher, Deconfined quantum critical points, Science 303, 1490 (2004).
  * Hermele  _et al._ [2005] M. Hermele, T. Senthil, and M. P. Fisher, Algebraic spin liquid as the mother of many competing orders, Physical Review B 72, 104404 (2005).
  * Song  _et al._ [2019] X.-Y. Song, C. Wang, A. Vishwanath, and Y.-C. He, Unifying description of competing orders in two-dimensional quantum magnets, Nature communications 10, 1 (2019).
  * Nambiar  _et al._ [2022] G. Nambiar, D. Bulmash, and V. Galitski, Monopole josephson effects in a dirac spin liquid (2022).
  * Marciano and Pagels [1978] W. Marciano and H. Pagels, Quantum chromodynamics, Phys. Rep. 36, 137 (1978).
  * Meyer-Ortmanns [1996] H. Meyer-Ortmanns, Phase transitions in quantum chromodynamics, Rev. Mod. Phys. 68, 473 (1996).
  * Kogut [1983] J. B. Kogut, The lattice gauge theory approach to quantum chromodynamics, Rev. Mod. Phys. 55, 775 (1983).
  * O’Raifeartaigh and Straumann [2000] L. O’Raifeartaigh and N. Straumann, Gauge theory: Historical origins and some modern developments, Rev. Mod. Phys. 72, 1 (2000).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
  * Kitaev [2003] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Ann. Phys. (N. Y.) 303, 2 (2003).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. D. Sarma, Non-abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Knolle and Moessner [2019] J. Knolle and R. Moessner, A field guide to spin liquids, Annu. Rev. Condens. Matter Phys. 10, 451 (2019).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, A rydberg quantum simulator, Nat. Phys. 6, 382 (2010).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, _et al._ , Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227 (2021).
  * Wu  _et al._ [2021] X. Wu, X. Liang, Y. Tian, F. Yang, C. Chen, Y.-C. Liu, M. K. Tey, and L. You, A concise review of rydberg atom based quantum computation and quantum simulation, Chin. Phys. B 30, 020305 (2021).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with Rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Lukin  _et al._ [2001] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole Blockade and Quantum Information Processing in Mesoscopic Atomic Ensembles, Phys. Rev. Lett. 87, 037901 (2001).
  * Labuhn  _et al._ [2016] H. Labuhn, D. Barredo, S. Ravets, S. de Léséleuc, T. Macrì, T. Lahaye, and A. Browaeys, Tunable two-dimensional arrays of single Rydberg atoms for realizing quantum Ising models, Nature 534, 667 (2016).
  * Jaksch  _et al._ [2000] D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Fast Quantum Gates for Neutral Atoms, Phys. Rev. Lett. 85, 2208 (2000).
  * Browaeys  _et al._ [2016] A. Browaeys, D. Barredo, and T. Lahaye, Experimental investigations of dipole–dipole interactions between a few Rydberg atoms, J. Phys. B 49, 152001 (2016).
  * Barredo  _et al._ [2016] D. Barredo, S. de Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, An atom-by-atom assembler of defect-free arbitrary two-dimensional atomic arrays, Science 354, 1021 (2016).
  * Béguin  _et al._ [2013] L. Béguin, A. Vernier, R. Chicireanu, T. Lahaye, and A. Browaeys, Direct Measurement of the van der Waals Interaction between Two Rydberg Atoms, Phys. Rev. Lett. 110, 263201 (2013).
  * Barredo  _et al._ [2018] D. Barredo, V. Lienhard, S. de Léséleuc, T. Lahaye, and A. Browaeys, Synthetic three-dimensional atomic structures assembled atom by atom, Nature 561, 79 (2018).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletic, and M. D. Lukin, Probing Topological Spin Liquids on a Programmable Quantum Simulator, Science 374, 1242 (2021).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Bluvstein  _et al._ [2022] D. Bluvstein, H. Levine, G. Semeghini, T. T. Wang, S. Ebadi, M. Kalinowski, A. Keesling, N. Maskara, H. Pichler, M. Greiner, V. Vuletić, and M. D. Lukin, A quantum processor based on coherent transport of entangled atom arrays, Nature 604, 451 (2022).
  * Guardado-Sanchez  _et al._ [2018] E. Guardado-Sanchez, P. T. Brown, D. Mitra, T. Devakul, D. A. Huse, P. Schauß, and W. S. Bakr, Probing the Quench Dynamics of Antiferromagnetic Correlations in a 2D Quantum Ising Spin System, Phys. Rev. X 8, 021069 (2018).
  * Keesling  _et al._ [2019] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum Kibble–Zurek mechanism and critical dynamics on a programmable Rydberg simulator, Nature 568, 207 (2019).
  * Levine  _et al._ [2019] H. Levine, A. Keesling, G. Semeghini, A. Omran, T. T. Wang, S. Ebadi, H. Bernien, M. Greiner, V. Vuletić, H. Pichler, and M. D. Lukin, Parallel Implementation of High-Fidelity Multiqubit Gates with Neutral Atoms, Phys. Rev. Lett. 123, 170503 (2019).
  * Lienhard  _et al._ [2018] V. Lienhard, S. de Léséleuc, D. Barredo, T. Lahaye, A. Browaeys, M. Schuler, L.-P. Henry, and A. M. Läuchli, Observing the Space- and Time-Dependent Growth of Correlations in Dynamically Tuned Synthetic Ising Models with Antiferromagnetic Interactions, Phys. Rev. X 8, 021070 (2018).
  * Madjarov  _et al._ [2020] I. S. Madjarov, J. P. Covey, A. L. Shaw, J. Choi, A. Kale, A. Cooper, H. Pichler, V. Schkolnik, J. R. Williams, and M. Endres, High-fidelity entanglement and detection of alkaline-earth Rydberg atoms, Nat. Phys. 16, 857 (2020).
  * Omran  _et al._ [2019] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Generation and manipulation of Schrödinger cat states in Rydberg atom arrays, Science 365, 570 (2019).
  * Schauß  _et al._ [2015] P. Schauß, J. Zeiher, T. Fukuhara, S. Hild, M. Cheneau, T. Macrì, T. Pohl, I. Bloch, and C. Gross, Crystallization in Ising quantum magnets, Science 347, 1455 (2015).
  * Song  _et al._ [2021] Y. Song, M. Kim, H. Hwang, W. Lee, and J. Ahn, Quantum annealing of Cayley-tree Ising spins at small scales, Phys. Rev. Res. 3, 013286 (2021).
  * Verresen  _et al._ [2021a] R. Verresen, M. D. Lukin, and A. Vishwanath, Prediction of Toric Code Topological Order from Rydberg Blockade, Phys. Rev. X 11, 031005 (2021a).
  * Samajdar  _et al._ [2021] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Quantum phases of Rydberg atoms on a kagome lattice, Proc. Natl. Acad. Sci. U.S.A. 118, e2015785118 (2021).
  * Myerson-Jain  _et al._ [2022] N. E. Myerson-Jain, S. Yan, D. Weld, and C. Xu, Construction of fractal order and phase transition with rydberg atoms, Phys. Rev. Lett. 128, 017601 (2022).
  * Slagle  _et al._ [2022] K. Slagle, Y. Liu, D. Aasen, H. Pichler, R. S. K. Mong, X. Chen, M. Endres, and J. Alicea, Quantum spin liquids bootstrapped from ising criticality in rydberg arrays, Phys. Rev. B 106, 115122 (2022).
  * Kornjača  _et al._ [2022] M. Kornjača, R. Samajdar, T. Macrì, N. Gemelke, S.-T. Wang, and F. Liu, Trimer quantum spin liquid in a honeycomb array of rydberg atoms, arXiv preprint arXiv:2211.00653 (2022).
  * Zhou  _et al._ [2022] Z. Zhou, Z. Yan, C. Liu, Y. Chen, and X.-F. Zhang, Quantum simulation of two-dimensional u(1) gauge theory in rydberg atom arrays, arXiv preprint arXiv:2212.10863 (2022).
  * Liu  _et al._ [2022] Y.-J. Liu, K. Shtengel, A. Smith, and F. Pollmann, Methods for simulating string-net states and anyons on a digital quantum computer, PRX Quantum 3, 040315 (2022).
  * Tantivasadakarn  _et al._ [2022] N. Tantivasadakarn, R. Verresen, and A. Vishwanath, The shortest route to non-abelian topological order on a quantum processor, arXiv preprint arXiv:2209.03964 (2022).
  * Verresen  _et al._ [2021b] R. Verresen, N. Tantivasadakarn, and A. Vishwanath, Efficiently preparing ghz, topological and fracton states by measuring cold atoms, arXiv preprint arXiv:2112.03061 (2021b).
  * Polyakov [1975] A. M. Polyakov, Compact gauge fields and the infrared catastrophe, Phys. Lett. B 59, 82 (1975).
  * Polyakov [1977] A. Polyakov, Quark confinement and topology of gauge theories, Nucl. Phys. B 120, 429 (1977).
  * Hermele  _et al._ [2004] M. Hermele, M. P. A. Fisher, and L. Balents, Pyrochlore photons: The U(1) spin liquid in a S = 1/2 three-dimensional frustrated magnet, Phys. Rev. B 69, 064404 (2004).
  * Huse  _et al._ [2003] D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi, Coulomb and Liquid Dimer Models in Three Dimensions, Phys. Rev. Lett. 91, 167004 (2003).
  * Banerjee  _et al._ [2008] A. Banerjee, S. V. Isakov, K. Damle, and Y. B. Kim, Unusual Liquid State of Hard-Core Bosons on the Pyrochlore Lattice, Phys. Rev. Lett. 100, 047208 (2008).
  * Shannon  _et al._ [2012] N. Shannon, O. Sikora, F. Pollmann, K. Penc, and P. Fulde, Quantum Ice: A Quantum Monte Carlo Study, Phys. Rev. Lett. 108, 067204 (2012).
  * Bramwell and Harris [2020] S. T. Bramwell and M. J. Harris, The history of spin ice, J. Condens. Matter Phys. 32, 374010 (2020).
  * Gingras [2011] M. J. Gingras, Spin ice, Introduction to frustrated magnetism , 293 (2011).
  * Udagawa and Jaubert [2021] M. Udagawa and L. Jaubert,  _Spin Ice_ (Springer, 2021).
  * Skjærvø  _et al._ [2020] S. H. Skjærvø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, Advances in artificial spin ice, Nat. Rev. Phys. 2, 13 (2020).
  * Nisoli  _et al._ [2013] C. Nisoli, R. Moessner, and P. Schiffer, Colloquium: Artificial spin ice: Designing and imaging magnetic frustration, Rev. Mod. Phys. 85, 1473 (2013).
  * Pauling [1935] L. Pauling, The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic Arrangement, J. Am. Chem. Soc. 57, 2680 (1935).
  * Gingras and McClarty [2014] M. J. Gingras and P. A. McClarty, Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets, Rep. Prog. Phys. 77, 056501 (2014).
  * Tewari  _et al._ [2006] S. Tewari, V. W. Scarola, T. Senthil, and S. D. Sarma, Emergence of artificial photons in an optical lattice, Phys. Rev. Lett. 97, 200401 (2006).
  * Glaetzle  _et al._ [2014] A. W. Glaetzle, M. Dalmonte, R. Nath, I. Rousochatzakis, R. Moessner, and P. Zoller, Quantum spin-ice and dimer models with rydberg atoms, Phys. Rev. X 4, 041037 (2014).
  * Celi  _et al._ [2020a] A. Celi, B. Vermersch, O. Viyuela, H. Pichler, M. D. Lukin, and P. Zoller, Emerging two-dimensional gauge theories in rydberg configurable arrays, Phys. Rev. X 10, 021057 (2020a).
  * Emonts and Wessel [2018] P. Emonts and S. Wessel, Monte carlo study of the discontinuous quantum phase transition in the transverse-field ising model on the pyrochlore lattice, Phys. Rev. B 98, 174433 (2018).
  * Sikora  _et al._ [2009] O. Sikora, F. Pollmann, N. Shannon, K. Penc, and P. Fulde, Quantum Liquid with Deconfined Fractional Excitations in Three Dimensions, Phys. Rev. Lett. 103, 247001 (2009).
  * Anderson [1956] P. W. Anderson, Ordering and antiferromagnetism in ferrites, Phys. Rev. 102, 1008 (1956).
  * Rokhsar and Kivelson [1988] D. S. Rokhsar and S. A. Kivelson, Superconductivity and the Quantum Hard-Core Dimer Gas, Phys. Rev. Lett. 61, 2376 (1988).
  * Balasubramanian  _et al._ [2022a] S. Balasubramanian, D. Bulmash, V. Galitski, and A. Vishwanath, Exact wavefunction dualities and phase diagrams of 3d quantum vertex models (2022a).
  * Balasubramanian  _et al._ [2022b] S. Balasubramanian, V. Galitski, and A. Vishwanath, Classical vertex model dualities in a family of two-dimensional frustrated quantum antiferromagnets, Phys. Rev. B 106, 195127 (2022b).
  * Savary and Balents [2017b] L. Savary and L. Balents, Disorder-Induced Quantum Spin Liquid in Spin Ice Pyrochlores, Phys. Rev. Lett. 118, 087203 (2017b).
  * Chern  _et al._ [2010] C.-H. Chern, C.-N. Liao, and Y.-Z. Chou, Disorder from disorder and confinement in the quantum Ising model in the pyrochlore lattice, arXiv:1003.4204 [cond-mat] (2010).
  * Wilson [1974] K. G. Wilson, Confinement of quarks, Phys. Rev. D 10, 2445 (1974).
  * Polyakov [1978] A. M. Polyakov, Thermal Properties of Gauge fields and Quark Liberation (1978).
  * Kogut [1979] J. B. Kogut, An introduction to lattice gauge theory and spin systems, Rev. Mod. Phys. 51, 659 (1979).
  * Celi  _et al._ [2020b] A. Celi, B. Vermersch, O. Viyuela, H. Pichler, M. D. Lukin, and P. Zoller, Emerging Two-Dimensional Gauge Theories in Rydberg Configurable Arrays, Phys. Rev. X 10, 021057 (2020b).
  * Röchner  _et al._ [2016] J. Röchner, L. Balents, and K. P. Schmidt, Spin liquid and quantum phase transition without symmetry breaking in a frustrated three-dimensional ising model, Phys. Rev. B 94, 201111 (2016).
  * Pace  _et al._ [2021] S. D. Pace, S. C. Morampudi, R. Moessner, and C. R. Laumann, Emergent Fine Structure Constant of Quantum Spin Ice Is Large, Phys. Rev. Lett. 127, 117205 (2021).
  * Chen [2016] G. Chen, “Magnetic monopole” condensation of the pyrochlore ice $U(1)$ quantum spin liquid: Application to ${\mathrm{Pr}}_{2}{\mathrm{Ir}}_{2}{\mathrm{O}}_{7}$ and ${\mathrm{Yb}}_{2}{\mathrm{Ti}}_{2}{\mathrm{O}}_{7}$, Phys. Rev. B 94, 205107 (2016).
  * Slagle and Kim [2017] K. Slagle and Y. B. Kim, Fracton Topological Order from Nearest-Neighbor Two-Spin Interactions and Dualities, Phys. Rev. B 96, 165106 (2017).
  * Melko  _et al._ [2001] R. G. Melko, B. C. den Hertog, and M. J. P. Gingras, Long-range order at low temperatures in dipolar spin ice, Phys. Rev. Lett. 87, 067203 (2001).
  * Melko and Gingras [2004] R. G. Melko and M. J. Gingras, Monte carlo studies of the dipolar spin ice model, J. Condens. Matter Phys. 16, R1277 (2004).
  * Savary and Balents [2012] L. Savary and L. Balents, Coulombic quantum liquids in spin-$1/2$ pyrochlores, Phys. Rev. Lett. 108, 037202 (2012).
  * Sahay  _et al._ [2022] R. Sahay, A. Vishwanath, and R. Verresen, Quantum spin puddles and lakes: Nisq-era spin liquids from non-equilibrium dynamics, arXiv preprint arXiv:2211.01381 (2022).
  * Zurek [1985] W. H. Zurek, Cosmological experiments in superfluid helium?, Nature 317, 505 (1985).
  * Kibble [1976] T. W. Kibble, Topology of cosmic domains and strings, J. Phys. A 9, 1387 (1976).
  * Del Campo and Zurek [2014] A. Del Campo and W. H. Zurek, Universality of phase transition dynamics: Topological defects from symmetry breaking, Int. J. Mod. Phys. A 29, 1430018 (2014).
  * Svetitsky [1986] B. Svetitsky, Symmetry aspects of finite-temperature confinement transitions, Phys. Rep. 132, 1 (1986).
  * Fradkin [2013] E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013).
  * Shankar [2017] R. Shankar,  _Quantum field theory and condensed matter: An introduction_ (Cambridge University Press, 2017).
  * Fredenhagen and Marcu [1983] K. Fredenhagen and M. Marcu, Charged states in $Z_{2}$ gauge theories, Commun. Math. Phys. 92, 81 (1983).
  * Fredenhagen and Marcu [1986] K. Fredenhagen and M. Marcu, Confinement criterion for qcd with dynamical quarks, Phys. Rev. Lett. 56, 223 (1986).
  * Fredenhagen and Marcu [1988] K. Fredenhagen and M. Marcu, Dual interpretation of order parameters for lattice gauge theories with matter fields, Nucl. Phys. B-Proceedings Supplements 4, 352 (1988).
  * Marcu [1986] M. Marcu, (uses of) an order parameter for lattice gauge theories with matter fields, in  _Lattice Gauge Theory: A Challenge in Large-Scale Computing_, edited by B. Bunk, K. H. Mütter, and K. Schilling (Springer US, Boston, MA, 1986) pp. 267–278.
  * Bricmont and Frölich [1983] J. Bricmont and J. Frölich, An order parameter distinguishing between different phases of lattice gauge theories with matter fields, Phys. Lett. B 122, 73 (1983).
  * Gregor  _et al._ [2011] K. Gregor, D. A. Huse, R. Moessner, and S. L. Sondhi, Diagnosing deconfinement and topological order, New J. Phys. 13, 025009 (2011).
  * de Léséleuc  _et al._ [2017] S. de Léséleuc, D. Barredo, V. Lienhard, A. Browaeys, and T. Lahaye, Optical control of the resonant dipole-dipole interaction between rydberg atoms, Phys. Rev. Lett. 119, 053202 (2017).
  * Glaetzle  _et al._ [2015] A. W. Glaetzle, M. Dalmonte, R. Nath, C. Gross, I. Bloch, and P. Zoller, Designing frustrated quantum magnets with laser-dressed rydberg atoms, Phys. Rev. Lett. 114, 173002 (2015).
  * van Bijnen and Pohl [2015] R. M. W. van Bijnen and T. Pohl, Quantum magnetism and topological ordering via rydberg dressing near förster resonances, Phys. Rev. Lett. 114, 243002 (2015).
  * Petrosyan and Mølmer [2014] D. Petrosyan and K. Mølmer, Binding potentials and interaction gates between microwave-dressed rydberg atoms, Phys. Rev. Lett. 113, 123003 (2014).
  * Graß  _et al._ [2018] T. Graß, P. Bienias, M. J. Gullans, R. Lundgren, J. Maciejko, and A. V. Gorshkov, Fractional quantum hall phases of bosons with tunable interactions: From the laughlin liquid to a fractional wigner crystal, Phys. Rev. Lett. 121, 253403 (2018).
  * Young  _et al._ [2021] J. T. Young, P. Bienias, R. Belyansky, A. M. Kaufman, and A. V. Gorshkov, Asymmetric blockade and multiqubit gates via dipole-dipole interactions, Phys. Rev. Lett. 127, 120501 (2021).
  * Lavasani  _et al._ [2022] A. Lavasani, Z.-X. Luo, and S. Vijay, Monitored quantum dynamics and the kitaev spin liquid, arXiv preprint arXiv:2207.02877 (2022).

##  Appendix A Gauge mean field theory

In this appendix, we first provide details of the gauge mean field theory
calculation sketched in Sec. III.2.1, with a focus on capturing the Higgs
transition. Then, we attempt to use the same technique in the small-$\Omega$
limit to obtain the confinement-deconfinement transition. We find that, in
this limit, the technique is fraught with a serious limitation stemming from
neglecting gauge fluctuations.

Starting from Eq. (44) of Sec. III.2.1 and performing the mean-field
decoupling, we get

$\displaystyle\hat{H}_{\text{MF}}$
$\displaystyle=\hat{H}_{\Phi}+\hat{H}_{\mathsf{s}}+\hat{H}_{\text{c}},\text{
where}$ (67) $\displaystyle\hat{H}_{\Phi}=$
$\displaystyle\frac{V}{2}\sum_{\mathbf{r}\in
A,B}\hat{Q}^{2}_{\mathbf{r}}-\frac{\Omega}{2}\sum_{(\mathbf{r}\in
A),\mu}\left(\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}}\expectationvalue{\hat{\mathsf{s}}^{+}_{\mathbf{r},\mu}}+\text{
H.c.}\right),$ $\displaystyle\hat{H}_{\mathsf{s}}=$
$\displaystyle-\frac{\Omega}{2}\sum_{(\mathbf{r}\in
A),\mu}\left(\expectationvalue{\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}}}\hat{\mathsf{s}}^{+}_{\mathbf{r},\mu}+\text{
H.c.}\right)$ $\displaystyle+\sum_{(\mathbf{r}\in
A),\mu}\hat{\mathsf{s}}^{z}_{\mathbf{r},\mu}\sum_{(\mathbf{r}^{\prime}\in
A),\nu}\left(V_{\mu\nu}(\mathbf{r}-\mathbf{r}^{\prime})\expectationvalue{\hat{\mathsf{s}}^{z}_{\mathbf{r}^{\prime},\nu}}\right),$
$\displaystyle\hat{H}_{\text{c}}=$
$\displaystyle\frac{\Omega}{2}\sum_{(\mathbf{r}\in
A),\mu}\left(\expectationvalue{\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}}}\expectationvalue{\hat{\mathsf{s}}^{+}_{\mathbf{r},\mu}}+\text{
H.c.}\right)$ $\displaystyle-\frac{1}{2}\sum_{(\mathbf{r}\in
A),\mu}\hat{\mathsf{s}}^{z}_{\mathbf{r},\mu}\sum_{(\mathbf{r}^{\prime}\in
A),\nu}\left(V_{\mu\nu}(\mathbf{r}-\mathbf{r}^{\prime})\expectationvalue{\hat{\mathsf{s}}^{z}_{\mathbf{r}^{\prime},\nu}}\right).$

$\hat{H}_{\text{c}}$ is a constant, and
$V_{\mu\nu}(\mathbf{r}-\mathbf{r^{\prime}})$ was defined in Sec. III.2.1.
$\hat{H}_{\mathsf{s}}$ above is of the form $-\sum_{(\mathbf{r}\in
A),\mu}\left(h^{x}_{\mathbf{r},\mu}\hat{\mathsf{s}}^{x}_{\mathbf{r},\mu}+h^{z}_{\mathbf{r},\mu}\hat{\mathsf{s}}^{z}_{\mathbf{r},\mu}\right)$,
where

$\displaystyle h^{x}_{\mathbf{r},\mu}$
$\displaystyle=\Omega\expectationvalue{\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}}},$
(68) $\displaystyle h^{z}_{\mathbf{r},\mu}$
$\displaystyle=-\sum_{(\mathbf{r}^{\prime}\in
A),\nu}\left(V_{\mu\nu}(\mathbf{r}-\mathbf{r}^{\prime})\expectationvalue{\hat{\mathsf{s}}^{z}_{\mathbf{r}^{\prime},\nu}}\right),$

and
$\expectationvalue{\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}}}$
is calculated in the ground state of $\hat{H}_{\Phi}$, which in turn depends
on $\expectationvalue{\hat{\mathsf{s}}^{+}}$. (We have implicitly assumed here
that
$\expectationvalue{\hat{\Phi}_{\mathbf{r}}^{\dagger}\hat{\Phi}_{\mathbf{r}+\mathbf{e}_{\mu}}}$
is real, which we will show can be assumed self-consistently.) This implies
that, in the ground state,

$\expectationvalue{\hat{\mathsf{s}}^{i}_{\mathbf{r},\mu}}=\frac{h^{i}_{\mathbf{r},\mu}}{2\absolutevalue{\mathbf{h}_{\mathbf{r},\mu}}}\text{
for }i=x,z.$ (69)

Our goal is to self-consistently minimize the ground-state energy of the mean-
field Hamiltonian subject to the constraints in Eqs. (42) and (43). We showed
in Sec. III.1.2 that the ordered ground state at $\Omega=0$ has momentum
$\mathbf{k}=\bm{0}$. Also, the TFP state in the large-$\Omega$ limit is a
$\mathbf{k}=\bm{0}$ state. So we start with a mean-field ansatz with full
translation symmetry (similar to Ref. [80]):

$\displaystyle\expectationvalue{\mathsf{s}^{+}_{\mathbf{r},\mu}}$
$\displaystyle=\frac{1}{2}\cos\theta,$ (70)
$\displaystyle\expectationvalue{\mathsf{s}^{z}_{\mathbf{r},\mu}}$
$\displaystyle=\frac{1}{2}\varepsilon_{\mu}\sin\theta,$

where $\varepsilon_{\mu}=1,1,-1,-1$ for $\mu=0,1,2,3$, respectively. To solve
the matter sector, it is convenient to deal with the Lagrangian instead of the
Hamiltonian. The imaginary-time Lagrangian for the matter sector is

$\displaystyle\mathcal{L}$ $\displaystyle=\frac{1}{2V}\sum_{\mathbf{r}\in
A,B}\absolutevalue{(\partial_{\tau}-iv_{\mathbf{r}})\Phi_{\mathbf{r}}}^{2}$
(71) $\displaystyle-\frac{\Omega\cos\theta}{4}\sum_{(\mathbf{r}\in
A),\mu}(\Phi_{\mathbf{r}}^{*}\Phi_{\mathbf{r}+\mathbf{e}_{\mu}}e^{ia_{\mathbf{r},\mu}}+\text{
c.c.})$ $\displaystyle-i\sum_{\mathbf{r}\in
A,B}\left[\eta_{\mathbf{r}}v_{\mathbf{r}}\left(\sum_{\mu}\mathsf{s}^{z}_{\mathbf{r}+\eta_{\mathbf{r}}\mathbf{e}_{\mu}/2}\right)+\tilde{\lambda}_{\mathbf{r}}(\absolutevalue{\Phi_{\mathbf{r}}}^{2}-1)\right],$

where the Lagrange multiplier $\tilde{\lambda}_{\mathbf{r}}$ (which gets
integrated over) enforces the constraint
$\absolutevalue{\Phi_{\mathbf{r}}}^{2}=1$. The Lagrange multiplier
$v_{\mathbf{r}}$ enforces the constraint (43). To zeroth order, we ignore the
gauge fluctuation $a_{\mathbf{r},\mu}$. The matter Lagrangian alone, despite
being quadratic in the rotor variables, is nevertheless interacting because a
quadratic term in rotor operators is nonlinear in terms of canonical bosons
(in other words, it is a cosine term in the phase of the rotor.) In order to
make progress, Ref. [80] assumes that, at the saddle point,
$\tilde{\lambda}_{\mathbf{r}}$ takes on a spatially uniform and purely
imaginary value $i\lambda$, and also implicitly assumes that $v_{\mathbf{r}}$
is 0 at the saddle point. Here, we will follow suit while acknowledging that
these approximations are uncontrolled. Making these simplifications, we obtain

$\displaystyle\mathcal{L}=$
$\displaystyle\frac{1}{2V}\sum_{\mathbf{r}}\absolutevalue{\partial_{\tau}\Phi_{\mathbf{r}}}^{2}-\frac{\Omega\cos\theta}{4}\sum_{(\mathbf{r}\in
A),\mu}\left(\Phi_{\mathbf{r}}^{*}\Phi_{\mathbf{r}+\mathbf{e}_{\mu}}+\text{
c.c.}\right)$ (72)
$\displaystyle+\lambda\sum_{\mathbf{r}}(\absolutevalue{\Phi_{\mathbf{r}}}^{2}-1).$

The constraints now simplify to

$\displaystyle\expectationvalue{\Phi^{\dagger}_{\mathbf{r}}\Phi_{\mathbf{r}}}=1,$
(73) $\displaystyle
h^{x}=\Omega\expectationvalue{\Phi^{\dagger}_{\mathbf{r}}\Phi_{\mathbf{r}+\mathbf{e}_{\mu}}}.$
(74)

Now, we have a quadratic Lagrangian, which we solve by Fourier transformation.
Our Fourier transformation convention is (for $\alpha\in\\{A,B\\}$)

$\Phi_{\mathbf{r},\alpha}(\tau)=T\sum_{\omega_{n}}\sum_{\mathbf{k}\in
BZ}\Phi_{\mathbf{k},\alpha}(\omega_{n})e^{i(\mathbf{k}\cdot\mathbf{r}-\omega_{n}\tau)},$
(75)

where $T$ is the temperature, $\omega_{n}$ are Matsubara frequencies and we
eventually take the limit $T\to 0$. Eq. (72) becomes

$\mathcal{L}=T\sum_{\mathbf{k},\omega_{n}}\begin{pmatrix}\Phi^{*}_{\mathbf{k},A}(\omega_{n})&\Phi^{*}_{\mathbf{k},B}(\omega_{n})\end{pmatrix}\mathcal{G}^{-1}_{\mathbf{k}}(\omega_{n})\begin{pmatrix}\Phi_{\mathbf{k},A}(\omega_{n})\\\
\Phi_{\mathbf{k},B}(\omega_{n})\end{pmatrix},$ (76)

where

$\mathcal{G}^{-1}_{\mathbf{k}}(\omega_{n})=\begin{pmatrix}\frac{\omega_{n}^{2}}{2V}+\lambda&-\frac{\Omega\cos\theta}{4}f_{\mathbf{k}}\\\
-\frac{\Omega\cos\theta}{4}f^{*}_{\mathbf{k}}&\frac{\omega_{n}^{2}}{2V}+\lambda\end{pmatrix}.$
(77)

Here,

$f_{\mathbf{k}}=1+e^{-ik_{1}}+e^{-ik_{2}}+e^{-ik_{3}},$ (78)

where $\mathbf{k}\equiv
k_{1}\mathbf{b}_{1}+k_{2}\mathbf{b}_{2}+k_{3}\mathbf{b}_{3}$, and
$\mathbf{b}_{1}$, $\mathbf{b}_{2}$ and $\mathbf{b}_{3}$ are reciprocal lattice
vectors of the FCC lattice satisfying
$\mathbf{a}_{i}\cdot\mathbf{b}_{j}=\delta_{ij}$.

Upon inverting the above matrix, we find that the eigenvalues of
$\mathcal{G}_{\mathbf{k}}(\omega_{n})$ are
$\frac{2V}{\omega_{n}^{2}+\left(\omega^{\pm}_{\mathbf{k}}(\lambda,\theta)\right)^{2}}$,
where the dispersion of the two bosonic bands is

$\omega^{\pm}_{\mathbf{k}}(\lambda,\theta)=\sqrt{2V\left(\lambda\pm\frac{\Omega\cos\theta}{4}\absolutevalue{f_{\mathbf{k}}}\right)}.$
(79)

As long as the spinon dispersion is gapped, spinons will not condense. From
the dispersion above, we see that the dispersion becomes gapless when
$\lambda=\Omega\cos\theta$. However, as we will see below, for fixed $\theta$
and $\Omega$, $\lambda$ is determined by the constraint in Eq. (73). Therefore
the condition $\lambda=\Omega\cos\theta$ is met for a specific
$\Omega=\Omega^{\text{MF}}_{H}$, which we will calculate below. Before that,
will go through a few intermediate steps. First, the matrix form of
$\mathcal{G}_{\mathbf{k}}(\omega_{n})$ is (assuming $\Omega>0$)

$\begin{aligned}
\mathcal{G}_{\mathbf{k}}(\omega_{n})=V\begin{pmatrix}\frac{1}{\omega_{n}^{2}+(\omega^{+}_{\mathbf{k}})^{2}}+\frac{1}{\omega_{n}^{2}+(\omega^{-}_{\mathbf{k}})^{2}}&g_{\mathbf{k}}\left(\frac{1}{\omega_{n}^{2}+(\omega^{+}_{\mathbf{k}})^{2}}-\frac{1}{\omega_{n}^{2}+(\omega^{-}_{\mathbf{k}})^{2}}\right)\\\
g_{\mathbf{k}}^{*}\left(\frac{1}{\omega_{n}^{2}+(\omega^{+}_{\mathbf{k}})^{2}}-\frac{1}{\omega_{n}^{2}+(\omega^{-}_{\mathbf{k}})^{2}}\right)&\frac{1}{\omega_{n}^{2}+(\omega^{+}_{\mathbf{k}})^{2}}+\frac{1}{\omega_{n}^{2}+(\omega^{-}_{\mathbf{k}})^{2}}\end{pmatrix}\end{aligned},$
(80)

where

$g_{\mathbf{k}}=\begin{cases}-\frac{f_{\mathbf{k}}}{\absolutevalue{f_{\mathbf{k}}}}&\text{
when }0\leq\theta<\pi/2,\\\ 0&\text{ when }\theta=\pi/2.\end{cases}$ (81)

With the Green’s function in hand, we are now ready to impose the constraints,
Eq. (73) and Eq. (74). First, we calculate equal-time correlation functions of
$\Phi$ (by performing the Matsubara sum on the Green’s function). Using these,
the constraints in Eq. (73) and Eq. (74) become, respectively,

$\displaystyle F_{1}(\lambda,\theta)$
$\displaystyle\equiv\frac{V}{2N_{\text{u.c.}}}\sum_{\mathbf{k}}\left(\frac{1}{\absolutevalue{\omega^{+}_{\mathbf{k}}}}+\frac{1}{\absolutevalue{\omega^{-}_{\mathbf{k}}}}\right)=1,$
(82) $\displaystyle\Omega F_{2}(\lambda,\theta)$
$\displaystyle\equiv\Omega\frac{V}{2N_{\text{u.c.}}}\sum_{\mathbf{k}}g_{\mathbf{k}}\left(\frac{1}{\absolutevalue{\omega^{-}_{\mathbf{k}}}}-\frac{1}{\absolutevalue{\omega^{+}_{\mathbf{k}}}}\right)=h^{x}.$
(83)

Next, by imposing Eq. (69) with the help of Eq. (68), we get

$h^{z}=-\frac{\mathcal{B}\sin\theta}{2},\text{ where
}\mathcal{B}=\frac{\sin\theta}{2}\sum_{(\mathbf{r}^{\prime}\in
A),\nu}V_{0,\nu}(-\mathbf{r}^{\prime})\varepsilon_{\nu}.$ (84)

For a given $\theta$, Eq. (82) determines $\lambda$. We see that there are
three self-consistent solutions for $\theta$:

$\theta=\begin{cases}0,\\\ \pi/2,\\\ \cos^{-1}\left(\frac{2\Omega
F_{2}(\lambda,\theta)}{\mathcal{B}}\right).\end{cases}$ (85)

Within gMFT (gauge mean field theory), these three solutions correspond to a
QSL, a “Coulomb ferromagnet” (spin liquid with nonzero ice ferromagnetic order
parameter), and an ice ferromagnet, respectively [80]. For a fixed parameter
$\Omega$, the true solution depends on which of the three solutions above has
lower energy with respect to the mean-field Hamiltonian (67). Suppose that,
for large enough $\Omega$, one is in the QSL phase, i.e., $\theta=0$ and
$\expectationvalue{\hat{\Phi}_{\mathbf{r}}}=0$. Now, the bosons will condense
when their dispersion becomes gapless, i.e., $\lambda=\Omega$. Using
constraint (82), we find that this transition point is
$\Omega^{H}_{\text{MF}}\approx 0.7V$, as also found in Ref. [68]. For
$\Omega>\Omega^{H}_{\text{MF}}$, the ground state is in the TFP phase.

Having identified the Higgs transition point, we now attempt to identify the
confinement-deconfinement transition for low $\Omega$, i.e., find $\Omega$ at
which $\theta=0$ becomes the lowest-energy saddle-point. Using Eq. (67), we
get the following expression for the mean-field energy:

$E_{\text{MF}}=K-N_{\text{u.c.}}\left(2\Omega
F_{2}(\lambda,\theta)\cos\theta+\frac{\mathcal{B}}{2}\sin^{2}\theta\right),$
(86)

where $K$ is the total kinetic energy of the bosons and can be calculated to
be

$K=\frac{1}{2}\sum_{\mathbf{k}}\left(\omega^{+}_{\mathbf{k}}+\omega^{-}_{\mathbf{k}}\right).$
(87)

In Fig. 16,

Figure 16: The energy per unit cell (in units of $V$) of saddle points
$\theta=0$ (QSL) and $\theta=\pi/2$ (ice ferromagnet) given by Eq. (86) up to
an overall additive constant that is the same for $\theta=0$ and
$\theta=\pi/2$. We also plot $-\frac{\Omega^{2}}{4V^{2}}$ arising from trivial
spin-flip pairs: this plot almost overlaps with the energy of the $\theta=0$
state.

we plot the energy $E_{\text{MF}}$ for $\theta=0$ (QSL) and $\theta=\pi/2$
(ice ferromagnet), and find a transition at $\Omega\approx 0.13V$. (The third
solution for $\theta$ becomes the lowest-energy solution only in a minuscule
window around $\Omega\approx 0.13V$, so we ignore it.) However, we will now
argue that this result is misleading.

In gMFT, the energy reduction in the QSL phase with respect to the ordered
phase (ice ferromagnet) arises from the minimization of kinetic energy of the
bosonic charges $\hat{\Phi}_{\mathbf{r}}$ that are allowed to hop. When
$\theta=0$, the hopping coefficient is maximized, while, for $\theta=\pi/2$,
the hopping coefficient is 0. However, microscopically, this hopping
corresponds to a single spin-flip. A pair of spin-flips at the same site leads
to a constant reduction of energy coming from second order perturbation
theory, given by $-\Omega^{2}N_{\text{u.c.}}/V$. It is constant in the sense
that this reduction is obtained for any state including the QSL and the ice
ferromagnet. The mean-field calculation, however, unfairly assigns this
reduction to the QSL but not to the ordered state. In fact, in Fig. 16, we
have also plotted $-\Omega^{2}/(4V)$ (the factor of 1/4 can perhaps be
attributed to using spin-$1/2$ and classical spins at the same time). As can
be seen, this plot almost completely overlaps with the energy of the QSL
calculated within gMFT. So it is clear that, within gMFT, the difference
between the energies of the QSL and the confined phase is quadratic in
$\Omega$ to leading order even though we know from perturbation theory that
the leading order term should be proportional to $\Omega^{6}$. Hence, gMFT
cannot be used in the vicinity of the confinement-deconfinement transition
unless gauge-fluctuations are properly taken into consideration.

##  Appendix B Corrections in the measurement protocol of the plaquette X
correlator

In the protocol described in Sec. IV.1.1, we equated the expectation value of
$\prod_{i=1}^{12}(2\hat{S}_{i}^{x})$ with $\hat{X}_{P}\hat{X}_{P^{\prime}}$ in
the ground state (more generally, the state prepared in experiment) assuming
that the state is supported entirely on the ice manifold. Corrections arise
because this is satisfied only approximately. However, we show in this
appendix that these corrections are of sixth order in $\Omega/V$.

Let $\ket{\Psi_{g}}$ be the ground state of the system. This implies that
$\ket{\Psi_{0}}=\hat{U}_{S}\ket{\Psi_{g}}$ is in the ice manifold, where
$\hat{U}_{S}$ is the unitary operator that implements the Schrieffer-Wolff
transformation (see Sec. III.1.1). First we show that the first order
correction in $\Omega/V$ is zero in the limit $\Omega\ll V$.

The quantity measured in an experiment implementing the protocol of Sec.
IV.1.1, assuming that the experiment prepares the ground state, is

$\displaystyle C_{X}^{\text{exp}}$
$\displaystyle=\langle\Psi_{g}|e^{i\hat{H}_{\text{Y}}t_{\text{Y}}}\prod_{i=1}^{12}(2\hat{S}^{z}_{i})e^{-i\hat{H}_{\text{Y}}t_{\text{Y}}}|\Psi_{g}\rangle$
$\displaystyle\approx\langle\Psi_{0}|\hat{U}_{S}\hat{A}\hat{U}_{S}^{\dagger}|\Psi_{0}\rangle,$

where, $\hat{A}=\prod_{i=1}^{12}2\hat{S}_{i}^{x}$. Now,

$\displaystyle\hat{U}_{S}\hat{A}\hat{U}_{S}^{\dagger}=$
$\displaystyle\hat{A}+\left[\hat{S},\hat{A}\right]+\frac{1}{2!}\left[\hat{S},\left[\hat{S},\hat{A}\right]\right]$
(88)
$\displaystyle+\frac{1}{3!}\left[\hat{S},\left[\hat{S},\left[\hat{S},\hat{A}\right]\right]\right]+\cdots.$

Since $\ket{\Psi_{0}}$ is in the ice manifold, the expectation value of the
first term above is the required plaquette X correlator,
$\langle\Psi_{0}|\hat{A}|\Psi_{0}\rangle=\langle\Psi_{0}|\hat{X}_{P}\hat{X}_{P^{\prime}}|\Psi_{0}\rangle$,
as explained in Sec. IV.1.1. The expectation value of the remaining terms is
the error. $\hat{S}$ can be written as $\sum_{i=1}^{\infty}\hat{S}_{i}$ where
$\hat{S}_{i}$ is of order $(\Omega/V)^{i}$. First, we will show that the first
order error,
$\mathcal{E}_{1}=\langle\Psi_{0}|\left[\hat{S}_{1},\hat{A}\right]|\Psi_{0}\rangle$,
is zero using the expressions for $\hat{S}$ up to first order in $\Omega/V$
from Appendix B of [40]:

$\hat{S}_{1}=\hat{\mathcal{P}}\hat{H}_{\Omega}\hat{\mathcal{{D}}}-\hat{\mathcal{D}}\hat{H}_{\Omega}\hat{\mathcal{P}},\;\text{where
}\hat{\mathcal{D}}=-\frac{\hat{1}-\hat{\mathcal{P}}}{\hat{H}_{0}-E_{0}}.$ (89)

Here, $\hat{\mathcal{P}}$ is the projector into the ice manifold, $E_{0}$ is
the energy of the ice manifold, and $\hat{H}_{0}$ and $\hat{H}_{\Omega}$ are
given in Eq. (6). Plugging in the expression for $\hat{S}_{1}$ in
$\mathcal{E}_{1}$ gives

$\displaystyle\mathcal{E}_{1}$
$\displaystyle=\langle\hat{\mathcal{P}}\hat{H}_{\Omega}\hat{\mathcal{{D}}}\hat{A}-\hat{A}\hat{\mathcal{P}}\hat{H}_{\Omega}\hat{\mathcal{{D}}}-\hat{\mathcal{D}}\hat{H}_{\Omega}\hat{\mathcal{P}}\hat{A}+\hat{A}\hat{\mathcal{D}}\hat{H}_{\Omega}\hat{\mathcal{P}}\rangle_{0}$
(90)
$\displaystyle=\langle\Psi_{0}|\hat{H}_{\Omega}\hat{\mathcal{{D}}}\hat{A}|\Psi_{0}\rangle+\langle\Psi_{0}|\hat{A}\hat{\mathcal{D}}\hat{H}_{\Omega}|\Psi_{0}\rangle.$

We have used $\hat{\mathcal{D}}\ket{\Psi_{0}}=0$ and
$\hat{\mathcal{P}}\ket{\Psi_{0}}=\ket{\Psi_{0}}$ above, and the expectation
values are taken in $\ket{\Psi_{0}}$. Since $\hat{A}$ flips the 12 spins on
$P$ and $P^{\prime}$, while $\hat{H}_{\Omega}$ flips only one spin, it is not
possible to go back to the ice manifold after applying
$\hat{H}_{\Omega}\hat{\mathcal{D}}\hat{A}$ on $\ket{\Psi_{0}}$. Thus
$\langle\hat{H}_{\Omega}\hat{\mathcal{D}}\hat{A}\rangle_{0}=0$, implying
$\mathcal{E}_{1}=0$.

The corrections at order $i$ involves expectation value in $\ket{\Psi_{0}}$ of
operators that consist of $i$ factors of $\hat{H}_{\Omega}$ and one factor of
$\hat{A}$. For example, the second-order error is

$\mathcal{E}_{2}=\expectationvalue{\left[\hat{S}_{2},\hat{A}\right]+\frac{1}{2!}\left[\hat{S}_{1},\left[\hat{S}_{1},\hat{A}\right]\right]}{\Psi_{0}},$
(91)

where

$\displaystyle\hat{S}_{2}=$
$\displaystyle\mathcal{P}\hat{H}_{\Omega}\hat{\mathcal{D}}\hat{H}_{\Omega}\hat{\mathcal{D}}-\hat{\mathcal{D}}\hat{H}_{\Omega}\hat{\mathcal{D}}\hat{H}_{\Omega}\hat{\mathcal{P}}-\hat{\mathcal{P}}\hat{H}_{\Omega}\hat{\mathcal{P}}\hat{H}_{\Omega}\hat{\mathcal{D}}^{2}$
(92)
$\displaystyle+\hat{\mathcal{D}}^{2}\hat{H}_{\Omega}\hat{\mathcal{P}}\hat{H}_{\Omega}\hat{\mathcal{P}}.$

It can be seen that, at order $i=1,2,3,4,$ and $5$, the error will be zero
because $\hat{A}$ would flip 12 spins, while $i$ factors of $\hat{H}_{\Omega}$
would flip only $i$ spins. At sixth order, however, the error can be nonzero
if $\ket{\Psi_{0}}$ has support on a configuration where one of the two
plaquettes, say $P^{\prime}$, is flippable. In such a situation, $\hat{A}$
will map $P^{\prime}$ to the complementary flippable configuration and will
also flip the spins on $P$. The six factors of $\hat{H}_{\Omega}$ can bring
$P$ back to the original configuration, giving a state within the ice manifold
which could have nonzero overlap with $\ket{\Psi_{0}}$.

In conclusion, we find that the experimentally measured correlator,
$C_{X}^{\text{exp}}$, is the same as the theoretically needed correlator,
$\langle\Psi_{0}|\hat{X}_{P}\hat{X}_{P^{\prime}}|\Psi_{0}\rangle$, with
corrections of order $(\Omega/V)^{6}$:

$C_{X}^{\text{exp}}=\langle\Psi_{0}|\hat{X}_{P}\hat{X}_{P^{\prime}}|\Psi_{0}\rangle+\Theta\left((\Omega/V)^{6}\right).$
(93)