# Self-Dual Criticality in Three-Dimensional $Z_2$ Gauge Theory with Matter

Andrés Somoza, Pablo Serna, Adam Nahum

► To cite this version:

Andrés Somoza, Pablo Serna, Adam Nahum. Self-Dual Criticality in Three-Dimensional $Z_2$ Gauge Theory with Matter. Physical Review X, 2021, 11 (4), pp.041008. 10.1103/PhysRevX.11.041008. hal-03456267

HAL Id: hal-03456267

https://hal.science/hal-03456267v1

Submitted on 11 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Self-Dual Criticality in Three-Dimensional $Z_2$ Gauge Theory with Matter

Andrés M. Somoza,¹ Pablo Serna,¹,² and Adam Nahum²,³

¹Departamento de Física—CIOyN, Universidad de Murcia, Murcia 30.071, Spain

²Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, 75005 Paris, France

³Theoretical Physics, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom

<figure><img src="image_2.png" /></figure>

(Received 28 March 2021; revised 14 June 2021; accepted 20 July 2021; published 12 October 2021)

The simplest topologically ordered phase in 2 + 1D is the deconfined phase of $Z_2$ gauge theory (realized in the toric code, for example). This phase permits a duality that exchanges electric and magnetic excitations ("e" and "m" particles). The phase transition where one of these particles condenses, while the other remains gapped, has 3D Ising exponents. But the transition out of the deconfined phase when self-duality symmetry is *preserved* is more mysterious. It has so far been unclear whether this transition is continuous, but if continuous, it may be the simplest critical point for which a useful continuum Lagrangian is still lacking. These questions are relevant to soft matter, too, since the gauge theory also describes classical membranes in 3D. Here, we study the self-dual transition with Monte Carlo simulations of the $Z_2$ gauge-Higgs model on cubic lattices of linear size $L \le 96$. Our results indicate a continuous transition, for example via a striking parameter-free scaling collapse. We use duality symmetry to distinguish the leading duality-odd and duality-even scaling operators $A$ and $S$. We explain why standard techniques for locating the critical point are ineffective, and we develop an alternative using "renormalization group trajectories" of cumulants. We check that two- and three-point functions are scale invariant, with scaling dimensions $x_A$ and $x_S$ (autocorrelations in the Monte Carlo dynamics also yield a dynamical exponent $z$). Separately, we propose a general picture for emergent 1-form symmetries, in terms of "patching" of membranes or world surfaces. We relate this to the percolation of anyon worldlines in spacetime. The latter yields a fourth exponent for the self-dual transition. We propose variations of the model for further investigation.

DOI: 10.1103/PhysRevX.11.041008

Subject Areas: Particles and Fields, Statistical Physics
Strongly Correlated Materials

## I. INTRODUCTION

Continuum field theory provides a language for a huge range of classical and quantum phase transitions [1,2], including many cases for which a simple Landau-Ginsburg formulation is insufficient [3–10]. For example, a wide range of topological phase transitions, lacking any local order parameter [3], may be brought under some measure of analytical control using the language of continuum gauge theory, together with various kinds of perturbative expansion (e expansions, large N expansions, etc.). However, despite the wild success of the field theory approach to critical phenomena, there exist phase transitions in simple and natural models that still remain out of reach of field theory tools. This paper characterizes what we suggest is the paradigmatic example of these mysterious transitions, which is the “self-dual” phase transition

between confined and deconfined phases of $Z_2$ gauge theory in three dimensions [11–14].

The $Z_2$-gauge-Higgs model [3,11,15] has a stable deconfined phase, as well as a trivial phase, in three dimensions. In the context of quantum systems in 2 + 1D (the model also has a 3D classical interpretation that we discuss below), the deconfined phase is the simplest $Z_2$ spin liquid [16–21]: the phase of matter realized, for example, by the toric code [22]. The anyon excitations of this phase include quasiparticles denoted $e$ and $m$, with nontrivial mutual statistics, which correspond to the charge and flux in the gauge theory.

The simplest lattice field theory formulation of the $Z_2$-gauge-Higgs model has a two-dimensional parameter space [3,11–13,15,23]. In the quantum language, these two couplings allow us to separately control the masses of both $e$ and $m$ excitations. In the language of the lattice gauge theory, one of the couplings controls the “stiffness” associated with fluctuations of the matter field, and the other, the stiffness of the gauge field: see the schematic in Fig. 1.

While there are only two stable phases in Fig. 1, there are various possibilities for the transition between

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Topology of the phase diagram of $Z_2$ gauge theory with matter. The shaded region is the trivial phase. The double line represents a first-order line, ending at a standard critical endpoint. The Higgs and confinement transitions have Ising exponents. The question mark shows the region studied in this paper. We give evidence for a scale-invariant, self-dual critical point.</figcaption></figure>

them [6,11–14]. The Higgs and confinement transitions correspond to condensation of the *e* particle and of the *m* particle, respectively. These two lines of transitions are, in fact, completely equivalent, as they are mapped into each other by the crucial duality transformation, which exchanges the two kinds of particles. They are subtle phase transitions with no local order parameter [3]. Nevertheless, they *are* amenable to field theory tools. For example, gauge fluctuations are, in fact, irrelevant at the Higgs transition [11]. Its universal scaling is therefore the same as in the limit of infinite gauge stiffness, where the partition function is simply related to that of the standard Ising model (with a sum over boundary conditions): In a sense, we can define a “fictitious” Ising order parameter. In the language of anyons, the reason that this transition (where *e* condenses) is relatively conventional is because *m* remains massive, which ensures that nontrivial braiding processes are not important at low energies.

By contrast, the nature of the transition out of the deconfined phase on the self-dual line, where there is a symmetry between *e* and *m*, has not been understood. Previous Monte Carlo [13] and series expansion [14] studies gave some evidence for a multicritical point here, but the order of the transition, and the structure of the phase diagram close to this “corner” of the deconfined phase, have not been definitively resolved [13,24]. The argument above, which relates the Higgs transition to a simple Landau theory, no longer applies, so we really have to confront the issue of coarse-graining a discrete gauge field, whose low-lying excitations have nontrivial mutual statistics.

The basic challenge can also be understood in geometrical terms. The gauge-Higgs model describes various phases of fluctuating membranes in three spatial

dimensions [6]. This system is fascinating in its own right, relevant to experiments on amphiphilic membranes, where the deconfined phase is known as the “symmetric sponge” phase [25–31]. We argue below that, if the “holes” in these membranes are “small,” and disappear under coarse-graining, the membranes are effectively *closed* surfaces. Mapping them to Ising domain walls is then one way to think about the fictitious Ising order parameter described above. But, as the self-dual point is approached, the holes become large (we demonstrate this explicitly), so this way of thinking breaks down.

In this paper, we determine many of the properties of the self-dual transition using extensive Monte Carlo simulations and arguments based on the renormalization group (RG) and symmetry. Our numerical results include the first demonstration of scale invariance at this transition, via scaling collapse of numerous observables. Our results for exponents also raise intriguing theoretical questions about how to understand this transition.

First, we give strong evidence that the transition is governed by a scale-invariant fixed point, for example, via a striking scaling collapse that does not require any fitting parameters. We classify the leading local operators $S(r)$ and $A(r)$ as even and odd under duality symmetry, respectively, and estimate their scaling dimensions $x_S$ and $x_A$ using scaling collapses and two-point functions. We check that three-point functions are compatible with conformal invariance.

We also address some fundamental aspects of the anyon condensation transitions away from the self-dual line. As noted above, a key feature of the Higgs and confinement transitions is the possibility of using a Landau theory for a “fictitious” Ising order parameter. (These are sometimes referred to as “Ising*” transitions [32].) The emergence of this order parameter may be related to the question of where in the phase diagram certain “1-form” symmetries [33–36] emerge under coarse-graining.

We propose an explicit construction of the fictitious Ising field (and of the string operators of the 1-form symmetry). This construction is based on “repairing” or “patching” the membranes that appear in a geometrical representation of the partition function. We relate the feasibility of this patching operation to the question of whether *e* and *m* worldlines “percolate” in spacetime, and we obtain the phase diagram for this percolation [6] numerically. This approach shows that the fictitious Ising fields can be constructed on the Ising* transition lines, but not at the self-dual critical point. However, we find that scale invariance at the self-dual transition can be diagnosed via the percolation of worldlines, and we compute their universal fractal dimension $d_f$. (The result hints at a possible relation between exponents.)

We discuss the role of self-duality symmetry, arguing that it becomes an emergent internal symmetry in the IR. While our numerical analysis here is for the standard
gauge-Higgs model, we also propose a modified lattice model, with a simpler action of duality, which would be interesting to study further.

The dynamics of the Monte Carlo algorithm (in Monte Carlo time) correspond to a physically sensible universality class for the stochastic dynamics of membranes in 3D. We find that the dynamical exponent for this universality class is $z \approx 2.48$ (not to be confused with the dynamical exponent $z_{QM} = 1$ of the zero-temperature quantum dynamics in the $2+1$D interpretation), and we show that two-time correlation functions are another way to obtain $x_{A,S}$. The fact that the dynamical exponent is large is one of the challenges in simulating this model: Unlike in many simple ordering transitions [37], no efficient nonlocal Monte Carlo update, which reduces $z$ to a small value, is known for this problem.

Various features of the fixed point make standard approaches to determining the precise location of the phase transition point, and the order of the transition, ineffective. These features include the structure of Binder cumulants close to the transition, the lack of any local operator with a small scaling dimension, and the fact that $x_S$ is very close to 1.5. (This is the threshold separating the divergence and convergence of the heat capacity, and the proximity to this threshold leads to large scaling corrections in this quantity.) These features were a key challenge for our initial attempts at data analysis. We describe how they may be overcome, for example by focusing on appropriate dimensionless observables that allow a parameter-free scaling collapse.

Our numerical estimates for the exponents $x_S$ and $x_A$ turn out to be close to certain exponents in the *XY* universality class. This is remarkable in view of the mutual statistics of the condensing quasiparticles [13,14,38–40], which mean that we do not expect a relationship with the *XY* fixed point (Sec. X D). The fixed point studied here is certainly distinct from the *XY* fixed point, as implied, for example, by the very different universal properties of the adjacent phases. On the other hand, it is not hard to find examples of pairs of 3D fixed points with exponents that are fairly close but distinct. This issue is discussed further in Sec. X. There are also many variations of the present model that remain to be studied (Secs. X B, X C, and XI). Figure 2 is a schematic of some of the topics discussed.

Textbook discussions of critical phenomena sometimes give the impression that studying universality in phase transitions is synonymous with studying Lagrangian quantum field theory. Therefore, it is important to remember that there are critical points for which we so far lack any useful continuum Lagrangian (Sec. XI). Given that the self-dual transition is second order, as previously suspected [13,14] and as the numerical evidence presented here shows, then it is perhaps the simplest example of one of these untamed beasts.

However, a rich variety of other topological transitions, with distinct (nontrivially braiding) anyons simultaneously

<figure><img src="image_4.png" /><figcaption>FIG. 2. Some of the topics discussed.</figcaption></figure>

becoming massless, are possible, with other discrete gauge theories providing further simple examples. Models with U(1) symmetry are also interesting [38,41–43], though they are more closely connected to continuum U(1) gauge theory (perhaps with Chern-Simons terms). A systematic program to understand all of these transitions would be valuable. Past results on the formulation of field theories for deconfined phase transitions [10,44–55], where mutually nonlocal fields and Berry phases connecting different gapless degrees of freedom also play a key role, may provide some tools.

The $\mathbb{Z}_2$ deconfined phase is adjacent to another family of critical “quantum loop models” [56–60] with no known Lagrangian description [60]. Interestingly, while these critical points may again be viewed in terms of membranes in spacetime, the obstacle to a continuum description is different there: a topological constraint on the dynamics, rather than the existence of massless particles with nontrivial braiding.

These different kinds of examples suggest that statistical ensembles of membranes [61] in three and four dimensions (elementary “string field theories” [62–64]) still hold many lessons for critical phenomena.

## II. ISING GAUGE MODEL

The gauge-Higgs model has many guises. We begin by reviewing several equivalent formulations of the partition
function we study and the basic features of the phase diagram. Readers should skip topics with which they are familiar.

### A. As a lattice gauge theory

This theory is the standard formulation of $Z_2$ gauge theory, with matter, on a cubic lattice [3,15]. ($Z_2$ gauge theory is also referred to as “Ising” gauge theory.) The degrees of freedom are classical Ising matter fields $\tau_i = \pm 1$ on the sites $i$ of the lattice and gauge fields $\sigma_{ij} = \pm 1$ on the links $\langle ij \rangle$. The action includes a stiffness $K$ for the gauge fluctuations and a coupling $J$ for the matter field. If $\square$ denotes a square plaquette, the partition function is

$$Z \propto \sum_{\{\sigma\}, \{\tau\}} \exp \left( K \sum_{\square} \left( \prod_{\langle ij \rangle \in \square} \sigma_{ij} \right) + J \sum_{\langle ij \rangle} \tau_i \sigma_{ij} \tau_j \right). \quad (1)$$

Throughout, we work on an $L \times L \times L$ torus [65]. The action is invariant under the $Z_2$ gauge transformation $\tau_i \to \tau_i \chi_i$, $\sigma_{ij} \to \chi_i \sigma_{ij} \chi_j$ with $\chi_i = \pm 1$. If desired, we can choose the gauge $\tau_i = 1$, leaving a lattice model for the $\sigma$ spins on the links only, with terms $J\sigma$ on the links, which emphasizes that Eq. (1) is a lattice model with no internal global symmetries [3]. However, along a line in the phase diagram, it has a self-duality symmetry, as discussed below. In parts of the phase diagram, the model also has 1-form symmetries, either microscopic or emergent, which we discuss in Sec. IX.

It will be convenient to define [3,23]

$$x = \tanh K, \quad y = \tanh J. \quad (2)$$

The phase diagram in this parametrization is shown in the main panel of Fig. 3. The dashed line is where the model is self-dual. The approximately rectangular region in the bottom-right corner, at large gauge stiffness $K$ and small matter field coupling $J$, is the deconfined phase supporting deconfined anyons (Sec. II D).

### B. As a model of membranes

The model can be mapped to a statistical ensemble of “membranes” on the cubic lattice [6,66]. In this picture, the parameters $x$ and $y$ control the microscopic surface tension for the membrane, and the microscopic line tension for the membrane boundary, respectively. The partition function is (see the Appendix A for details)

$$Z(x, y) = \sum_{\mathcal{M}} x^{|\mathcal{M}|} y^{|\partial \mathcal{M}|}. \quad (3)$$

Here, a membrane configuration $\mathcal{M}$ is simply any set of plaquettes of the cubic lattice: We call the plaquettes in $\mathcal{M}$ “occupied.” The energy of a configuration depends on the total number $|\mathcal{M}|$ of occupied plaquettes in the

<figure><img src="image_5.png" /><figcaption>FIG. 3. Sketch of the phase diagram in the $(x, y)$ plane. The gauge stiffness $K = \tanh^{-1} x$ increases to the right, and the matter field coupling $J = \tanh^{-1} y$ increases upwards (self-dual line, shown as a dashed line). Inset: same phase diagram in the $(y, y')$ coordinates, where duality acts as a reflection. (The exponents below imply that the $e$ and $m$ condensation lines are asymptotically parallel as they approach the self-dual critical point, though the curvature is not visible at this scale.)</figcaption></figure>

configuration and on the total length of the membrane “boundary,” $|\partial \mathcal{M}|$. This length is the number of links where an odd number of occupied plaquettes meet. We refer to these as occupied links.

Note that the deconfined phase occurs in the regime where the membrane surface is cheap but the membrane boundary is expensive. The extreme limit of the deconfined phase is $x = 1$, $y = 0$, where we have an ensemble of closed membranes with vanishing surface tension. We may exit the deconfined phase either by suppressing the membrane area (decreasing $x$ sufficiently) or by tearing holes in the membranes (increasing $y$ sufficiently) [6].

Figure 4 shows part of a membrane configuration taken from a simulation close to the self-dual critical point that we study. Plaquettes in $\mathcal{M}$ have been colored in (arbitrarily), and the boundary links in $\partial \mathcal{M}$ have been marked in red. The membrane representation suggests investigating “geometrical” (percolation-like) observables close to the critical point, as well as thermodynamic ones [6]. We discuss this in Sec. IX, showing that the loops in $\partial \mathcal{M}$ form a scale-invariant ensemble at the self-dual critical point.

The membrane picture is one way to see the duality property of the model [3]. In Eq. (3), the partition function is expressed as a sum over membrane configurations on the original cubic lattice. An alternative graphical representation yields an ensemble of precisely the same form, but for membranes on the dual cubic lattice (Appendix A), with dual values of the plaquette and link fugacities:

$$x' \equiv \frac{1 - y}{1 + y}, \quad y' \equiv \frac{1 - x}{1 + x}. \quad (4)$$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Membranes with boundaries: part of a configuration, close to the critical point on the self-dual line. Occupied plaquettes are shown in color (colors have no meaning). Occupied links, where an odd number of plaquettes meet, are shown in thick red lines.</figcaption></figure>

This pair of mappings shows that $Z(x, y)$ is equal, up to a trivial constant, to $Z(x', y')$ [67]. Below, we show that duality can also be thought of as a conventional symmetry operation. This symmetry is not manifest in the formulations above but may be made apparent in an alternate representation of the partition sum in terms of worldlines of $e$ and $m$ particles (Sec. II C).

Note that, in view of Eq. (4), we are free to choose $(y, y')$ as coordinates for the phase diagram, as in the inset to Fig. 3. The line $y = y'$ is then the self-dual line, where the Boltzmann weights are invariant under duality symmetry.

### C. Manifestly self-dual loop representation

The quantum interpretation reviewed in the next subsection motivates yet another representation of the path integral, in terms of two species of “loops,” which represent worldlines of both $e$ and $m$ particles. This representation makes self-duality manifest. The price is minus signs in the Boltzmann weight, which encode the mutual semion statistics of the anyons $e$ and $m$.

The partition function can be written as (see Appendix A 3 for details)

$$Z \propto 4 \sum_{C_e, C_m} (-1)^{X(C_e, C_m)} y^{|C_e|} y'^{|C_m|}. \quad (5)$$

The electric and magnetic worldline configurations $C_e$ and $C_m$, which we refer to as loop configurations, are sets of “occupied” links on the original and dual lattices, respectively. See Fig. 5. Any even number of occupied links may be adjacent to each site, so the term “loops” is used loosely (see endnote [68] for details). Note that $C_e$ may be identified with

<figure><img src="image_7.png" /><figcaption>FIG. 5. Cubes of the original and dual lattices, with $e$ worldlines (red) and $m$ worldlines (blue) on the former and the latter, respectively. (This configuration has linking number $X = 1$.)</figcaption></figure>

the membrane boundaries $\partial M$ in the previous representation.

The crucial feature in Eq. (5) is the topological factor $(-1)^X$, which gives a factor of $-1$ for each linking between an $e$ worldline and an $m$ worldline. [That is, $X(C_e, C_m) = 0$, $1$ is the $\mathbb{Z}_2$ linking number of the two worldline configurations. It can be computed, for example, by introducing an arbitrary membrane configuration $\tilde{M}$ such that $\partial \tilde{M} = C_e$ and counting the number of intersections between $\tilde{M}$ and $C_m$ modulo 2.] The values of the dual-link fugacities $y$ and $y'$ are defined in Eq. (2).

An interesting model of $U(1)$ (oriented) flux lines with a linking sign has been studied in Ref. [38] (see also variations in Refs. [38,41–43]). That model has many features in common with Eq. (5) and also describes a problem of condensation of anyons with mutual statistics. However, it also has significant differences as a result of a global $U(1) \times U(1)$ symmetry. We expect the $U(1) \times U(1)$ model to be described, at least in principle, by a continuum Chern-Simons gauge theory. The “$\mathbb{Z}_2 \times \mathbb{Z}_2$” loop model in Eq. (5) is also closely related to a quantum wave function in $3+1D$ that sustains $2+1D \mathbb{Z}_2$ topological order on its boundary [69–71].

Returning to Eq. (5), it is possible to sum over $C_m$ exactly (Appendix A 3). We then return to the membrane expression (3) for the partition function, with $C_e = \partial M$. Similarly, integrating out $C_e$ gives the dual membrane picture on the dual lattice. (Note that the line tension $y'$ or $y$ of the species that is integrated out determines the surface tension of the membranes.)

However, the representation (5) makes the duality symmetry that exists when $y = y'$ (i.e., when $x = x'$) manifest. We can think of the symmetry operation as a translation by $(1/2, 1/2, 1/2)$, where the lattice spacing of each cubic lattice is unity. This translation exchanges the cubic lattice with its dual, so it exchanges $e$ and $m$ worldlines. Microscopically, this is not an internal symmetry (since it involves translation), but we will argue in Sec. III that self-duality becomes an internal $\mathbb{Z}_2$ symmetry of the IR theory.

Section X A presents an alternative loop model in which the $e$ and $m$ loops share the same lattice.
## D. Anyons and the toric code in a field

The 3D gauge theory is expected to capture the universal physics of a wide range of 2 + 1D quantum models with a $Z_2$ spin liquid phase. An anisotropic limit of the 3D theory, where the z direction becomes a continuous imaginary time coordinate, maps exactly to the partition function for such a Hamiltonian on the square lattice. See Refs. [9,15] for details of such mappings. Here, we review the basic excitations of the deconfined phase, as well as how they relate to the geometrical pictures above, in qualitative terms.

It is convenient to start with the toric code [22], a particularly simple model lying in the deconfined phase. The degrees of freedom are spin-1/2 s on the links of the square lattice. The Hamiltonian includes a plaquette term and a vertex term:

$$H = -V \left( \sum_{\square} XXXX + \sum_{+} ZZZZ \right). \quad (6)$$

The first product is shorthand for the Pauli-X operators on the four links making up a given plaquette, and the second is for the four Pauli-Z operators on the links touching a given site. Here, V is a coupling constant, which we have chosen to be equal for the two terms to ensure self-duality symmetry. Noting that we can equally well view spins either as living at the midpoints of bonds on the original square lattice or at the midpoints of bonds on the dual square lattice, the duality symmetry operation may be viewed as a translation by (1/2, 1/2), together with an exchange of X and Z [72].

In ground states [74], all the plaquette and vertex terms in Eq. (6) are equal to 1, and they are superpositions of “strings”—either on the original lattice (if we use the Z basis, a link with $Z = -1$ is regarded as part of a loop) or the dual lattice (if we use the X basis) [22]. There are two fundamental types of excitation, related by duality [75]. A vertex where $ZZZZ = -1$ is an e particle, and a plaquette where $XXXX = -1$ is an m particle. These are distinct types of anyons. Each is a boson, but adiabatically braiding an e around an m changes the sign of the wave function. (In other words, they are mutual semions. The combination of an e particle and an m particle forms another type of anyon whose topological sector is denoted “e”: This also has -1 statistics with e and m, but it is a fermion [22,76].) In the Z basis, an e excitation is the endpoint of a string (Fig. 6).

The toric code is a fine-tuned limit in which the e and m particles are nondynamical. Critical phenomena are possible when the model is perturbed so that pair creation and annihilation of these particles becomes possible. This result may be achieved by adding magnetic fields in both the X and Z directions [13,14]. (For example, adding the operator X to the Hamiltonian allows both hopping and pair creation or annihilation of bare e particles on a given link.) The resulting model has been intensely studied [13,14,24,77] (including as an effective theory for an anisotropic

<figure><img src="image_8.png" /><figcaption>FIG. 6. Relation between e particles and membranes. Left diagram: In the Z basis, the toric-code wave function is a superposition of strings of occupied links (bold) representing $Z = -1$, and an e excitation is a vertex where an odd number of occupied links meet. Right diagram: Constructing a path integral (for a generic perturbed model) in this basis, world surfaces of strings become membranes $M$ (grey), and worldlines of e particles form the membrane boundary, or “loops,” $\partial M = C_e$. (The figure on the right is schematic: In the model we study, membrane configurations look like Fig. 4.)</figcaption></figure>

Kitaev-Heisenberg-Γ magnet [40]). Duality exchanges the two magnetic fields, so the line $h_X = h_Z$ preserves duality symmetry.

The phase diagram of the toric code in X and Z fields is expected to be equivalent to that discussed in the previous sections, up to nonuniversal constants [13]. The dimensionless field $h_X/V$, which can induce condensation of the e particle, plays the role of the vertical coordinate in Fig. 1, and $h_Z/V$ plays the role of the horizontal coordinate.

The connection with the geometrical pictures above arises from writing the imaginary-time partition function in various choices of basis. We describe this only in qualitative terms below.

In the Z basis, the wave function is a superposition of terms like that illustrated in Fig. 6 (left) with e particles (at sites of the square lattice) forming endpoints of strings. Constructing the sum over Feynman trajectories using this basis, the world surfaces of strings form a set of membranes $M$, and the worldlines of e particles form a set of loops that are the boundaries $C_e = \partial M$ of these membranes. (In the limit $h_X = 0$, there are no bare e particles in the ground state, and correspondingly, the membranes are closed surfaces.) This picture is a continuous-time version of that in Sec. II B. The dual membrane picture is obtained by working in the X basis, where the worldlines of m particles are manifest.

Alternately, we may pick a basis in which both the plaquette products XXXX and the vertex products ZZZZ are diagonal, which is possible since all of these terms commute. The Feynman trajectory sum is then over world-line configurations $C_e$ and $C_m$ for both e and m particles, and it is a continuous-time version of the loop model in Eq. (5).

## E. Ising* and first-order lines

To conclude this overview of the model, we recap some features of the phase transition lines in Fig. 1 or, equivalently, Fig. 3.
Starting in the deconfined phase (in the lower right-hand corner of Fig. 3), we may exit it in three ways, two of which are related by duality [3,11,13]. Condensing the *e* particle while keeping the mass of *m* finite corresponds to the upper boundary of the deconfined region in Fig. 3 (main panel); this is the Higgs transition in the lattice gauge theory. Condensing *m* while *e* remains massive is equivalent by duality, and it is the left-hand boundary of the deconfined region in Fig. 3 (main panel). This case is the confinement transition in lattice gauge theory.

These transitions are continuous with Ising exponents, at least sufficiently close to the boundaries of the phase diagram [15], as we now rapidly review. These transitions, described by a weakly gauged Landau theory, are sometimes referred to as Ising* transitions (see, for example, Ref. [32]) to denote the fact that, because of gauging, only the $Z_2$-even operators of the Ising CFT survive as local operators.

Consider, first, the Higgs transition in the limit $x = 1$, i.e., infinite gauge stiffness $K = \infty$. The freezing of gauge fluctuations in this limit gives an exact mapping to a standard cubic lattice Ising model. But Ising exponents are retained at least along some part of the phase transition line, for finite $K$, which can be argued by fixing the gauge and deriving an effective longer-range Ising Hamiltonian perturbatively in $e^{-K}$ [3,11]. In the quantum language, the point is that, as long as the $m$ and $e$ anyons are gapped, we can, for many purposes, neglect the fact that the *e* particle that is condensing is an anyon, rather than a local excitation [39]. By duality, equivalent points hold for the confinement transition.

A more intuitive way to understand the relation to Ising is developed in Sec. IX and Ref. [78]. Let us start with the membrane representation in the limit where the membrane boundary is completely suppressed [$y = 0$ in Eq. (3), or $y' = 0$ in the dual membrane picture]. Again, this corresponds to the transitions on the boundaries of the phase diagram. In this limit, the membranes form closed surfaces, so they may be mapped exactly to domain walls in a nearest-neighbor Ising model (with a sum over boundary conditions [79]).

Now, when we increase $y$ slightly, the membranes acquire holes in them, which means that there is no longer an unambiguous mapping to an Ising model. But if the holes are sufficiently small, we might expect this ambiguity to be unimportant on large scales, so we can again think in terms of an ordering transition for a fictitious Ising order parameter.

We make this idea of a fictitious Ising order parameter precise using an explicit construction, based on the idea of “repairing” or “patching” the membranes in the representation (3) of the partition function (Sec. IX and Ref. [78]). We argue that this construction can be performed all the way along the Ising* critical lines but not at the self-dual critical point, where a different universality class takes over.

For completeness, let us note that we can also think of the Ising* transition in the loop model representation, Eq. (5). Loosely speaking, when one species of loops has a finite typical size, coarse-graining beyond this scale gives a loop model for a single species of unoriented loops. Models of this kind are a standard representation of the Ising universality class, in terms of worldlines of the Ising quanta.

The self-dual transition point [11–14], where the two Ising* lines meet, will be discussed in the rest of the text.

The line of first-order transitions occurs within the trivial phase, so it is relatively conventional. It is also a line where self-duality symmetry is spontaneously broken. The natural expectation is that the critical endpoint of this line is in the Ising universality class, with the Ising order parameter being the anti-self-dual operator defined below. Ising universality for this critical endpoint is consistent with a very rough estimate of the universal crossing value of the Binder cumulant, as shown in Appendix B.

In addition to these thermodynamic transitions, we may also define geometrical transitions [6] using the geometry of the membranes in Eq. (3) (Appendix D).

### III. SELF-DUALITY AS A SYMMETRY

We anticipate that, for any scale-invariant critical point on the self-dual line, self-duality becomes an internal $Z_2$ symmetry of the IR theory.

One way to argue for this is via the manifestly self-dual representation of the partition function in Eq. (5) with $y = y'$, which has a translation symmetry by $(1/2, 1/2, 1/2)$, exchanging *e* and *m* worldlines. Correspondingly, the 2D quantum model in Sec. II D has a symmetry involving translation by $(1/2, 1/2)$, which exchanges *e* and *m* particles.

The simplest assumption is that, at a scale-invariant critical point, this microscopic symmetry gives rise to an internal $Z_2$ symmetry of the IR fixed point theory. Loosely speaking, the action of the translation on the rescaled spatial coordinate of the coarse-grained theory disappears in the IR limit, so the microscopic symmetry transformation should map to a purely internal symmetry transformation on the operators of the IR theory [80].

We can motivate this further by noting that alternative models for the deconfined phase can be constructed in which the duality symmetry, exchanging *e* and *m*, is an internal symmetry even at the lattice level. References [82,83] give exactly solvable 2D string-net Hamiltonians for the deconfined phase with this property [84]. We may also define a variant of the 3D loop model (5) in which the *e* and *m* loops live on the same lattice, with a $Z_2$ symmetry exchanging them. This model is defined in Sec. X A. It is plausible that by varying the interactions in either of these models we could access the same self-dual fixed point, at the corner of the deconfined phase, as in the original model.
<figure><img src="image_9.png" /><figcaption>FIG. 7. Phase diagram on the self-dual line, i.e., on the line $y = (1-x)/(1+x)$ (where $x = x'$ and $y = y'$).</figcaption></figure>

The phase diagram of the gauge-Higgs model, restricted to the self-dual line, is shown schematically in Fig. 7. The first-order line in Fig. 1 corresponds to spontaneous breaking of self-duality symmetry, as discussed below.

### A. Defining (anti)symmetric operators

Duality acts on the phase diagram as $y \leftrightarrow y'$. We now define lattice operators $S$ and $A$, which are conjugate to the self-dual and anti-self-dual couplings, namely, $y + y'$ and $y - y'$, respectively.

We continue to use the language of membranes (Sec. II B). First, define “face” and “edge” operators, $\mathcal{F}(p)$ and $\mathcal{E}(\ell)$, respectively, which are equal to either 0 or 1 and which measure whether a given plaquette $p$ or link $\ell$ of the cubic lattice is occupied in the membrane configuration $M$. In other words, $\mathcal{F}(p) = 1$ if $p \in M$ and $\mathcal{F}(p) = 0$ if $p \notin M$; $\mathcal{E}(\ell) = 1$ if $\ell \in \partial M$ and $\mathcal{E}(\ell) = 0$ if $\ell \notin \partial M$. (See Appendix A for expressions for $\mathcal{F}$ and $\mathcal{E}$ in terms of $\sigma$ and $\tau$.)

The duality transformation maps these operators to operators on the dual lattice which play an equivalent role in the dual membrane ensemble. From the expressions in Appendix A, or, alternately, by extending the duality transformation to the case of spatially varying couplings $x$ and $y$, we see that this mapping is

$$ \mathcal{F} \rightarrow -\frac{2x}{1-x^2} \mathcal{E} + \frac{x}{1+x}, \quad (7) $$

$$ \mathcal{E} \rightarrow -\frac{2y}{1-y^2} \mathcal{F} + \frac{y}{1+y}. \quad (8) $$

We have suppressed plaquette and link indices to avoid clutter. The transformed operator on the rhs is located at the link or plaquette that is dual to the plaquette or link of the operator on the rhs.

Next, let us symmetrize these operators with respect to the lattice point group, which naturally leads to operators that are centered either on a cube of the lattice or on a vertex. We use $\mathcal{F}_{\text{cube}}(c)$ to denote the sum of $\mathcal{F}$ over the six plaquettes of a cube $c$, and $\mathcal{F}_{\text{vertex}}(v)$ to denote (one-half times) the sum of $\mathcal{F}$ over the 12 plaquettes that touch a vertex $v$. Similarly, $\mathcal{E}_{\text{cube}}(c)$ is (one-half times) the sum over the 12 links in a cube, and $\mathcal{E}_{\text{vertex}}(v)$ is the sum over the six links touching a vertex. (We include the factors of $1/2$ so that the expectation values of $\mathcal{E}_{\text{cube}}$ and $\mathcal{E}_{\text{vertex}}$ are equal, and similarly for $\mathcal{F}_{\text{cube}}$ and $\mathcal{F}_{\text{vertex}}$.)

Finally, specializing to the self-dual line, we define

$$ A_{\text{cube}} = \mathcal{F}_{\text{cube}} + \frac{2x}{1-x^2} \mathcal{E}_{\text{cube}} - \frac{6x}{1+x}, \quad (9) $$

$$ S_{\text{cube}} = \mathcal{F}_{\text{cube}} - \frac{2x}{1-x^2} \mathcal{E}_{\text{cube}} + \frac{6x}{1+x}, \quad (10) $$

and analogously for operators $A_{\text{vertex}}$ and $S_{\text{vertex}}$ at the vertices.

These operators transform simply under duality:

$$ A_{\text{cube}} \leftrightarrow -A_{\text{vertex}}, \quad (11) $$

$$ S_{\text{cube}} \leftrightarrow +S_{\text{vertex}}. \quad (12) $$

In addition,

$$ \sum_c A_{\text{cube}}(c) = \sum_v A_{\text{vertex}}(v), \quad (13) $$

$$ \sum_c S_{\text{cube}}(c) = \sum_v S_{\text{vertex}}(v). \quad (14) $$

These “integrated” operators, which can be written either as sums over cubes or vertices, are the anti-self-dual and self-dual perturbations of the self-dual line.

Now, we expand the lattice operators above in terms of continuum operators of a putative IR fixed point. Denote the leading $\mathbb{Z}_2$-odd and $\mathbb{Z}_2$-even scalar continuum operators by $A(r)$ and $S(r)$, respectively, with no subscript. We also write $A_{\text{cube}}(r)$, $A_{\text{vertex}}(r)$, etc., for lattice operators, where $r$ is the location of the appropriate cube or vertex.

To be consistent with Eqs. (11) and (14), the operator $S_{\text{cube}}$ must be of the form

$$ S_{\text{cube}}(r) = (\text{self-dual operators}) + (\text{derivatives of anti-self-dual operators}), $$

and analogously for the other lattice $A$ and $S$ operators, in order that their integrated versions have well-defined symmetry under duality. Taking into account point-group symmetry, some of the allowed terms in $S_{\text{cube}}$ and $S_{\text{vertex}}$ are

$$ S_{\text{cube}}(r) = \alpha S(r) + \beta \nabla^2 A(r) + \gamma \nabla^2 S(r) + \dots, \quad (15) $$

$$ S_{\text{vertex}}(r) = \alpha S(r) - \beta \nabla^2 A(r) + \gamma \nabla^2 S(r) + \dots, \quad (16) $$

Here, $\alpha, \beta$, and $\gamma$ are nonuniversal constants. The sign of the $\mathbb{Z}_2$-odd term is reversed in the second line, in order that mixed correlators of lattice operators are consistent with Eq. (11). Equivalent formulas apply for the lattice $A$ operators, with $A$ and $S$ exchanged, and separate nonuniversal constants.

We use the operators $A_{\text{cube}}$ and $S_{\text{cube}}$ in our simulations. We see that these lattice operators may be identified (up to
derivative operators and other operators that are expected to be highly irrelevant) with the leading self-dual and anti-self-dual continuum operators.

From now on, we denote the lattice operators simply as $A(r)$ and $S(r)$, where $r$ is the coordinate of a cube. We use $A$ or $S$, without an argument, to denote the spatially averaged quantity—for example,

$$A = \frac{1}{L^3} \sum_r A(r). \qquad (17)$$

We write $x_A$ and $x_S$ for the scaling dimensions of the two operators.

## B. Spontaneous breaking of duality symmetry

The phase diagram on the self-dual line is shown in Fig. 7. Importantly, $A$ in Eq. (17) is an order parameter for the symmetry breaking that occurs when we exit the deconfined phase. By self-duality symmetry, its average vanishes,

$$\langle A \rangle = 0, \qquad (18)$$

but in the duality-broken phase, its magnitude $\sqrt{\langle A^2 \rangle}$ remains nonzero in the thermodynamic limit.

Raw data for this quantity are shown in Fig. 8, close to the critical point of interest. In all plots, we parametrize the position along the self-dual line with $x$, so the deconfined phase corresponds to the right-hand side of the figure. At first glance, Fig. 8 is consistent with the order parameter becoming nonzero in a continuous fashion below some $x_c$ (whose estimation we discuss below).

The operator $S$, whose average is shown in Fig. 9, is analogous to the “energy” operator at a conventional classical transition since it does not break symmetry: For a continuous transition, the correlation length exponent is

<figure><img src="image_10.png" /><figcaption>FIG. 8. Duality-breaking order parameter $\sqrt{\langle A^2 \rangle}$ as a function of $x$ on the self-dual line, for various system sizes (indicated in the legend). The lines are just a guide to the eye. The deconfined phase is at larger $x$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. $\langle S \rangle$ as a function of $x$ for different system sizes. The legend indicates system size. The lines are just a guide to the eye.</figcaption></figure>

$$\nu = \frac{1}{3 - x_S}. \qquad (19)$$

The data in Figs. 8 and 9 may be restated in terms of the average occupation of plaquettes and links in the membrane picture. On the section of the self-dual line where duality is spontaneously broken, there are two coexisting equilibria with different plaquette and link densities. We plot these two solutions explicitly in Appendix B. In the critical regime of interest here, the average occupation number of links is relatively small, approximately 2.5%, but despite this feature, they make up a scale-invariant ensemble of loops (Sec. IX).

We now discuss how to establish the universal properties of the transition.

# IV. SCALE INVARIANCE

## A. Initial obstacles

One standard means of locating a phase transition is to analyze the specific heat, which, for many simple ordering transitions, diverges at the critical point. If so, data for different system sizes can typically be scaled, allowing the critical point and correlation length exponent to be determined. Here, the variable $S(r)$ plays the role of an energy, as discussed in the previous section, and $L^3 \text{var}(S)$ is analogous to a specific heat. Values for different system sizes are shown in Fig. 10.

At first sight, the behavior is the expected one: Curves show a peak. But on closer inspection, it is unclear whether the peak diverges at large $L$ or tends to a constant. It also becomes clear that variation of the width and height of the peaks does not follow the simple scaling form

$$\text{Var}(S) = L^{-2x_S} f(z), \qquad (20)$$

where $z = (x - x_c)L^{1/\nu}$ and $\nu = 1/(3 - x_S)$. At first glance, it looks like this transition will be plagued by large finite-size effects, and it will be difficult to see any sign of scale invariance. In fact, this is not the case, which will
<figure><img src="image_12.png" /><figcaption>FIG. 10. Heat capacity: variance of S (multiplied by L³) as a function of x for different system sizes. The legend indicates system size. The lines are B-spline fits and are just a guide to the eye.</figcaption></figure>

become clear after analyzing the behavior of the variable A. (We will return to the specific heat below.)

Another standard tool to determine the location of a critical point is the Binder cumulant for the order parameter [88]. Here, A is our order parameter, and we define a rescaled [89] version of the Binder parameter:

$$b_4(A) = - \frac{\kappa_4(A)}{2\kappa_2(A)^2}, \qquad (21)$$

where $\kappa_n(A)$ is the nth order cumulant. With this normalization, $b_4(A)$ becomes 0 in the deconfined phase (where A is disordered and has a Gaussian distribution) and tends to 1 in the first-order coexistence region (where A has a two-delta distribution).

At a conventional second-order symmetry-breaking transition (e.g., Ising), the Binder parameter varies monotonically from 0 to 1, and different system sizes show a crossing that allows accurate location of the critical point. This is not the case here, as shown in Fig. 11. Rather than crossing, the curves present a minimum near $x = 0.6367$. However, the curves do tend to touch here, consistent with

<figure><img src="image_13.png" /><figcaption>FIG. 11. We show $b_4(A) = -(1/2)\kappa_4(A)/\text{Var}(A)^2$ as a function of x for different system sizes. The legend indicates system size. The lines are B-spline fits and are just a guide to the eye.</figcaption></figure>

scale invariance [$b_4(A)$ is a dimensionless quantity, which should be asymptotically L independent at a critical point], although some finite-size effects can be appreciated. A previous estimate of $x_c \approx 0.6359$ [13] is not consistent with the location of the minimum (we will give a more accurate estimate below).

So, after a first look at these two standard quantities, it is hard to assess whether the data obey scaling collapse, and it seems, at first sight, that accurate estimation of $x_c$ will be more troublesome than for other systems and may be plagued by finite-size effects. Having reached this point, a key step for our understanding was analyzing a parameter-free scaling collapse, which we describe next.

## B. Parameter-free scaling collapse; RG trajectories

We advocate using a parameter-free procedure to determine the quality of scaling collapse of the data near a critical point. We construct a parametric plot using as coordinates two dimensionless quantities, $b_4(A)$ and $b_1(A)$ (defined below).

In the scaling region,

$$b_4(A) = f(z), \qquad (22)$$

where $z = (x - x_c)L^{1/\nu}$. Other dimensionless ratios of cumulants are candidates for the second dimensionless quantity. Binder defined a ratio based on the sixth-order cumulant: $V_L = \kappa_6(A)/(30\kappa_2(A)^3)$ [88]. High-order cumulants are sensitive to the tails of the distribution and can be difficult to estimate accurately. Therefore, we instead advocate using the ratio $\langle |A| \rangle / \langle A^2 \rangle^{1/2}$:

$$b_1(A) = \frac{1}{1 - \sqrt{2/\pi}} \left( \frac{\langle |A| \rangle}{\kappa_2(A)^{1/2}} - \sqrt{\frac{2}{\pi}} \right). \qquad (23)$$

The coefficients have again been chosen so that $b_1(A)$ tends to 0 in the deconfined phase and to 1 in the coexistence region. Note that $b_1(A)$ behaves qualitatively like $b_4(A)$ in Fig. 11, and its expected scaling form is as in Eq. (22), with a different scaling function. For a standard Ising transition, $b_1$ goes monotonically from 0 to 1 with a crossing for different system sizes: There, it can be used to determine the critical temperature, with the advantage of being slightly easier to estimate than $b_4(A)$.

By plotting $b_4(A)$ versus $b_1(A)$, we obtain a parametric plot where $z$ is the parameter; see Fig. 12. If scaling is obeyed, points with different $x$ and $L$, but the same $z$, must overlap. This approach is a fair test of scale invariance because we do not have to fix or fit any parameters by hand and instead just plot raw data.

The data trace a trajectory from the point (0,0) to the point (1,1), showing very good overlap, except near the region $(b_1(A), b_4(A)) \approx (-0.109, -0.285)$ where we see some finite-size effects. However, these finite-size effects
<figure><img src="image_14.png" /><figcaption>FIG. 12. Parameter-free scaling collapse for $b_4(A)$ as a function of $b_1(A)$ for several system sizes (colored). Black stars mark the two phases: the deconfined phase at (0,0) and the broken-duality phase (i.e., first-order coexistence) at (1,1). The inset shows a zoom of the lower-left corner of the main panel. Two selected $x$ values (0.6364 and 0.6370) are highlighted, and arrows are drawn in between consecutive system sizes. Close to the critical point, $x \approx 0.6367$, values remain in a small region around $(b_1(A), b_4(A)) \approx (-0.109, -0.285)$.</figcaption></figure>

become small for $L > 32$. This figure represents, on its own, strong evidence that the multicritical point is a second-order phase transition.

This figure can also be used to estimate the critical point. We construct RG trajectories in the $(b_1(A), b_4(A))$ plane by following the points for a fixed value of $x$ as $L$ is increased. The points representing a system (for a fixed generic $x$) should flow along the universal line towards either the (0,0) or the (1,1) fixed point. The inset of Fig. 12 shows this flow for two different $x$ values. For the system sizes used, this simple procedure already determines the critical point with four digits of precision. We observe that for $x > 0.6367$, the system flows towards (0,0), while for $x < 0.6366$, it flows to (1,1). The repulsive fixed point is located (within the precision of the procedure) at the lower-left extreme of the universal curve.

We see that the data in Fig. 11 should approximately scale for $L \gtrsim 32$ (fits are given below). A similar figure using ratios involving $S$ [e.g., $b_4(S)$ or $\kappa_3(S)/\kappa_2(S)^{3/2}$] does not show good overlap, as expected from the discussion of Fig. 10. It would be strange to have very large finite-size effects in quantities depending on $S$ but not on those depending on $A$. The explanation turns out to be very simple. The exponent $x_S$ is very near 1.5, where the regular contribution to $\text{Var}(S)$ cannot be neglected. When this is taken into account, quantities depending on $S$ also obey scaling (Sec. V B).

## V. CRITICAL EXPONENTS

We turn to scaling fits in order to determine $x_c$ and the scaling dimensions of $A$ and $S$ ($x_A$ and $x_S$, respectively).

<table><caption>TABLE I. Results of fits. Errors shown are purely statistical.</caption><thead><tr><th>Variable</th><th>x<sub>c</sub></th><th>x<sub>S</sub></th><th>x<sub>A</sub></th><th>&chi;<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>b<sub>1</sub>(A)</td><td>0.636660(16)</td><td>1.446(56)</td><td></td><td>49.53</td><td>46</td></tr><tr><td>b<sub>4</sub>(A)</td><td>0.636670(14)</td><td>1.445(62)</td><td></td><td>65.8</td><td>46</td></tr><tr><td>&radic;(&lang;A<sup>2</sup>&rang;)</td><td>0.636702(20)</td><td>1.502(43)</td><td>1.222(16)</td><td>44.9</td><td>40</td></tr><tr><td>(&lang;|A&rang;)</td><td>0.636702(22)</td><td>1.510(48)</td><td>1.221(16)</td><td>43.1</td><td>40</td></tr><tr><td>Var(S)</td><td>0.636661(14)</td><td>1.5(fixed)</td><td></td><td>88.6</td><td>81</td></tr><tr><td>&kappa;<sub>3</sub>(S)</td><td>0.636651(18)</td><td>1.506(9)</td><td></td><td>68.6</td><td>66</td></tr></tbody></table>

Details of how fits were constructed may be found in Appendix C. The results of the various fits are summarized in Table I.

### A. Scaling collapse for A

The scaling form for dimensionless quantities such as $b_1(A)$ involves, in addition to the scaling function, the parameters $x_c$ and $\nu = 1/(3 - x_S)$. Figure 13 shows the scaling collapse of the data for $b_1(A)$ and the fitted scaling function. The critical coupling $x_c$ obtained (Table I) is very near the initial estimation made in Sec. III B, and $x_S$ is near 1.5, as noted above. A fit of $b_4(A)$ gives very similar results (Table I).

In order to obtain the exponent $x_A$, we fit

$$ \sqrt{\langle A^2 \rangle} = L^{-x_A} g(z). \qquad (24) $$

The resulting scaling function $g(z)$ is shown in Fig. 14, and the fitted parameters are indicated in Table I. Fitting $\langle |A| \rangle$ yields similar results.

No finite-size-scaling corrections are included in these fits, although for $L = 32$ these corrections are still non-negligible (compared to the error bars). If data for $L = 32$ are excluded from the fits for $b_1(A)$ and $b_4(A)$, the estimates of $x_S$ increase.

<figure><img src="image_15.png" /><figcaption>FIG. 13. Scaling collapse of $b_1(A)$ versus scaling variable $z = (x - x_c)L^{1/\nu}$, where $1/\nu = 3 - x_S$. The blue line corresponds to the fitted scaling function using B-splines with 12 degrees of freedom. The legend indicates the different system sizes.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 14. Scaling collapse of $L^{x_A} \sqrt{\langle A^2 \rangle}$ versus scaling variable $z = (x - x_c)L^{1/\nu}$, where $1/\nu = 3 - x_S$. The blue line is the fitted scaling function using B-splines with 10 degrees of freedom. The legend indicates system sizes.</figcaption></figure>

## B. Scaling collapse for S

We have suggested above that the failure of a straightforward scaling collapse for $Var(S)$ is due to $x_S$ being very close to 3/2, the threshold where the regular contribution becomes comparable with the scaling contribution. Fortunately, a simple modification of the scaling ansatz should be accurate when $|x_S - 3/2| \ll 1/\log L$:

$$L^3 \mathrm{Var}(S) \approx f(z) + 4\pi\alpha_S^2 \log(L). \quad (25)$$

Here, $\alpha_S^2$ is the normalization constant for the two-point function of S (Sec. VI). The function $f(z)$ includes a z-independent constant contribution, which arises from nonuniversal short-distance correlations [90].

We have performed fits to this form, keeping $\nu = 2/3$ fixed (in line with the approximation above), so only the critical coupling $x_c$ and the coefficient $C = 4\pi\alpha_S^2$ of the logarithm can be adjusted to obtain scaling collapse. The scaling function is shown in Fig. 15. We obtain a good fit, even when data from small system sizes are included. The estimated $x_c$ is again very similar to previous estimates. Also, the constant $C = 10.05(23)$ obtained is consistent with our calculation of the correlation function in the next section. In summary, the fit to $Var(S)$ is consistent with the correlator of S obeying scaling, with $x_S$ very close to 3/2 (indeed, allowing $\nu$ to be free in this fit, instead of fixed to 2/3, did not improve the fit quality).

An alternative way to avoid dealing with the regular contribution is to analyze higher-order cumulants. The singular contribution near a critical point scales as $\kappa_n(S) = L^{-nx_S}f(z)$, while the regular contribution should scale as $\kappa_n(S)_{\text{regular}} \propto L^{-d(n-1)}$. For $x_S \approx 3/2$ and $n = 2$, both contributions scale in the same way, but for $n = 3$, the singular contribution should dominate. Indeed, the data for $\kappa_3(S)$ can be collapsed, as shown in Fig. 16. We obtain $x_c$ and $x_S$ values fully consistent with previous results

<figure><img src="image_17.png" /><figcaption>FIG. 15. Scaling collapse for $L^3(\mathrm{Var}(S) - C \log(L))$ versus the scaling variable $z = (x - x_c)L^{1/\nu}$, where $1/\nu = 3 - x_S$. The blue line corresponds to the fitted scaling function using B-splines with 10 degrees of freedom. The legend indicates different system sizes. In this fit, only $x_S$ has been fixed to 3/2; see text for explanation of the scaling ansatz.</figcaption></figure>

(Table I), although it is worth noting that the statistical error of $x_S$ is much smaller.

## C. Summary of exponents from the fits

We have provided clear evidence that the Ising gauge-Higgs model has a scale-invariant multicritical point. Simulations are inevitably restricted to finite length scales, so they can never rigorously exclude an extremely weak first-order transition; but all of the observables we have examined exhibit good scaling collapse, with fairly modest finite-size effects.

As there are some finite-size effects, we consider a reasonable confidence interval for the critical point to be $x_c \in [0.63665, 0.6367]$. For the study of correlation functions, in the next sections we round to four digits and consider critical behavior at $x_c \approx 0.6367$.

For the exponent $x_S$, the value obtained from $\kappa_3(S)$ (Table I) has the smallest statistical error, and we take it as a

<figure><img src="image_18.png" /><figcaption>FIG. 16. Scaling collapse for $L^{3x_S}\kappa_3(S)$ versus the scaling variable $z = (x - x_c)L^{1/\nu}$, where $1/\nu = 3 - x_S$. The blue line corresponds to the fitted scaling function using B-splines with 10 degrees of freedom.</figcaption></figure>
reference in the following. Our statistical error bars do not take into account possible systematic errors related, for example, to finite-size effects, so a realistic confidence interval should be larger; however, we have verified that dropping the smaller system sizes in the fit only very slightly increases $x_S$, remaining within the statistical error bars.

Our estimate of $x_S$ leads to $\nu = 0.669(4)$. This is not too far from estimates $\nu \approx 0.7$ [14] and $\nu \approx 0.69$ [77] based on the calculation of the gap in the toric code using small-field series expansions (up to eighth order). However, a basic issue with the series expansion method is that it cannot detect a first-order transition [14]; i.e., it must assume a continuous transition rather than demonstrating one. Reference [77] attempts to rectify this by comparing estimates of the ground-state energy from series expansion with a variational wave function, but we expect that the accuracy with which this method could detect a weak first-order transition is severely limited by the accuracy of the variational wave function.

A realistic confidence interval for $x_A$ should again be larger than the statistical one in Table I. We note that the $x_c$ estimates obtained from $\sqrt{\langle A^2 \rangle}$ and $\langle |A| \rangle$ are slightly larger than for the other fits; if we estimate $x_A$ keeping $x_c = 0.63666$ fixed, then the value drops to 1.20, slightly below the statistical confidence interval.

Standard scaling relations [91] imply that the order parameter exponent $\beta$, defined by $\sqrt{\langle A^2 \rangle} \sim (x_c - x)^\beta$ in the infinite system, is $\beta = x_A/(3 - x_S)$ (compare Fig. 8). Asymptotically close to the self-dual critical point, the shape of the Higgs and confinement lines in the inset to Fig. 3 should be $y_- \sim \pm |\delta y_+|^{(3-x_A)/(3-x_S)}$, where $y_\pm = (y \pm y')$ gives the self-dual and anti-self-dual couplings. Since $(3 - x_A)/(3 - x_S)$ is a little larger than one, the Higgs and confinement lines are asymptotically parallel as they approach the critical point.

The values obtained for the critical exponents clearly differ from Ising values, but they are surprisingly close to certain exponents in the XY model. This point will be discussed in Sec. X.

## VI. TWO-POINT CORRELATORS

We now show that two-point functions of the local operators $A(r)$ and $S(r)$ are consistent with scale invariance,

$$ \langle A(0)A(r) \rangle = \frac{\alpha_A^2}{r^{2x_A}}, \quad \langle S(0)S(r) \rangle_{\text{conn}} = \frac{\alpha_S^2}{r^{2x_S}}. \quad (26) $$

Figure 17 compares data for the critical two-point functions to such power-law fits, giving good agreement at larger separations. The exponents $x_A$ and $x_S$ in the fits have been fixed to the values 1.224 and 1.506, respectively (see Table I), while the nonuniversal constants $\alpha_{A,S}$, which we will require in Sec. VII, have been left free.

<figure><img src="image_19.png" /><figcaption>FIG. 17. Two-point correlators for the operators $A$ (top) and $S$ (bottom). The displacement between the two operators is taken parallel to a lattice direction. Dashed lines are fits of the $L = 64$ data in the range $r \in [10, 15]$ to the forms in Eq. (26), with $x_A = 1.224$, $x_S = 1.506$ fixed and $\alpha_{A,S}$ free, giving $\alpha_A^2 = 0.72$ and $\alpha_S^2 = 0.77$.</figcaption></figure>

The simulations also give access to dynamical correlation functions in Monte Carlo time, which we analyze in Sec. VIII. These also encode the exponents $x_A$ and $x_S$, together with a dynamical exponent $z$.

## VII. THREE-POINT FUNCTION AND CONFORMAL INVARIANCE

Conformal invariance fixes the three-point functions in terms of the fields' scaling dimensions and operator product expansion (OPE) coefficients [91]. Conversely, data for three-point functions allow a direct numerical test of conformal invariance.

The OPE coefficients for the fields $A$ and $S$ that are allowed by duality symmetry to be nonzero are $C_{AAS}$ and $C_{SSS}$. Here, we examine the three-point function $\langle A(0)A(\mathbf{r})S(\mathbf{r}') \rangle_{\text{conn}}$ and give a very rough estimate of the corresponding OPE coefficient $C_{AAS}$. Data for $\langle S(0)S(\mathbf{r})S(\mathbf{r}') \rangle_{\text{conn}}$ were too noisy for a similar analysis.

The form dictated by conformal invariance for the three-point function is

$$ \langle A(0)A(\mathbf{r})S(\mathbf{r}') \rangle_{\text{conn}} = \frac{C_{AAS} \times \alpha_A^2 \alpha_S}{|\mathbf{r}|^{2x_A - x_S} |\mathbf{r}'|^{x_S} |\mathbf{r} - \mathbf{r}'|^{x_S}}, \quad (27) $$
where $\alpha_{A,S}$ are the same operator normalization constants that appear in the two-point functions (26).

We consider four possible spatial arrangements for the three points in the correlator, lying either on a line (L) or on the vertices of a right triangle (D):

$$L_{AAS}(r) = \langle A(0, 0, 0)A(r, 0, 0)S(2r, 0, 0) \rangle_{\text{conn}}, \quad (28)$$

$$L_{ASA}(r) = \langle A(0, 0, 0)S(r, 0, 0)A(2r, 0, 0) \rangle_{\text{conn}}, \quad (29)$$

$$D_{AAS}(r) = \langle A(0, 0, 0)A(r, 0, 0)S(r, r, 0) \rangle_{\text{conn}}, \quad (30)$$

$$D_{ASA}(r) = \langle A(0, 0, 0)S(r, 0, 0)A(r, r, 0) \rangle_{\text{conn}}. \quad (31)$$

We use ratios of these three-point functions to test for conformal invariance. The CFT prediction depends only on $x_A$ and $x_S$ (and on the arrangement of points in the correlator), so this test does not require an independent estimate of the nonuniversal constants $\alpha_{A,S}$ in Eq. (27).

Figure 18 compares each of the three independent three-point function ratios with the CFT prediction (for $x_A = 1.224$, $x_S = 1.506$), which is marked with a dashed line. Modulo uncertainty in the exponent estimates, the data should converge to these lines at large $r$. Statistical errors limit us to small $r$ because of the rapid decay of the three-point functions with $r$. Even so, there is agreement, within errors, with the CFT prediction once $r \gtrsim 3$.

Motivated by this consistency, we make a very preliminary estimate of the universal constant $C_{AAS}$. Figure 19 shows finite-$r$ estimates obtained from Eq. (27) (for each geometry of the three-point functions). The data suggest that $C_{AAS} \sim 1.5$. The uncertainty is large because of the very small range of $r$ and because the uncertainty in $\alpha_{A,S}$

<figure><img src="image_20.png" /><figcaption>FIG. 18. Test of conformal invariance: ratios of the three-point correlators defined in the text. Assuming exponent values $x_A = 1.224$, $x_S = 1.506$, conformal invariance requires these ratios to converge at large $r$ to the values indicated by dashed lines (error bars are from variation between six samples). We find agreement with the predicted value, within error bars, once $r \gtrsim 3$. Error bars become too large for a useful comparison once $r > 4$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIG. 19. Main panel: finite-$r$ estimates of the OPE coefficient $C_{AAS}$ using three-point functions with four different geometries, using data from system size $L = 48$ (error bars are from variation between six samples). Inset: average of the four estimates for $L = 48$ and also for $L = 32$.</figcaption></figure>

(obtained from the two-point function in Sec. VI) is hard to estimate. It would be worthwhile to improve this estimate. References [92,93] discuss methods for numerical estimation of OPE coefficients.

## VIII. STOCHASTIC DYNAMICS OF MEMBRANES

So far, we have discussed the gauge theory as a problem of equilibrium statistical mechanics in either 2 + 1 or 3 + 0 dimensions. But our simulations involve a fourth coordinate, which is Monte Carlo time (denoted $t$). The Monte Carlo dynamics may be interpreted physically as a model for stochastic thermal motion of classical membranes (Sec. II B), or, alternately, of classical spins in 3D (Sec. II A). These dynamics contain additional universal data beyond the data in static correlations: most importantly, the dynamical exponent $z$ that dictates how the typical relaxational timescale $\tau$ scales with system size $L$ at the critical point, $\tau \sim L^z$ [94]. Two-time correlation functions in this dynamics are also an alternative means of determining the exponents $x_A$ and $x_S$, as shown below.

### A. Universal dynamics and duality

There is great freedom in the microscopic definition of the stochastic dynamics, i.e., the choice of update for our Monte Carlo Markov chain. But we expect to find a dynamical fixed point that embraces a large class of microscopic updates that are local and preserve detailed balance (our updates are local and are described in Sec. VIII C below). This case is analogous to, say, the critical 3D Ising model, which shows a robust universality class for spin-flip dynamics with no conservation laws (the universality class of “Model A” [95–110]).

As in the Ising model, the dynamical universality class may change if we introduce conservation laws
[96]. For example, dynamics that conserve the total membrane area (the total number of occupied plaquettes) could be relevant to some experimental settings. The dynamical universality class may also change if we include nonlocal updates in the Monte Carlo simulations: Finding a nonlocal update that speeds up simulations by reducing z is a challenging open problem (Sec. XI).

The present dynamical critical point has one subtlety that arises from self-duality. We have argued that self-duality is a $\mathbb{Z}_2$ symmetry of the 3D fixed point, allowing us to classify scaling operators as $\mathbb{Z}_2$ even or odd (S and A, respectively). The mixed correlator $\langle AS \rangle$ therefore vanishes in the equilibrium ensemble [111]. But the Monte Carlo dynamics itself is not $\mathbb{Z}_2$ symmetric [12]. To define the dynamics, we had to choose one of the two dual representations, either on the original cubic lattice or on its dual, breaking the symmetry between them. As a result, the mixed correlator $\langle AS \rangle$ can be nonzero for nonequal times.

Assuming that the scaling operators $A(r)$ and $S(r)$ of the three-dimensional theory are lifted to scaling operators $A(r, t)$ and $S(r, t)$ in the dynamical theory, standard dynamical scaling [96] at large $|r|$ and $t$ gives [112]

$$ \langle A(r, t)A(0, 0) \rangle = t^{-2x_A/z} F_{AA}(t/|r|^z, t/L^z), \quad (32) $$

$$ \langle S(r, t)S(0, 0) \rangle = t^{-2x_S/z} F_{SS}(t/|r|^z, t/L^z), \quad (33) $$

$$ \langle A(r, t)S(0, 0) \rangle = t^{-(x_A+x_S)/z} F_{AS}(t/|r|^z, t/L^z). \quad (34) $$

The $\mathbb{Z}_2$ symmetry of the equilibrium critical point ensures that the last line vanishes at equal time.

## B. Dynamical scaling collapse

First, we obtain the dynamical exponent from the typical relaxation timescale $\tau(L)$ of a sample of size $L$. We estimate this timescale from the exponential decay of various two-time correlators, in particular, those of S and A (inset to Fig. 20). The various estimates are consistent with each other [113], and fitting $\tau(L)$ to a power law gives (the error bar is purely statistical)

$$ z = 2.48(10). \quad (35) $$

For comparison, this is larger than the dynamical exponent for spin-flip dynamics in the 3D Ising model [97–110], for which a recent estimate is $z = 2.0245(15)$ [110]. The exponent is closer to that for Metropolis dynamics at the confinement transition of pure $\mathbb{Z}_2$ gauge theory (corresponding, by duality, to plaquette flip dynamics in the ensemble with $x = 1$), which was estimated as $z = 2.70(3)$ [114].

The main panel of Fig. 20 demonstrates scaling collapse for the temporal correlators of the spatially averaged operators A and S, using this exponent. [The relevant scaling forms are given by integrating Eqs. (32)–(34).] Results are consistent with expectations from Sec. VIII A,

<figure><img src="image_22.png" /><figcaption>FIG. 20. Main panels: scaling collapse of autocorrelation functions for $\langle A(0)A(t) \rangle$ (top), $\langle S(0)S(t) \rangle$ (center), and $\langle S(0)A(t) \rangle$ (bottom) as a function of $t/L^z$ (using $x_A = 1.224$ and $x_S = 1.506$). Inset: autocorrelation time as a function of system size for four different correlators—$\langle E(0)E(t) \rangle$ (blue triangle), $\langle A(0)A(t) \rangle$ (orange triangle), $\langle S(0)S(t) \rangle$ (green pentagon), and $\langle S(0)A(t) \rangle$ (red star). The straight line fits all the data points to a power law $\tau = AL^z$, with $z = 2.48$.</figcaption></figure>

including the continuous vanishing of the scaling function for $\langle S(0)A(t) \rangle$ as $t \to 0$.

## C. Monte Carlo updates

The simplest Monte Carlo update is one that flips the state of a single plaquette with the appropriate Metropolis probability. However, a notable feature of configurations close to the multicritical point is that only a very small fraction ($\approx 2.5\%$) of links are occupied. When occupied links are rare, an attempted plaquette update has a high chance of creating four new occupied links, significantly increasing the energy, and therefore a high chance of being rejected.

This fact suggests that while plaquette updates are necessary for allowing occupied links to move, they are inefficient at moving surfaces around. To speed up the equilibration of surfaces, we therefore combine plaquette updates with a second update that flips the state of all six surfaces of a cube. Since this move never changes the number of occupied links, it does not face the problem above, and since it is still a local update, we do not expect it to change z, as we have confirmed numerically. See Appendix C for further details, including the scheme for parallelization.

# IX. MEMBRANE PATCHING, EMERGENT 1-FORM SYMMETRY, AND WORLDLINE PERCOLATION

In this section, we widen our focus to more general transitions out of the deconfined phase, including those
away from the self-dual line. We give a construction for the “fictitious” Ising order parameters (which are the key feature of Ising* transitions) on the Higgs and confinement lines. We find that these fictitious order parameters can be constructed all the way along the Ising* lines but not at the self-dual critical point. The disappearance of the fictitious order parameter at the self-dual critical point is associated with the emergence of a scale-invariant ensemble of loops there.

Studying this ensemble of loops with percolation-like observables [6,91] allows another numerical test of scale invariance at the self-dual critical point, and it gives another critical exponent with which to characterize it (Secs. IX B and IX C).

### A. Patching membranes

The fact that the *e* and *m* condensation transitions have Ising exponents (away from the self-dual line) is easy to understand at the boundaries of the phase diagram (Fig. 3), as reviewed in Sec. II E [3,11]. In these limits, the partition function can be written as a sum over closed membrane configurations. Mapping these closed membranes to Ising domain walls gives the relation to Ising.

Moving away from this extreme limit, the membranes acquire “holes” [6]. (We use the term “hole” loosely: More precisely, we mean any connected cluster of links in $∂M$, as defined in Sec. II B.) But it is natural to think that, if these holes have a finite typical size, coarse-graining beyond this size may restore a picture in terms of closed membranes that can be interpreted as Ising domain walls [115]. This is a heuristic explanation for why an effective Landau theory is useful even some distance away from the boundary of the phase diagram.

Here, we first want to give an explicit construction of these emergent degrees of freedom. Our approach is simply to “repair” the membranes *M* in a given configuration. We will be schematic, deferring further details and a numerical demonstration to Ref. [78]. Second, we want to understand what happens to the fictitious order parameter when we move along the Ising* transition line towards the self-dual critical point. This issue is closely connected to the question of where in the phase diagram emergent 1-form symmetries [33–36] exist.

We consider the membrane ensemble in Eq. (3), which is convenient for describing one of the two dual 1-form symmetries. By duality, analogous considerations apply for the dual symmetry. (The Ising* transition that we discuss below is the *m* condensation line.)

Let us briefly make the connection with a formal point of view. The fictitious Ising order parameter will make sense if (perhaps after coarse-graining) we can consistently define string operators $V_P = \pm 1$, supported on arbitrary paths *P* in spacetime, that count the parity of the number of membranes that intersect *P*. Let us assume that we can define such operators, which are functions of the membrane

configuration *M* and whose value is unchanged if the shape of the path *P* is deformed (while preserving the locations of the endpoints, if *P* is not a closed loop). We can then define a coarse-grained Ising variable $\phi_r$, up to a global $\mathbb{Z}_2$ ambiguity, by identifying $\phi_r\phi_{r'}$ with $V_P$ for any path *P* between *r* and *r'*. (Here, for simplicity, we consider an infinite system [117].)

Such string operators, obeying an invariance under deformations, define a $\mathbb{Z}_2$ 1-form symmetry (see Refs. [35,36], and Refs. [33,118,119], where the alternate term “gauge-like symmetry” is used, for definitions) [120]. Switching briefly to the language of 2D quantum states, the analogous quantum operators in the toric code are simply the familiar topological string operators [22,33], which can be used to create pairs of *m* anyons at their endpoints (a similar dual operator creates pairs of *e* anyons). Perturbing away from the solvable limit of the toric code, dressed versions of these string operators are expected to exist, in principle, so long as the other anyons, which braid non-trivially with *m*, remain gapped [36,121–123].

Returning to the membrane picture, how would we explicitly define such string operators, or the $\phi_r$ configuration, in a simulation? A natural approach is to start with the membrane configuration, and try to “patch up” the holes, to give closed membranes. This approach is not a strictly local process since holes can be of any size. We also have some freedom in the convention, or algorithm, for constructing the patching surfaces (one possible convention is given in endnote [125]). But, if holes have a finite typical size, and large holes are exponentially rare, we expect the nonlocality in the patching operation to be mild. Each finite “loop” (cluster of links) in $∂M$ may be patched by attaching a finite surface of comparable size. This is illustrated in Fig. 21.

Once this is done, we may define $\phi_r$ (again with a global $\mathbb{Z}_2$ freedom). Because of the nonlocality of the patching operation, this effective field is only likely to be a useful concept on length scales larger than the typical size of a

<figure><img src="image_23.png" /><figcaption>FIG. 21. Membrane patching process (schematic). A configuration *M* of membranes, with a nonempty boundary $∂M$ made up of finite clusters of occupied links (shown in red) is “patched” by attaching finite surfaces (blue) to the components of $∂M$. In the resulting ensemble of *closed* membranes, line operators may be defined, specifying an emergent 1-form symmetry. A fictitious scalar field $\phi = \pm 1$ may also be introduced, whose domain walls are the patched membranes. The Ising* transition is the ordering transition for this fictitious field.</figcaption></figure>
loop. (It will therefore be most useful as an effective field in a Lagrangian when it has a correlation length that is parametrically larger than this loop size, or infinite.) This construction allows us to compute the “two-point” correlation function of $\phi_r$ in a simulation and to extract the corresponding anomalous dimension, despite the fact that $\phi_r$ is not a local gauge-invariant quantity [78].

We may also define thickened string operators. Let $P$ be, say, a straight path of length $\ell \gg 1$. In order to determine how many domain walls $P$ pass through after patching, we must check how many loops $P$ thread in the unpatched configuration $\mathcal{M}$. This check requires us to examine a cigar-shaped region around $P$, wide enough to contain (with probability close to one) all the loops that $P$ threads. The operator $V_P$ is therefore a function of the degrees of freedom within this cigar-shaped region. When large loops are exponentially suppressed, the largest loop that $P$ threads will typically be of size of order $\ln \ell$ (due to rare large loops), so the typical width of the cigar should be of this order. However, close to the endpoints of $P$, it is sufficient for the width to be only somewhat larger than the typical loop size (i.e., $\ell$ independent).

The above pertains to the case where the “holes” have a finite typical size $\xi_h$. If, on the other hand, $\xi_h$ diverges, so that samples of arbitrarily large size $L$ contain holes of size comparable with $L$, then this procedure for defining the string operator and effective Ising order parameter is liable to fail (since $\phi_r$ and $V_P$ can become highly nonlocal). In other words, a sufficient [126] condition for this procedure to work is that the appropriate set of worldlines, $\partial \mathcal{M}$, is in the nonpercolating phase when viewed as a bond percolation configuration.

For this reason, it is interesting to revisit the question of where in the phase diagram these worldlines percolate [6], which we do next. For example, as we move along the confinement transition line, starting at $y=0$ (where there are no worldlines) and moving towards the self-dual critical point, where does the fictitious order parameter $\phi_r$—as defined by the above simple algorithm—stop making sense? The results below indicate that it makes sense all the way along the confinement line but not at the self-dual critical point where that line terminates. They are consistent with the simplest expectation: that the 1-form symmetry that exists at $y=0$ persists as an emergent symmetry all the way along the confinement transition line but disappears at the self-dual critical point. Therefore, the RG flow from the self-dual fixed point to the Ising* fixed point involves the emergence of the 1-form symmetry. Similar considerations apply for the dual 1-form symmetry along the Higgs line.

## B. Percolation summary

Our result for the percolation phase diagram is shown in Fig. 22 and explained in Sec. IX C. Within our numerical precision, the percolation phase boundary matches the thermodynamic boundary of the deconfined phase along the entire Higgs transition line (to the right of the self-dual

<figure><img src="image_24.png" /><figcaption>FIG. 22. Phase diagram showing the phase boundary for percolation of $e$ worldlines [defined as clusters of links in $\partial \mathcal{M}$; see Eq. (3)]. The deconfined phase lies within the non-percolating phase. The self-dual line is shown as a dashed line.</figcaption></figure>

line) and passes through the self-dual critical point. (The percolation transition line also lies very close to the first-order line, though closer examination indicates that these two lines do not entirely coincide; see Appendix D.)

The fact that the percolation line passes through the self-dual critical point agrees with the scenario in Ref. [6]. A far as we are aware, however, this result is not guaranteed *a priori*: The geometrical percolation transition could have separated from the Higgs transition at some point along the Higgs line, with the multicritical point lying in the interior of the percolating phase (see Ref. [127]).

It is also striking that the self-dual critical point lies on the percolation phase boundary despite having a very low fraction of occupied links, around 2.5%. Despite their low density, these links make up a scale-invariant ensemble of clusters. Figure 23 shows the loops $\partial \mathcal{M}$ in an example configuration.

This scale invariance allows us to define a new exponent at the self-dual critical point, namely, the fractal dimension $d_f$ of the critical loops. *A priori*, this exponent is independent of the scaling dimensions of local operators discussed above: $d_f = 3 - x_{\text{conn}}$ is determined by the scaling dimension $x_{\text{conn}}$ of a nonlocal geometrical operator of the type familiar from percolation [91,129]. Interestingly, though, our numerical result for $d_f$ below (Sec. IX C),

$$d_f = 1.77(2), \qquad (36)$$

is consistent with $x_{\text{conn}} \stackrel{?}{=} x_A$, perhaps hinting at additional hidden symmetry structure at this critical point. (See Sec. X A for an argument that $x_A \le x_{\text{conn}}$.)

## C. Percolation observables

We locate the boundary between percolating and short-loop (nonpercolating) phases using the spanning
<figure><img src="image_25.png" /></figure>

FIG. 23. Loops in a sample of system size $L = 32$ at $x = 0.6367$. A spanning loop is highlighted in red (note periodic B.C.s). Our definition of a loop allows for junctions where four or six occupied links meet, but as can be seen, they are relatively rare at the critical point.

probability $P_s$, which is the probability that the sample contains a loop that spans the sample in a given axis direction [130]. In the thermodynamic limit, this quantity converges to 0 and to 1 in the nonpercolating and percolating phases, respectively, and it is expected to take a universal value in between 0 and 1 at a continuous transition between the two phases.

We estimate the percolation phase boundary from crossings in $P_s$, plotted as a function of $y$, using small system sizes $L = 8, 12, \text{and } 16$ (data not shown). We use larger sizes to analyze the transitions at $x = 1$, at $y = 1$, at the multicritical point, and in the region around it. We may also obtain the correlation length exponent from a scaling collapse of $P_s$.

The phase diagram in Fig. 22 shows three different phases. The deconfined phase has short loops, while the thermodynamically trivial phase splits into a percolating and a nonpercolating phase (this is possible because the percolation transition need not have any thermodynamic signature). Note that here we are considering percolation of *e* worldlines: The phase diagram for percolation of *m* worldlines (in the dual membrane representation) may be obtained by duality.

As a check, we first examine the percolation transitions on the boundary of the phase diagram, where we expect to see standard universality classes (data in Appendix D). At $x = 1$, results are as expected from the Ising mapping, with a fractal dimension consistent with the known value for critical Ising worldlines [131,132] and a correlation length exponent consistent with the Ising value. At $y = 1$, where the percolation transition is purely geometrical (has no thermodynamic signature), exponents are consistent with the standard 3D percolation universality class.

<figure><img src="image_26.png" /></figure>

FIG. 24. Scaling collapse of $P_s$ as a function of $(x - x_c)L^{1/\nu}$ with $x_c = 0.63666$ fixed. Inset: flow of $P_s$ as a function of the system size for several $x$ values close to the multicritical point, showing approximate scale invariance at $x = 0.6367$.

Figure 24 shows data on the self-dual line, close to the self-dual critical point. The data are compatible with a critical point very close to $x = 0.6367$ (inset of Fig. 24), i.e., with geometrical criticality coinciding with the self-dual critical point at the corner of the deconfined phase. A scaling collapse of $P_s$ as a function of $(x - x_c)L^{1/\nu}$, leaving $x_c$ and $\nu$ free (not shown), gives $x_c = 0.63664(10)$ and $\nu = 0.69(6)$, compatible with our best estimates for the self-dual multicritical point. In Fig. 24, we show the scaling collapse when $x_c$ is fixed to our previous best estimate $x_c = 0.636660$ (Sec. VA). We use B-splines with five knots and obtain $\nu = 0.70(6)$ for a fit that gives $\chi^2 = 27$ for 24 degrees of freedom. Using $\nu^{-1} = 3 - x_S$, this result for $\nu$ is consistent with our previous estimate of $x_S$, though with lower precision.

At the self-dual critical point, loops are fractal and exist on all scales (Fig. 23). The fractal dimension $d_f$ can be estimated from fitting the total mass of the largest loop to a power law in $L$, Fig. 25. The straight line fits the whole

<figure><img src="image_27.png" /></figure>

FIG. 25. Mass of the largest cluster (number of links) as a function of the system size at $x = 0.6367$. The dashed line fits system sizes as a power law $AL^{d_f}$, with $d_f = 1.77(2)$.
range of system sizes from 8 to 64, providing the estimate $d_f = 1.77(2)$ quoted above.

## X. RELATED MODELS

The previous section concludes our analysis of numerical data. We now consider some variations of the model and relations to other models. Section XA connects the self-dual critical point to another partition function, for a “topologically constrained” ensemble of loops, which may be interesting to study further. Sections XB and XC discuss perturbations and crossovers in the gauge-Higgs model. Section XD discusses our numerical observation that the exponents $x_A$ and $x_S$ are close to exponents in the XY model, as well as issues related to continuum field theory.

### A. An unusual self-dual loop model

In Sec. II C, we considered a representation of the gauge-Higgs partition function as a loop model, for two species of loops [133], with a topological sign factor $(-1)^{\text{linking}}$ in the Boltzmann weight. In that model, the two species of loops live on distinct cubic lattices. Here, we consider a modified loop model in which the loops live on the same cubic lattice. This change allows the partition function to be reexpressed in a form involving a topological constraint rather than a topological sign factor.

Let $C_e$ and $C_m$ be two species of loops on the cubic lattice. Here, we define the allowed loop configurations differently than those in Sec. II C: Now, we insist that the loops are strictly self-avoiding and mutually avoiding (a loop may visit a given site once at most, for example). With this definition, the linking number $\hat{X}$ is well defined.

The partition function is

$$Z_{\text{mod}} = 4 \sum_{C_e, C_m} y^{|C_e|} y^{|C_m|} (-1)^{\hat{X}(C_e, C_m)}. \quad (37)$$

($|C_e|$ is the number of occupied links in $C_e$, etc., and the loops in $C_e \cup C_m$ are mutually avoiding.) We do not yet know the full phase diagram of this new model, but it is plausible that it may also show a self-dual critical point, in the same universality class as the lattice gauge theory studied above.

As an aside, we note that the model could be varied in many ways. We could allow the loops to be clusters [as in the previous model, Eq. (5)], by allowing the number of occupied links adjacent to a site to be any even number [rather than just 0 or 2 as in Eq. (37)]. With this choice, the model maps onto the original gauge-Higgs model in the limit $y \to 0$ and in the limit $y' \to 0$. This choice may have advantages for simulations (as may other choices of lattice as noted below). These changes do not affect the points we make here, so we consider the more easily visualized ensemble of strictly self-avoiding loops.

Let $C$ denote the full loop configuration, without regard to species labels, and let us specialize to the self-dual line where $y = y'$:

$$Z_{\text{mod}} = 4 \sum_C y^{|C|} \sum_{\substack{\text{species} \\ \text{labels}}} (-1)^{\hat{X}(C_e, C_m)}. \quad (38)$$

The final sum is over assignments of the loops in $C$ to species $e$ or $m$, i.e., over splittings of $C$ into $C_e$ and $C_m$. For simplicity, let us choose (nonperiodic) boundary conditions such that loops cannot end on the boundary or wind around the system.

We can sum over the species assignments explicitly, for a fixed $C$. The result is simple (see endnote [134] for the argument):

$$Z_{\text{mod}} = 4 \sum_C y^{|C|} \times 2^{(\text{loops in } C)} \times \chi_C. \quad (39)$$

Here, $\chi_C = 0, 1$ depends only on the topology of $C$ and simply imposes a restriction (constraint) on the allowed topologies: $\chi_C = 1$ so long as every loop in the configuration links with an even number of other loops, and $\chi_C = 0$ otherwise. This topological constraint is crucial. [Removing it, by removing the factor $\chi_C$ from Eq. (39), leaves the partition function for a version of the XY model [135].]

Strikingly, the expression in Eq. (39) is sign-free [unlike Eq. (38)] and could be sampled with Monte Carlo methods, using a local update that preserved the mod 2 total linking number of each loop. It would be interesting to know the phase diagram of this model or variants of it. (For an efficient numerical study, it might be useful to modify the lattice geometry of the model so that loops can form nontrivial links on a shorter length scale [136].)

This model also allows an interesting topological interpretation for correlation functions of the anti-self-dual operator.

In the ensemble (37), let us define the operator $\tilde{A}(r)$ at a site $r$ to take the value 0 if the site is not visited by a loop, 1 if the site is visited by an $e$ loop, and $-1$ if the site is visited by an $m$ loop. This operator is odd under duality, so it is analogous to the operator $A(r)$ defined for the gauge-Higgs model in Sec. III A.

By again explicitly summing over the loops' species labels, we may write correlators of $\tilde{A}$ in the formulation of Eq. (39). First, an insertion of $\tilde{A}(r)$ forces a loop to pass through $r$. Second, the $\tilde{A}(r)$ insertion forces the total linking number of this anchored loop (with other loops) to reverse its parity. In the original ensemble (39), every loop has even linking. In the presence of $\tilde{A}$ insertions, the linking number of a loop that passes through an odd number of $\tilde{A}$ operators must instead be odd (while the linking number of a loop that passes through an even number of $\tilde{A}$ operators remains even).
<figure><img src="image_28.png" /><figcaption>FIG. 26. Schematic: some of the configurations contributing to the correlator $\langle \tilde{A}(r) \tilde{A}(r') \rangle$ in the model (39). Circles indicate locations $r, r'$.</figcaption></figure>

In Fig. 26, we illustrate some of the configurations that contribute to $\langle \tilde{A}(r) \tilde{A}(r') \rangle$. For simplicity, we show the schematic situation at small $y$, where the loop length is suppressed. The first term involves a single loop (with zero linking) that passes through both insertions. The other terms involve a chain of linked loops, with the loops at the two ends of the chain, passing through $r$ and $r'$, having odd linking. (This picture shows that, in the limit of small $y$, the correlation length $\xi$ for this correlator is proportional to $|\ln y|^{-1}$ since we must pay a factor of $y$ for every unit of loop length.) Moving away from the regime of small $y$, it remains true that there are two types of terms: those where the two $\tilde{A}$ insertions lie on the same loop and those where they lie on distinct loops that (unlike all the other loops in the configuration) have odd linking.

Assuming this model has a self-dual critical point, then this topological picture for $\langle \tilde{A} \tilde{A} \rangle$ yields an inequality for the fractal dimension $d_f$ of loops, which sheds some light on the coincidence of exponents that we found numerically in Sec. IX. Recall that the connectivity correlator, $P_{\text{conn}}(r, r') \sim r^{-2x_{\text{conn}}}$, is the probability that two distant sites are connected by a loop. Above, we have shown that these connected configurations are a subset of the configurations that contribute to $\langle \tilde{A} \tilde{A} \rangle$. In other words, $\langle \tilde{A} \tilde{A} \rangle = P_{\text{conn}} + R$, where $R$ is the sum over the remaining configurations, which is positive. Therefore,

$$x_A \le x_{\text{conn}}, \quad (40)$$

or, equivalently, $d_f \le 3 - x_A$.

We expect that when the model above is perturbed away from the self-dual line, an Ising* transition can take place, as in the original gauge-Higgs model. In the language of Eq. (37), this occurs in the same manner discussed in Sec. II E: If one of the loop species has a finite typical size, it can be integrated out at large scales, leaving a simple ensemble of Ising-like worldlines (topologically unconstrained loops with a fugacity of 1 per loop). It is also possible to see this crossover in the language of Eq. (39): Adding the duality-breaking perturbation $\tilde{A}$ relieves the topological linking constraint in Eq. (39) and leads to an ensemble where large loops have a fugacity of 1 [137].

## B. Perturbations of the gauge-Higgs model

After this detour, we return to the standard gauge-Higgs model to discuss some remaining questions.

We have characterized the leading self-dual and anti-self-dual scalar (spin-zero) operators at the self-dual critical point numerically, but it remains to characterize the sub-leading operators in these sectors, as well as operators with higher spin. One motivation for this is to formally determine the number of relevant scaling parameters once duality is broken, as we explain below.

On the appropriate line in parameter space, self-duality is an exact property of the standard gauge-Higgs model. But in many settings where the gauge-Higgs model is a useful effective theory, exact self-duality will be broken in the ultraviolet by additional interactions. It is natural to conjecture that the phase diagram structure in Fig. 1 can nevertheless survive, with self-duality appearing as an emergent symmetry at the corner of the deconfined phase, where Higgs and confinement transitions meet. In order for this to be the case, $A$ and $S$ should be the only relevant scalar operators at the self-dual critical point.

At first glance, this case is demonstrated by the fact that we only had to tune two parameters to reach this critical point. However, this is not quite correct: The microscopic self-duality symmetry of the self-dual line forces all anti-self-dual perturbations to vanish there (not only the leading $A$ perturbation). Therefore, in principle, we should separately check whether the subleading duality-odd scalar operator is relevant or irrelevant. Since $A$ itself has a large scaling dimension, we might expect that this subleading operator will be irrelevant, but this should be checked.

The subleading duality-even operator is irrelevant, but a sufficiently large duality-even perturbation may yield a “self-dual tricritical” point with an additional relevant direction.

In Ref. [77], it was argued, using series expansions, that the toric code with $X, Y$, and $Z$ fields had a critical line, with varying exponents, in the $h_x = h_z$ plane, which will be interesting to investigate further, as continuously varying exponents in 3D are rare. However, it should be noted that, in the present language, the perturbation $h_y$ breaks both internal and spatiotemporal symmetry. The toric-code Hamiltonian with $h_x = h_z$, discussed in Sec. II D, has a duality symmetry $D$ that we may take to be $X \to T(Z)$, $Z \to T(X)$, $Y \to -T(Y)$, where $T$ represents a translation by $(1/2, 1/2)$. It also has an antiunitary time-reversal symmetry, which we may take to act as $X \to X$, $Z \to Z$, $Y \to -Y$, $i \to -i$. Adding the $h_y$ field breaks both of these symmetries. (It preserves their product.) It would be interesting to identify the leading continuum perturbation of the self-dual critical point that is induced by the $h_y$ coupling.

Recent work has demonstrated infinite-randomness scaling for a range of Higgs transitions in 2+1D quantum gauge theories with quenched disorder in the couplings [138]. It would be interesting to study the effect of disorder on the self-dual topological phase transition. The exponents $x_S$ and $x_A$ imply that spatially uncorrelated quenched
disorder is a strongly relevant perturbation of the self-dual
critical point in its 2 + 1D quantum manifestation, regard-
less of whether this disorder preserves duality or not [139].

C. Dimensional crossovers

Various dimensional crossover effects may also be worth studying. By making one of the three lattice directions finite and of width $1/T \gg 1$ (with periodic boundary conditions), we may study the effect of a low but nonzero temperature in the quantum problem [91]. Standard considerations show that on the boundaries of the phase diagram (at $y=0$ or $y'=0$), the 3D Ising* transitions give way to 2D Ising transitions, but in the interior of the phase diagram, these transitions become crossovers, with a finite correlation length [140]. This correlation length is exponentially large in $1/T$ at small $T$. In the worldline representation (5), this scaling is associated with closed worldlines of the massive anyon, which are of length $1/T$ and wrap around the temporal cycle [141]. (This exponential scaling may be why a numerical study of the $\mathbb{Z}_2$ gauge-Higgs model instead reported finite-temperature transitions [142].) A line of first-order transitions on the self-dual line will remain at small finite temperature. What happens to the self-dual critical point at nonzero temperature is less clear. The simplest possibility is that it becomes a conventional critical endpoint so that the interior of the phase diagram contains only a first-order line, bounded by two conventional critical endpoints.

Other boundary condition choices for the finite dimension give other phase diagrams. For example, consider the loop model (5) in a slab geometry of thickness $\ell$, with open boundary conditions in the finite direction, with loop strands forbidden from terminating on the boundary. Physically, this can be obtained by taking a 2D quantum system deep in the deconfined phase (corresponding to $y, y' \ll 1$) and then varying the couplings inside a strip of width $\ell$ to allow anyons to proliferate there. Isolated strands that span the finite direction are now forbidden, so the mechanism that rendered the correlation length finite in the previous quasi-2D geometry is removed. Instead, after coarse-graining to scales much larger than $\ell$, we may argue for an effective 2D loop model with two species of loops (and with no nontrivial topological sign factor). Away from the self-dual line, 2D Ising transitions are likely, which are associated with proliferation of a single loop species (a single anyon type). It may also be possible to have a gapless Ashkin-Teller-like regime on part of the self-dual line, where both species are critical.

If we think of the length-$\ell$ direction here as imaginary
time instead of space, this setup may be related to the
interesting finding of 2D Ashkin-Teller criticality in a
deformed toric-code wave function, for which equal-time
correlators map to correlators in a 2D classical model [143].
This deformed wave function is given by a finite-depth
nonunitary circuit acting on the toric-code wave function.
This can be visualized in a path integral representation. In

the zero-temperature path integral for the deconfined phase,
the evolution for imaginary times ($-\infty, 0$) can be viewed as
preparing the ground-state ket and the evolution for times
$(0, \infty)$ as preparing the corresponding bra. To obtain equal-
time correlators in the deformed wave function, we insert a
"slab" of finite temporal extent in between these two pieces,
representing the action of the nonunitary circuit on the bra
and ket. This is reminiscent of the setup for the quasi-2D
loop model above (since in the spacetime region outside the
slab, we set $y = y' = 0$, meaning that, in the loop model
picture, worldlines are forbidden except inside the slab).

It will also be interesting to characterize boundary
critical phenomena, and conformally invariant boundary
conditions, for the self-dual topological transition.

D. Comparison with XY\* and other field theories

A striking feature of our numerical results is that the
values for scaling dimensions are close to certain values for
the 3D XY (starred) fixed point. Below, we discuss why this
is surprising, as well as some of the challenges for a field
theory treatment of the model.

At first sight—however, see below—a relationship with
XY appears to be a natural guess [6], by analogy with
conventional ordering transitions, where two Ising critical
lines (together with a first-order line) can meet at an XY
critical point. In the present model, flow to the XY fixed
point may be ruled out immediately, for example, by the
fact that the adjacent phases and the low-lying operators do
not match. But, in addition, “starred” (orbifolded) versions
of the XY fixed point may be ruled out—these also fail to
reproduce the universal physics.

Given two conventional Ising-like order parameters $\varphi_x$ and $\varphi_y$, and an additional $\mathbb{Z}_2$ symmetry that exchanges them, XY criticality for $\varphi = (\varphi_x, \varphi_y)$ can arise by tuning one parameter because the symmetry-allowed “cubic” anisotropy $\varphi_x^4 + \varphi_y^4 - 6\varphi_x^2\varphi_y^2$ is a (weakly) irrelevant operator at the XY fixed point [144,145].

In the present model, we know that the Ising* transition
lines, away from the self-dual point, can be understood as
ordering transitions for “fictitious” (non-gauge-invariant)
Ising-like order parameters. Therefore, at first sight, it is
tempting to make an analogy with the above Landau theory,
which would mean identifying the operator $S$ with the
thermal operator $\varphi^2$ and the operator $A$ with the symmetry-
breaking mass operator $\varphi_x^2 - \varphi_y^2$. The scaling dimensions of
these operators in the XY model are $x_{\varphi^2} = 1.51136(22)$
and $x_{\varphi_x^2 - \varphi_y^2} = 1.23629(11)$ [146–148]. Strikingly, the
differences between these values and our results for $x_S$
and $x_A$ in Table I are small, comparable in size with the
(statistical) error bars quoted in the table [149].

The basic problem with this analogy is that it ignores the
nontrivial mutual statistics between *e* and *m* excitations
[13,14,38] that are the key feature of the transition. These
mutual statistics do not affect critical exponents on the
Ising* lines because only one of the two excitations is massless on these lines. But both excitations become massless at the self-dual critical point.

For example, any consistent description of the fixed point should correctly reproduce the spectrum of low-lying anyonic quasiparticles that exists when we perturb slightly away from the self-dual critical point into the deconfined phase. It is hard to see how this could be consistent with a mapping that related the operators of the present fixed point to those of the XY fixed point. For example, starred versions of the anisotropic XY theory, in which we gauge the order parameters $\phi_x, \phi_y$ with flat gauge fields, manifestly fail to capture these mutual statistics.

The argument against a connection with (starred) XY can also be made in the geometrical pictures. In the membrane picture, the possibility of defining a fictitious Ising order parameter is associated with the membranes being effectively closed on large scales, as discussed in Sec. IX. But at the self-dual critical point, we have holes in these membranes on all scales, as we have demonstrated explicitly. Therefore, the attempt to make a connection with a simple Landau theory, at least in this manner, fails at this critical point.

Therefore, it seems very likely that the exponents $x_A$ and $x_S$ at the self-dual critical point are numerically close to XY exponents but distinct from them. (As we have seen, if the exponents were the same as those of the XY fixed point, this relationship between the operator spectrum of a topological phase transition and that of a simple ordering transition would have to be of a fundamentally new kind.) It would be interesting to have an explanation for this closeness of exponents in distinct universality classes [150].

Let us comment on another approach to a putative field theory. This is to use a U(1) gauge theory, with symmetry-breaking perturbations, as a “UV” starting point, hoping to flow to the fixed point of interest. Progress on field theories for deconfined criticality [10,44–55], including explicitly self-dual formulations [51,52] (for instance, Ref. [52] constructs a symmetry-enriched $\mathbb{Z}_2$ topological order that includes a self-duality within a larger symmetry group), provides some insight. However, a basic challenge in this approach is the need to eliminate global U(1) symmetries associated with U(1) gauge flux conservation, which are unphysical in the present context. This challenge requires monopole perturbations in the Lagrangian. These operators are not local functions of the fields in the Lagrangian, so they are hard to treat. Concretely, one possibility is to start with the well-known representation of $\mathbb{Z}_2$ topological order as a $U(1) \times U(1)$ gauge theory with a mutual Chern-Simons term [151–154]. Each gauge field may be coupled to a bosonic matter field, whose quanta are $e$ and $m$ anyons, respectively, and transitions out of the deconfined phase may be induced by varying the masses of these particles. The problem is that, as it stands, this field theory has global $U(1)$ (flux conservation) symmetries that are not present in the

microscopic model we have been discussing. To faithfully capture the symmetries of the problem of interest we must include the appropriate symmetry-breaking terms in the Lagrangian, which include (gauge-invariant) monopole operators. As noted above, we are currently short of techniques for treating field theories with such terms in the Lagrangian. Therefore, while the above gives a formal UV Lagrangian that may flow to the fixed point of interest, it is more challenging to construct a theory that is tractable or predictive. We will discuss possible means of sidestepping these obstacles to a field theory description elsewhere.

XI. OUTLOOK

The three-dimensional $\mathbb{Z}_2$ gauge-Higgs model is the simplest nontrivial lattice gauge theory [3,11,13,14,23,24,155]. Its remarkable duality property allows for a self-dual topological phase transition whose properties have long been unresolved. We have given direct evidence for scale invariance at this transition, exploring system sizes up to 2 orders of magnitude larger than the lattice spacing. Exciting directions remain open, on the computational, experimental, and theoretical fronts.

First, there are many intriguing questions that could be addressed using further simulations. At the basic level, armed with the accurate estimate of $x_c$, further characterization of the critical point will be possible, examining the scaling dimensions of a wider range of operators (Sec. X B) and pinning down OPE coefficients more precisely (Sec. VII). It will also be interesting to study correlators involving extended operators [6,78,124,156].

We have also proposed new models that could be simulated. The loop model in Sec. X A has a simplified action of self-duality. It has a sign-free reformulation of a nonstandard kind, as an ensemble of loops with a simple topological constraint. (This connects, heuristically, to the longstanding question from polymer physics of how to think about the renormalization group for models with topological constraints [60,157–165].) This sign-free formulation could be exploited to determine the model’s phase diagram and may suggest a more general strategy for obtaining sign-free lattice models for topological transitions.

In the context of the standard lattice gauge theory, a range of perturbations and crossovers may be studied (Secs. X B and X C), for example, to search for self-dual tricriticality.

The self-dual topological phase transition can be viewed as a paradigmatic challenge for Monte Carlo algorithm design. Although it is Monte Carlo sign-free (unlike many other lattice gauge theories [166–169]), the lack of a nonlocal cluster update [37] for ensembles of membranes and the large dynamical exponent (Sec. VIII) make it expensive to simulate. Creative algorithmic improvements would be valuable. We might consider updates acting on larger finite clusters, perhaps optimized using machine learning [170,171].
If we are in the deconfined phase but close to the self-dual critical point, various features of the spectrum of massive quasiparticles [14,77] will be universal and could perhaps be examined using Monte Carlo [172], series expansion [14,77], or tensor network techniques [173,174]. For example, does the fermionic $e$ excitation exist as a stable bound state in this regime, or does it inevitably decay into an $e$ and an $m$?

Even away from the self-dual point, interesting questions remain. The existence of fictitious Ising order parameters on the Higgs and confinement transition lines is the key to the theoretical understanding of these transitions [3,6,11,66]. We have argued that we can construct these field configurations explicitly by a quasilocal patching process in the membrane representation of the partition function so that, for example, the Ising “two-point function” $G(r, r')$ can be computed numerically. Formally, this is the expectation value of a dressed string operator that extends from $r$ to $r'$ (Sec. IX). In a separate work, we will analyze the emergence of this structure in more detail [78].

The self-dual critical point may be accessible experimentally, either in its 3D classical or its 2+1D quantum manifestation. It would be exciting to see the full structure of the gauge-Higgs phase diagram, with the meeting of the two Ising* lines, in experiments on amphiphilic membranes (verifying a longstanding conjecture [6]). In order to access this point, the membranes must have free edges, i.e., a nonempty membrane boundary $\partial M$. However, the results in Appendix B suggest that the required density of free edges may be relatively small.

Strategies for quantum simulation of lattice gauge theories are under intensive development [175–184], so it may one day be possible to explore the self-dual critical point and its real-time quantum dynamics experimentally.

Perplexing theoretical questions remain. Why are our estimates for $x_A$ and $x_S$ so close to XY values (Sec. X D)? Further numerical characterizations of the critical point mentioned above may shed light on this. Significant input may also come from the conformal bootstrap [185–188], by exploring the space of theories with the requisite $\mathbb{Z}_2$ symmetry.

There remains the fundamental question that we started with: Can we formulate a useful continuum field theory for the self-dual topological transition? Criteria for “usefulness” could include the possibility of calculating exponents in a systematic expansion, as well as the possibility of deriving the structure of phase diagrams analytically. More generally, the time seems ripe for a numerical and theoretical attack on phase transitions where multiple species of anyons, with nontrivial statistics, simultaneously condense [13,14,38,39].

## ACKNOWLEDGMENTS

We thank Hans Evertz, Paul Fendley, David Huse, Nabil Iqbal, Jack Kemp, John McGreevy, Miguel Ortuño,

Siddharth Parameswaran, Yang Qi, T. Senthil, Tin Sulejmanpasic, Sagar Vijay, and Chong Wang for useful discussions. A. S. acknowledges support by AEI (Spain)/ FEDER (EU) Grant No. PID2019-104272RB-C52/AEI/ 10.13039/501100011033. A. S. and P.S acknowledges support by Fundación Séneca Grant No. 19907/GERM/15. A. N. acknowledges support from the Gordon and Betty Moore Foundation under the EPiQS initiative (Grant No. GBMF4303), from EPSRC Grant No. EP/N028678/1, and from a Royal Society University Research Fellowship.

## APPENDIX A: MEMBRANE REPRESENTATION OF Z

We now review the standard relationship between the Ising gauge theory partition function and partition functions for membranes on either the original cubic lattice or its dual [6]. In the interpretation as a 2D quantum system in imaginary time, these membranes are world surfaces of either electric or magnetic strings (cf. Fig. 6), depending on whether we use the original lattice or the dual lattice. In general, the strings (which live in a 2D spatial plane) can be open lines, terminating at $e$ or $m$ particles in the respective cases. Therefore, the membranes (which live in 3D spacetime) are not closed, in general, but rather have boundaries, which are the worldlines of $e$ or $m$ particles, respectively. The “action” of a given membrane configuration (the logarithm of the Boltzmann weight) is, up to constants, the area of the world surfaces plus the length of the worldlines.

### 1. Membranes on the original lattice

Using the fact that the variables take only the values $\pm 1$, $Z$ in Eq. (1) can be rewritten in a form convenient for a standard graphical expansion:

$$Z(x,y) \equiv \frac{1}{2^{4L^3}} \sum_{\{\sigma\}, \{\tau\}} \prod_{\square} \left(1+x\prod_{\sigma} \sigma\right) \prod_{\ell} (1+y\sigma\tau\tau), \quad (A1)$$

with $K = \frac{1}{2}\ln(1+x/1-x)$ and $J = \frac{1}{2}\ln(1+y/1-y)$. We expand out the products over (1) plaquettes $\square$ and (2) links $\ell$ in Eq. (A1), and represent a given term by coloring plaquettes of the lattice and highlighting links in bold, as in Fig. 4. A plaquette is colored (“occupied”) iff we pick the “$x\prod\sigma$” term for that plaquette, and similarly, a link is bold if we pick the “$y\sigma\tau\tau$” term. For a given term in the expansion, the collection of occupied plaquettes constitutes the membrane configuration $M$.

Now, for each term, we must sum over $\sigma$ and $\tau$. The term will vanish if there is any link $\ell$ where the terms we have chosen contribute $\sigma_\ell$ an odd number of times in total, which means that the set of bold links must coincide with the membrane boundary $\partial M$ to have a nonvanishing term ($\partial M$ is defined as the set of links where an odd number of colored plaquettes meet). If this is satisfied, then the sums
over $\sigma$ and $\tau$ are both nonvanishing, giving trivial factors $2^{3L^3}$ and $2^{L^3}$, respectively, that cancel the normalization term chosen in Eq. (A1). We are left with the partition function as a sum over membrane configurations weighted by $x^{|M|}y^{|\partial M|}$, Eq. (3).

From this expansion, we also see that the face and edge operators defined in Sec. III may be written as

$$ \mathcal{F}(p) = \frac{x}{1-x^2} \left( \prod_{\ell \in p} \sigma_{\ell} - x \right), \quad (A2) $$

$$ \mathcal{E}(\ell) = \frac{y}{1-y^2} (\sigma_{ij}\tau_i\tau_j - y), \quad (A3) $$

for a plaquette $p$ and link $\ell = \langle ij \rangle$, respectively. For example, $\mathcal{F}(p)(1-x\prod_{\ell \in p} \sigma_{\ell}) = x\prod_{\ell \in p} \sigma_{\ell}$, which means that inserting $\mathcal{F}(p)$ in a correlator has the effect of restricting the expansion to terms where the plaquette $p$ is occupied. Equivalently, $\mathcal{F}(p)$ is 1 if $p$ is occupied and 0 otherwise. Anticipating the next subsection, let us also write down the dual operators $\mathcal{E}_{\text{dual}}(\ell^*)$ and $\mathcal{F}_{\text{dual}}(p^*)$, where $\ell^*$ is a link of the dual lattice that is dual to (pierces) some plaquette $p$ of the original lattice, and $p^*$ is a plaquette of the dual lattice that is pierced by some link $\langle ij \rangle$ of the original lattice. By the definition of the dual membrane ensemble (see below),

$$ \mathcal{F}_{\text{dual}}(p^*) = \frac{1 - \sigma_{ij}\tau_i\tau_j}{2}, \quad (A4) $$

$$ \mathcal{E}_{\text{dual}}(\ell^*) = \frac{1 - \prod_{\ell \in p} \sigma_{\ell}}{2}. \quad (A5) $$

The duality relations given in Eqs. (7) and (8) of the main text follow from Eqs. (A2)–(A5) above.

The expansion above is the standard high-temperature expansion, meaning that terms are weighted by powers of the “fugacities” $x$ and $y$, which are small when $K$ and $J$ are small. Since the lattice is finite, the expansion may be done exactly, to all orders: i.e., one may think of it as a reformulation of the partition function and not as a perturbative series. It is a generalization of the high-temperature expansion of the Ising model, which would be obtained if the $\sigma$ field was absent and we just had loops associated with $y\tau\tau$.

## 2. Membranes on the dual lattice

Equation (A1) can be related to membranes on the dual lattice even more directly.

Let us choose the gauge $\tau = 1$ so that the partition function is a sum over only the $\sigma = \pm 1$ on each link. We can represent a given term by a collection of occupied links, where a link $\ell$ is occupied iff $\sigma_{\ell} = -1$. (Note that this notion of a link being occupied is unrelated to the one in the previous subsection.) Next, recall that plaquettes of the dual lattice are in 1:1 correspondence with links of the original lattice, so a configuration of occupied links is equivalent to a

configuration of occupied plaquettes $\tilde{\mathcal{M}}$ on the dual lattice. What is the Boltzmann weight of $\tilde{\mathcal{M}}$? Each occupied plaquette costs $x' \equiv (1-y)/(1+y)$ [from the ratio of the $1+y\sigma$ term in Eq. (A1) with $\sigma = -1$ to that with $\sigma = +1$]. Further, a link of the dual lattice where an odd number of occupied plaquettes meet means a square on the original lattice where $\prod \sigma = -1$. So each link in $\partial \tilde{\mathcal{M}}$ contributes $y' \equiv (1-x)/(1+x)$. Including the normalization,

$$ Z(x, y) = \frac{(1+x)^{3L^3} (1+y)^{3L^3}}{2^{3L^3}} \sum_{\tilde{\mathcal{M}}} x'^{|\tilde{\mathcal{M}}|} y'^{|\partial \tilde{\mathcal{M}}|}. \quad (A6) $$

## 3. Manifestly self-dual representation

Next, we demonstrate the reformulation in terms of two species of loops (or, more precisely, clusters), cf. Fig. 5.

In addition to the degrees of freedom $\sigma$ and $\tau$ on the links and sites (respectively) of the original lattice, let us add degrees of freedom $\tilde{\sigma}$ and $\tilde{\tau}$ on the links and sites of the dual lattice. Let us denote the links of the original lattice by $\mathcal{L}$ and those of the dual lattice by $\tilde{\mathcal{L}}$. Define

$$ Z' = \sum_{\sigma, \tau, \tilde{\sigma}, \tilde{\tau}} e^{-S_{\text{top}}[\sigma, \tilde{\sigma}]} \prod_{\ell \in \tilde{\mathcal{L}}} (1 + y\sigma\tau\tau) \prod_{\tilde{\ell} \in \tilde{\mathcal{L}}} (1 + y'\tilde{\sigma}\tilde{\tau}\tilde{\tau}). \quad (A7) $$

The “topological” action $S_{\text{top}}[\sigma, \tilde{\sigma}]$ is both gauge invariant and symmetric between $\sigma$ and $\tilde{\sigma}$: $e^{-S_{\text{top}}} = (-1)^{\hat{X}}$, where $\hat{X}$ is the $\mathbb{Z}_2$ linking number of the flux lines of $\sigma$ with those of $\tilde{\sigma}$. However, it is convenient here to define it as

$$ e^{-S_{\text{top}}} = \prod_{\tilde{\ell} \in \tilde{\mathcal{L}}} \left( \delta_{\tilde{\sigma}_{\tilde{\ell},1}} + \delta_{\tilde{\sigma}_{\tilde{\ell},-1}} \prod \sigma \right), \quad (A8) $$

where these properties are not manifest.

To see the equivalence to the original Ising gauge theory (A1), we simply pick the gauge $\tilde{\tau} = 1$ and do the sum on $\tilde{\sigma}$ separately for each link,

$$ \sum_{\tilde{\sigma}} \left( \delta_{\tilde{\sigma}_{\tilde{\ell},1}} + \delta_{\tilde{\sigma}_{\tilde{\ell},-1}} \prod \sigma \right) (1 + y'\tilde{\sigma}) \\ = (1 + y') \left( 1 + x \prod \sigma \right), \quad (A9) $$

so that

$$ Z' = 2^{L^3} (1 + y')^{3L^3} \sum_{\{\sigma\}, \{\tau\}} \prod_{\ell} (1 + x \prod \sigma) \prod_{\tilde{\ell}} (1 + y\sigma\tau\tau). \quad (A10) $$

To obtain the expression in terms of $C_e$ and $C_m$ in Eq. (5), we first perform the graphical expansion of the two products in Eq. (A7), giving the sum over loop
configurations $C_e$ and $C_m$ (Fig. 5). In addition to the fugacities $y$ and $y'$, these are weighted by

$$ \sum_{\sigma, \tilde{\sigma}} e^{-S_{\text{top}}} \left( \prod_{\ell \in C_e} \sigma_{\ell} \right) \left( \prod_{\tilde{\ell} \in C_m} \tilde{\sigma}_{\tilde{\ell}} \right) = 2^{(4L^3+2)} (-1)^{X(C_e, C_m)}. \quad (A11) $$

We can see this by using Eq. (A8) to make a graphical expansion of the left-hand side above in terms of a membrane configuration $M$ on the original lattice with boundary $\partial M = C_e$. For a given term in the expansion, the $\tilde{\sigma}$ are fixed by the Kronecker deltas, which dictates the sign of the product $\prod_{\tilde{\ell} \in C_m} \tilde{\sigma}_{\tilde{\ell}}$ on the rhs of Eq. (A11). There are, for periodic boundary conditions, $2^{L^3+2}$ choices of $M$ for fixed $\partial M = C_e$, but they all give the same sign. Altogether,

$$ Z' = 2^{(6L^3+2)} \sum_{C_e, C_m} (-1)^{X(C_e, C_m)} y^{|C_e|} y'^{|C_m|}. \quad (A12) $$

Summarizing, $Z'$ can be related to $Z$ in Eqs. (3) and (A1) by

$$ Z = 4c \sum_{C_e, C_m} (-1)^{X(C_e, C_m)} y^{|C_e|} y'^{|C_m|}, \quad (A13) $$

where $c = 2^{L^3}(1+y')^{-3L^3} = 2^{-2L^3}(1+x)^{3L^3}$, which gives Eq. (5) of the main text. The proportionality constant $c$ [omitted in Eq. (5)] is a trivial (nonuniversal) extensive contribution to the free energy, but the factor of 4 is universal and is the ground-state degeneracy of the 2D quantum system in its deconfined phase.

## APPENDIX B: FIRST-ORDER COEXISTENCE

Although we have concentrated our study on the vicinity of the multicritical point, we can extract from the data some information related to the first-order coexistence region along the self-dual line. Starting from the deconfined phase

<figure><img src="image_29.png" /><figcaption>FIG. 27. We show $b_4(A) = -(1/2)\kappa_4(A)/\text{Var}(A)^2$ as a function of $x$ for three system sizes. The lines are cubic polynomial fits.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIG. 28. Top panel: average plaquette occupation number versus $x$ obtained from $\langle |A| \rangle$ and $\langle S \rangle$ on the self-dual line (see text). The limits of these curves as $L \to \infty$ give the two equations of state for the two coexisting phases. The colors in the legend indicate different system sizes. The dashed black line is the average of the two coexisting phases, determined by $\langle A \rangle = 0$, and the black stars indicate the locations of $x_{\text{cep}}$ and $x_c$ ($x_{\text{cep}} < x_c$). Lines are just a guide to the eye. Bottom panel: similarly for the average link occupation number.</figcaption></figure>

(large $x$), this region starts at the multicritical point $x_c$ and ends at a critical endpoint $x_{\text{cep}}$. The estimate $x_c \approx 0.6367$ was obtained in Secs. IV and V. The location of $x_{\text{cep}}$ is, in principle, easier to determine because in this region, $b_4(A)$ [defined in Eq. (21)] behaves monotonically and presents a crossing, as shown in Fig. 27. From the figure, we roughly estimate $x_{\text{cep}} \approx 0.605$. Though this is a rough estimate, it is worth noting that the value of $b_4(A)$ at the crossing point is consistent with standard Ising universality (as for the liquid-gas critical endpoint), for which $b_4(A) \approx 0.7$ [189].

In between $x_c$ and $x_{\text{cep}}$, histograms of $A$ or of the total membrane area or membrane boundary length have two peaks, corresponding to the two coexisting phases. For large system sizes, our MC scheme will not properly sample both minima, so it could become hard to obtain equations of state for each phase. However, we can exploit the symmetry properties of $A$ and $S$. Denoting expectation values in the two equilibria by $\langle \dots \rangle_{\pm}$, in the thermodynamic
limit, we have $\langle A \rangle_{\pm} = \pm \langle |A| \rangle$. Therefore, by Eqs. (9) and (10),

$$ \frac{\langle \text{Area} \rangle_{\pm}}{3L^3} = \frac{\pm \langle |A| \rangle + \langle S \rangle}{12}, \qquad (B1) $$

$$ \frac{\langle \text{Length} \rangle_{\pm}}{3L^3} = \frac{1 - x^2 \pm \langle |A| \rangle - \langle S \rangle}{2x} + \frac{1 - x}{2}, \qquad (B2) $$

on the coexistence region of the self-dual line and in the thermodynamic limit. These equations give equations of state for each phase. The results are shown in Fig. 28.

## APPENDIX C: DETAILS OF MC SCHEME AND OF FITS

For most of our simulations, each MC step consists in updating all of the plaquettes (taking each of the three orientations in turn) and then updating all of the cubes. To allow parallelization, we divide plaquettes parallel to the (x, y) plane into two sublattices, and similarly for plaquettes in the (y, z) and (x, z) planes. We also divide cubes into two sublattices. We use one MC step as our unit of time. We study system sizes up to L = 96, and our longest simulations had $4 \times 10^9$ MC steps. Error bars for cumulants of A and S are calculated using bootstrap methods [37] (for this purpose, the correlation time is estimated as the time for the correlation to decay by a factor of 10).

For the fits in Sec. V, the scaling functions were described using B-splines with 8 to 12 degrees of freedom. The data used for the fits were restricted to $x \in [0.633, 0.640]$ and to scaling variable $z \in [-0.5, 0.5]$, although the particular intervals could slightly change from fit to fit. The system sizes included correspond to the best statistical fit, in the sense that the p-value (probability of getting a $\chi^2$ value

<figure><img src="image_31.png" /><figcaption>FIG. 29. Main panel: scaling collapse of $P_s$ as a function of $(y - y_{c,I})L^{1/\nu_I}$, for $x = 1$, with $y_{c,I} = 0.21809$ and $\nu_I = 0.63012$. Inset: $M_{\text{max}}$ as a function of the system size L at $x = 1$, $y = 0.218$. The straight line shows a power law using the fractal dimensions of Ising worldlines $d_{f,I} = 1.7349$.</figcaption></figure>

below the one obtained from the fit, for the degrees of freedom used) was maximized.

## APPENDIX D: FURTHER PERCOLATION DATA

### 1. Percolation at the Ising* transition ($x = 1$)

The transition at $x = 1$ maps to the 3D Ising critical point. Up to a difference in boundary conditions, the worldlines are those in a standard high-temperature expansion of the Ising model, and the percolation transition happens precisely at the Ising critical point for the cubic lattice, $y = 0.21809$. We indeed find that curves for the spanning probability $P_s$ cross close to this value and can be collapsed by plotting as a function of $(y - y_c)L^{1/\nu}$, using known Ising critical values, $y_{c,I} = 0.21809$ and $\nu_I = 0.63012$ [131] (see Fig. 29). We also check that the mass of the largest loop (number of links, $M_{\text{max}}$) follows a power law with a fractal dimension consistent with the known value $d_{f,I} = 1.7349(65)$ for Ising worldlines [131]. The inset of Fig. 29 shows $M_{\text{max}}$ as a function of the system size at $y = 0.218$.

### 2. Percolation on the $y=1$ boundary

When the percolation transition takes place within the thermodynamically trivial phase, we expect conventional percolation universality [190]. As an example, we consider the case $y = 1$. An attempt to obtain scaling collapse of $P_s$ suggests that finite-size effects are important for this range of system sizes. Figure 30 shows an attempt at scaling collapse using $\nu_P = 0.8762$ [192]. An estimate of the fractal dimension of the loops from $M_{\text{max}}$ (inset to Fig. 30) gives $d_f = 2.56$ (to be compared with 2.53 for the percolation universality class).

<figure><img src="image_32.png" /><figcaption>FIG. 30. Main panel: scaling collapse of $P_s$ as a function of $(x - x_{c,P})L^{1/\nu_P}$, for $y = 1$, with $x_{c,P} = 0.0865$ and $\nu_P = 0.8762$. Inset: $M_{\text{max}}$ as a function of the system size L at $x = 0.087$. The straight line shows the power law with fractal dimension of percolation universality class, $d_{f,P} = 2.53$.</figcaption></figure>
### 3. Percolation on the self-dual line

The phase diagram in Fig. 22 in the main text shows that we encounter several percolation transitions as we move along the self-dual line. We have shown data close to the self-dual critical point $x_c$ in the main text. For smaller $x$, we encounter the first-order line where two phases coexist, one with $A > 0$ and one with $A < 0$. To separate the properties of the two coexisting phases, we may average $P_s$ separately for configurations with $A > 0$ and $A < 0$. The phase with $A > 0$ appears to percolate throughout the entire range of the first-order line. Therefore, the phase with $A < 0$ must also percolate for some region of the first-order line close to the critical endpoint (since the two phases become identical there). One possibility (at first sight, the more natural) is that the phase with $A < 0$ undergoes a percolation transition at some intermediate $x$ lying on the interior of the first-order line. Another possibility is that this transition is pushed all the way to $x_c$, with the $A < 0$ phase having an extremely weak but nonzero percolation order parameter for $x \lesssim x_c$. Data for small sizes do not allow us to determine which of these occurs.

[1] M. E. Fisher, *The Renormalization Group in the Theory of Critical Behavior*, Rev. Mod. Phys. **46**, 597 (1974).

[2] S. Sachdev, *Quantum Phase Transitions, Handbook of Magnetism and Advanced Magnetic Materials* (Cambridge University Press, 2011).

[3] F. J. Wegner, *Duality in Generalized Ising Models and Phase Transitions without Local Order Parameters*, J. Math. Phys. (N.Y.) **12**, 2259 (1971).

[4] J. M. Kosterlitz and D. J. Thouless, *Ordering, Metastability and Phase Transitions in Two-Dimensional Systems*, J. Phys. C **6**, 1181 (1973).

[5] P.-G. de Gennes, *Exponents for the Excluded Volume Problem as Derived by the Wilson Method*, Phys. Lett. A **38**, 339 (1972).

[6] D. A. Huse and S. Leibler, *Are Sponge Phases of Membranes Experimental Gauge-Higgs Systems?*, Phys. Rev. Lett. **66**, 437 (1991).

[7] P. E. Lammert, D. S. Rokhsar, and J. Toner, *Topology and Nematic Ordering*, Phys. Rev. Lett. **70**, 1650 (1993).

[8] P. Francesco, P. Mathieu, and D. Sénéchal, *Conformal Field Theory* (Springer Science & Business Media, New York, 2012).

[9] E. Fradkin, *Field Theories of Condensed Matter Physics* (Cambridge University Press, Cambridge, England, 2013).

[10] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, *Deconfined Quantum Critical Points*, Science **303**, 1490 (2004).

[11] E. Fradkin and S. H. Shenker, *Phase Diagrams of Lattice Gauge Theories with Higgs Fields*, Phys. Rev. D **19**, 3682 (1979).

[12] G. A. Jongeward, J. D. Stack, and C. Jayaprakash, *Monte Carlo Calculations on $Z_2$ Gauge-Higgs Theories*, Phys. Rev. D **21**, 3360 (1980).

[13] I. S. Tupitsyn, A. Kitaev, N. V. Prokof'ev, and P. C. E. Stamp, *Topological Multicritical Point in the Phase Diagram of the Toric Code Model and Three-Dimensional Lattice Gauge Higgs Model*, Phys. Rev. B **82**, 085114 (2010).

[14] J. Vidal, S. Dusuel, and K. P. Schmidt, *Low-Energy Effective Theory of the Toric Code Model in a Parallel Magnetic Field*, Phys. Rev. B **79**, 033109 (2009).

[15] J. B. Kogut, *An Introduction to Lattice Gauge Theory and Spin Systems*, Rev. Mod. Phys. **51**, 659 (1979).

[16] N. Read and B. Chakraborty, *Statistics of the Excitations of the Resonating-Valence-Bond State*, Phys. Rev. B **40**, 7133 (1989).

[17] S. Kivelson, *Statistics of Holons in the Quantum Hard-Core Dimer Gas*, Phys. Rev. B **39**, 259 (1989).

[18] N. Read and S. Sachdev, *Large-N Expansion for Frustrated Quantum Antiferromagnets*, Phys. Rev. Lett. **66**, 1773 (1991).

[19] X. G. Wen, *Mean-Field Theory of Spin-Liquid States with Finite Energy Gap and Topological Orders*, Phys. Rev. B **44**, 2664 (1991).

[20] T. Senthil and M. P. A. Fisher, $Z_2$ Gauge Theory of Electron Fractionalization in Strongly Correlated Systems, Phys. Rev. B **62**, 7850 (2000).

[21] R. Moessner, S. L. Sondhi, and E. Fradkin, *Short-Ranged Resonating Valence Bond Physics, Quantum Dimer Models, and Ising Gauge Theories*, Phys. Rev. B **65**, 024504 (2001).

[22] A. Y. Kitaev, *Fault-Tolerant Quantum Computation by Anyons*, Ann. Phys. (Amsterdam) **303**, 2 (2003).

[23] R. Balian, J. M. Drouffe, and C. Itzykson, *Gauge Fields on a Lattice. I. General Outlook*, Phys. Rev. D **10**, 3376 (1974).

[24] F. Wu, Y. Deng, and N. Prokof'ev, *Phase Diagram of the Toric Code Model in a Parallel Magnetic Field*, Phys. Rev. B **85**, 195104 (2012).

[25] D. A. Huse and S. Leibler, *Phase Behaviour of an Ensemble of Nonintersecting Random Fluid Films*, J. Phys. (Paris) **49**, 605 (1988).

[26] M. E. Cates, D. Roux, D. Andelman, S. T. Milner, and S. A. Safran, *Random Surface Model for the $L_3$-Phase of Dilute Surfactant Solutions*, Europhys. Lett. **5**, 733 (1988).

[27] F. David, *O(n) Gauge Models and Self-Avoiding Random Surfaces in Three Dimensions*, Europhys. Lett. **9**, 575 (1989).

[28] D. Roux, M. E. Cates, U. Olsson, R. C. Ball, F. Nallet, and A. M. Bellocq, *Light Scattering from a Surfactant Sponge Phase: Evidence for a Hidden Symmetry*, Europhys. Lett. **11**, 229 (1990).

[29] D. Roux, C. Coulon, and M. E. Cates, *Sponge Phases in Surfactant Solutions*, J. Phys. Chem. **96**, 4174 (1992).

[30] D. Roux, *Sponge Phases: An Example of Random Surfaces*, Physica (Amsterdam) **213A**, 168 (1995).

[31] L. Peliti, *Amphiphilic Membranes*, in Fluctuating Geometries in Statistical Mechanics and Field Theory, edited by F. David, P. Ginsparg, and J. Zinn-Justin (Elsevier, Amsterdam, 1996), pp. 195-285.

[32] R. Nandkishore, M. A. Metlitski, and T. Senthil, *Orthogonal Metals: The Simplest Non-Fermi Liquids*, Phys. Rev. B **86**, 045128 (2012).
[33] Z. Nussinov and G. Ortiz, *A Symmetry Principle for Topological Quantum Order*, Ann. Phys. (Amsterdam) **324**, 977 (2009).

[34] A. Kapustin and R. Thorngren, *Higher Symmetry and Gapped Phases of Gauge Theories*, in *Algebra, Geometry, and Physics in the 21st Century* (Springer, New York, 2017), pp. 177–202.

[35] D. Gaiotto, A. Kapustin, N. Seiberg, and B. Willett, *Generalized Global Symmetries*, J. High Energy Phys. **02** (2015) 172.

[36] X.-G. Wen, *Emergent Anomalous Higher Symmetries from Topological Order and from Dynamical Electromagnetic Field in Condensed Matter Systems*, Phys. Rev. B **99**, 205139 (2019).

[37] M. E. J. Newman and G. T. Barkema, *Monte Carlo Methods in Statistical Physics* (Clarendon Press, Oxford, 1999).

[38] S. D. Geraedts and O. I. Motrunich, *Monte Carlo Study of a $U(1) \times U(1)$ System with $\pi$-Statistical Interaction*, Phys. Rev. B **85**, 045114 (2012).

[39] F. J. Burnell, *Anyon Condensation and Its Applications*, Annu. Rev. Condens. Matter Phys. **9**, 307 (2018).

[40] A. Nanda, K. Dhochak, and S. Bhattacharjee, *Phases and Quantum Phase Transitions in an Anisotropic Ferromagnetic Kitaev-Heisenberg-$\gamma$ Magnet*, Phys. Rev. B **102**, 235124 (2020).

[41] S. D. Geraedts and O. I. Motrunich, *Phases and Phase Transitions in a $U(1) \times U(1)$ System with $\theta = 2\pi/3$ Mutual Statistics*, Phys. Rev. B **86**, 045106 (2012).

[42] J. Y. Lee, S. Geraedts, and O. I. Motrunich, *Monte Carlo Study of Phase Transitions out of Symmetry-Enriched Topological Phases of Bosons in Two Dimensions*, Phys. Rev. B **93**, 035103 (2016).

[43] S. D. Geraedts and O. I. Motrunich, *Line of Continuous Phase Transitions in a Three-Dimensional $U(1)$ loop model with $1/r^2$ Current-Current Interactions*, Phys. Rev. B **85**, 144303 (2012).

[44] M. Levin and T. Senthil, *Deconfined Quantum Criticality and Néel Order via Dimer Disorder*, Phys. Rev. B **70**, 220403(R) (2004).

[45] O. I. Motrunich and A. Vishwanath, *Emergent Photons and Transitions in the $O(3)$ Sigma Model with Hedgehog Suppression*, Phys. Rev. B **70**, 075104 (2004).

[46] T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, *Deconfined Criticality Critically Defined*, J. Phys. Soc. Jpn. **74**, 1 (2005).

[47] A. Tanaka and X. Hu, *Many-Body Spin Berry Phases Emerging from the $\pi$-Flux State: Competition between Antiferromagnetism and the Valence-Bond-Solid State*, Phys. Rev. Lett. **95**, 036402 (2005).

[48] T. Senthil and M. P. A. Fisher, *Competing Orders, Non-linear Sigma Models, and Topological Terms in Quantum Magnets*, Phys. Rev. B **74**, 064405 (2006).

[49] A. W. Sandvik, *Evidence for Deconfined Quantum Criticality in a Two-Dimensional Heisenberg Model with Four-Spin Interactions*, Phys. Rev. Lett. **98**, 227202 (2007).

[50] S. Bhattacharjee, *Quantum Destruction of Spiral Order in Two-Dimensional Frustrated Magnets*, Phys. Rev. B **84**, 104430 (2011).

[51] C. Wang, A. Nahum, M. A. Metlitski, C. Xu, and T. Senthil, *Deconfined Quantum Critical Points: Symmetries and Dualities*, Phys. Rev. X **7**, 031051 (2017).

[52] C.-M. Jian, A. Rasmussen, Y.-Z. You, and C. Xu, *Emergent Symmetry and Tricritical Points near the Deconfined Quantum Critical Point*, arXiv:1708.03050.

[53] S. Gazit, F. F. Assaad, S. Sachdev, A. Vishwanath, and C. Wang, *Confinement Transition of 2 Gauge Theories Coupled to Massless Fermions: Emergent Quantum Chromodynamics and SO(5) Symmetry*, Proc. Natl. Acad. Sci. U.S.A. **115**, E6987 (2018).

[54] S. Jiang and O. Motrunich, *Ising Ferromagnet to Valence Bond Solid Transition in a One-Dimensional Spin Chain: Analogies to Deconfined Quantum Critical Points*, Phys. Rev. B **99**, 075103 (2019).

[55] B. Zhao, J. Takahashi, and A. W. Sandvik, *Multicritical Deconfined Quantum Criticality and Lifshitz Point of a Helical Valence-Bond Phase*, Phys. Rev. Lett. **125**, 257204 (2020).

[56] M. Freedman, C. Nayak, K. Shtengel, K. Walker, and Z. Wang, *A Class of P, T-Invariant Topological Phases of Interacting Electrons*, Ann. Phys. (Amsterdam) **310**, 428 (2004).

[57] M. Freedman, C. Nayak, and K. Shtengel, *Line of Critical Points in 2+1 Dimensions: Quantum Critical Loop Gases and Non-Abelian Gauge Theory*, Phys. Rev. Lett. **94**, 147205 (2005).

[58] M. Freedman, C. Nayak, and K. Shtengel, *Lieb-Schultz-Mattis Theorem for Quasitopological Systems*, Phys. Rev. B **78**, 174411 (2008).

[59] M. Troyer, S. Trebst, K. Shtengel, and C. Nayak, *Local Interactions and Non-Abelian Quantum Loop Gases*, Phys. Rev. Lett. **101**, 230401 (2008).

[60] Z. Dai and A. Nahum, *Quantum Criticality of Loops with Topologically Constrained Dynamics*, Phys. Rev. Research **2**, 033051 (2020).

[61] D. R. Nelson, T. Piran, and S. Weinberg, *Statistical Mechanics of Membranes and Surfaces* (World Scientific, Singapore, 2004).

[62] J. Polchinski, *String Theory: An Introduction to the Bosonic String* (Cambridge University Press, Cambridge, England, 1998), Vol. 1.

[63] V. I. Dotsenko and A. M. Polyakov, *Fermion Representations for the 2D and 3D Ising Models*, in *Conformal Field Theory and Solvable Lattice Models* (Mathematical Society of Japan, 1988), pp. 171–203.

[64] N. Iqbal and J. McGreevy, *Toward a 3d Ising Model with a Weakly-Coupled String Theory Dual*, SciPost Phys. **9**, 019 (2020).

[65] We have used a $\propto$ sign rather than an equality because we choose to absorb a trivial constant into the definition of Z.

[66] K. Gregor, D. A. Huse, R. Moessner, and S. L. Sondhi, *Diagnosing Deconfinement and Topological Order*, New J. Phys. **13**, 025009 (2011).

[67] Explicitly, $Z(x,y)=c^{3L^3}Z(x',y')$, with $c=[(1+x)(1+y)/2]$.

[68] The configuration $C_e$ obeys the same restrictions as $\partial M$: It must make sense to regard it as the boundary for a membrane configuration on the original lattice. The same holds for $C_m$ on the dual lattice. Specifically, each site on the original lattice lies on an even number of the links in $C_e$.
(possibly zero) and similarly for the dual lattice and $C_m$. Additionally, the $Z_2$ winding numbers of $C_e, C_m$ in each of the three directions are equal to zero modulo 2 (this is a requirement for a worldline configuration to be promoted to the boundary of a membrane configuration).

[69] C. Wang, A. Nahum, and T. Senthil, *Topological Para-magnetism in Frustrated Spin-1 Mott Insulators*, Phys. Rev. B **91**, 195131 (2015).

[70] S. Geraedts and O. Motrunich (unpublished).

[71] C. W. von Keyserlingk, F. J. Burnell, and S. H. Simon, *Three-Dimensional Topological Lattice Models with Surface Anyons*, Phys. Rev. B **87**, 045107 (2013).

[72] If a basis transformation is applied to the Hamiltonian, the symmetry becomes a simple translation [40,73].

[73] X.-G. Wen, *Quantum Orders in an Exact Soluble Model*, Phys. Rev. Lett. **90**, 016803 (2003).

[74] There is an order-1 degeneracy that depends on the system's topology.

[75] E. Cobanera, G. Ortiz, and Z. Nussinov, *The Bond-Algebraic Approach to Dualities*, Adv. Phys. **60**, 679 (2011).

[76] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, *Non-Abelian Anyons and Topological Quantum Computation*, Rev. Mod. Phys. **80**, 1083 (2008).

[77] S. Dusuel, M. Kamfor, R. Orús, K. P. Schmidt, and J. Vidal, *Robustness of a Perturbed Topological Phase*, Phys. Rev. Lett. **106**, 107203 (2011).

[78] A. Nahum, P. Serna, and A. M. Somoza (to be published).

[79] In this Ising model, we must sum over both periodic and antiperiodic boundary conditions for each direction to allow an odd number of domain walls to span the system perpendicular to each direction.

[80] See Ref. [81] for a careful discussion of the role of translation symmetries in the continuum in some other spin models.

[81] M. A. Metlitski and R. Thorngren, *Intrinsic and Emergent Anomalies at Deconfined Critical Points*, Phys. Rev. B **98**, 085140 (2018).

[82] C. Heinrich, F. Burnell, L. Fidkowski, and M. Levin, *Symmetry-Enriched String Nets: Exactly Solvable Models for Set Phases*, Phys. Rev. B **94**, 235136 (2016).

[83] M. Cheng, Z.-C. Gu, S. Jiang, and Y. Qi, *Exactly solvable models for symmetry-enriched topological phases*, Phys. Rev. B **96**, 115107 (2017).

[84] The existence of microscopic models for the deconfined phase with an “on-site” self-duality symmetry [82,83] suggests that this symmetry of the critical theory is not anomalous. In this respect, it is different from the Kramers-Wannier duality symmetry of the 1+1D Ising model, which cannot be realized as an on-site symmetry [85–87].

[85] D. Aasen, R. S. K. Mong, and P. Fendley, *Topological Defects on the Lattice: I. The Ising Model*, J. Phys. A **49**, 354001 (2016).

[86] A. Karch, D. Tong, and C. Turner, *A Web of 2D Dualities: $Z_2$ Gauge Fields and Arf Invariants*, SciPost Phys. **7**, 7 (2019).

[87] R. A. Jones and M. A. Metlitski, *1D Lattice Models for the Boundary of 2D “Majorana” Fermion SPTs: Kramers-Wannier Duality as an Exact $Z_2$ Symmetry*, arXiv:1902.05957.

[88] K. Binder, *Critical Properties from Monte Carlo Coarse Graining and Renormalization*, Phys. Rev. Lett. **47**, 693 (1981).

[89] The standard definition of the Binder cumulant is $U_L = 2b_4/3$.

[90] To see that Eq. (25) holds, recall that $L^3\text{Var}(S)$ is given by the integral over the connected two-point function of $S(r)$. If we neglect terms of size $(x_S - 3/2)\log L$, then it is sufficient to replace the power law $r^{-2x_S}$ occurring in this integral with $r^{-3}$:

$$L^3 \text{Var}(S) \approx \int_{[0,L]^3} d^3 r r^{-3} H\left(\frac{r}{L}, z\right) + B. \quad (\text{D1})$$

Here, $H$ is a scaling function for the correlator in the finite system, and $H(0, 0) = \alpha_S^2$. The integral is cut off at $r$ of order 1, and the constant $B$ represents a short-distance contribution. This integral gives Eq. (25), in which the nonuniversal constant $B$ has been absorbed into the function $f(z)$.

[91] J. Cardy, *Scaling and Renormalization in Statistical Physics* (Cambridge University Press, Cambridge, England, 1996), Vol. 5.

[92] M. Caselle, G. Costagliola, and N. Magnoli, *Numerical Determination of the Operator-Product-Expansion Coefficients in the 3D Ising Model from Off-Critical Correlators*, Phys. Rev. D **91**, 061901(R) (2015).

[93] M. Hasenbusch, *Two- and Three-Point Functions at Criticality: Monte Carlo Simulations of the Three-Dimensional ($q+1$)-State Clock Model*, Phys. Rev. B **102**, 224509 (2020).

[94] The dynamical exponent $z$ of the 3+1D stochastic dynamics should not be confused with the dynamical exponent $z_{QM} = 1$ of the 2+1D quantum system.

[95] R. J. Glauber, *Time-Dependent Statistics of the Ising Model*, J. Math. Phys. (N.Y.) **4**, 294 (1963).

[96] P. C. Hohenberg and B. I. Halperin, *Theory of Dynamic Critical Phenomena*, Rev. Mod. Phys. **49**, 435 (1977).

[97] S. Wansleben and D. Landau, *Dynamical Critical Exponent of the 3D Ising Model*, J. Appl. Phys. **61**, 3968 (1987).

[98] S. Wansleben and D. P. Landau, *Monte Carlo Investigation of Critical Dynamics in the Three-Dimensional Ising Model*, Phys. Rev. B **43**, 6006 (1991).

[99] C. Münkel, D. W. Heermann, J. Adler, M. Gofman, and D. Stauffer, *The Dynamical Critical Exponent of the Two-, Three- and Five-Dimensional Kinetic Ising Model*, Physica (Amsterdam) **193A**, 540 (1993).

[100] N. Ito, *Non-Equilibrium Critical Relaxation of the Three-Dimensional Ising Model*, Physica (Amsterdam) **192A**, 604 (1993).

[101] P. Grassberger, *Damage Spreading and Critical Exponents for “Model A” Ising Dynamics*, Physica (Amsterdam) **214A**, 547 (1995).

[102] A. Jaster, J. Mainville, L. Schülke, and B. Zheng, *Short-Time Critical Dynamics of the Three-Dimensional Ising Model*, J. Phys. A **32**, 1395 (1999).

[103] N. Ito, K. Hukushima, K. Ogawa, and Y. Ozeki, *Non-equilibrium Relaxation of Fluctuations of Physical Quantities*, J. Phys. Soc. Jpn. **69**, 1931 (2000).

[104] Y. Murase and N. Ito, *Dynamic Critical Exponents of Three-Dimensional Ising Models and Two-Dimensional Ising Model*, Phys. Rev. Lett. **106**, 031601 (2011).
*Three-States Potts Models*, J. Phys. Soc. Jpn. **77**, 014002 (2008).

[105] M. Collura, *Off-Equilibrium Relaxational Dynamics with an Improved Ising Hamiltonian*, J. Stat. Mech. (2010) P12036.

[106] D. Niermann, C. P. Grams, P. Becker, L. Bohatý, H. Schenck, and J. Hemberger, *Critical Slowing Down near the Multiferroic Phase Transition in MnWO₄*, Phys. Rev. Lett. **114**, 037204 (2015).

[107] D. Mesterhazy, J. H. Stockemer, and Y. Tanizaki, *From Quantum to Classical Dynamics: The Relativistic O(N) Model in the Framework of the Real-Time Functional Renormalization Group*, Phys. Rev. D **92**, 076001 (2015).

[108] C. Duclut and B. Delamotte, *Frequency Regulators for the Nonperturbative Renormalization Group: A General Study and the Model A as a Benchmark*, Phys. Rev. E **95**, 012107 (2017).

[109] L. T. Adzhemyan, E. V. Ivanova, M. V. Kompaniets, and S. Y. Vorobyeva, *Diagram Reduction in Problem of Critical Dynamics of Ferromagnets: 4-Loop Approximation*, J. Phys. A **51**, 155003 (2018).

[110] M. Hasenbusch, *Dynamic Critical Exponent z of the Three-Dimensional Ising Universality Class: Monte Carlo Simulations of the Improved Blume-Capel Model*, Phys. Rev. E **101**, 022126 (2020).

[111] As discussed in Sec. III A, the lattice operators are only self-dual or anti-self-dual up to derivative terms.

[112] Detailed balance implies that the correlator is invariant under $t \to -t$, so we take $t > 0$.

[113] At first glance, we might expect the relaxation times for duality-odd and duality-even operators to differ by a nontrivial order-one factor, due to coupling to eigenstates of the Markov matrix with different symmetry under duality. But the Monte Carlo dynamics is not, in fact, symmetric under duality, so all operators can couple to the lowest excited state, whose eigenvalue determines $\tau$.

[114] N. Xu, C. Castelnovo, R. G. Melko, C. Chamon, and A. W. Sandvik, *Dynamic Scaling of Topological Ordering in Classical Systems*, Phys. Rev. B **97**, 024432 (2018).

[115] Analogs of this phenomenon may also be found in 2D [116].

[116] Y. Shi, A. Lamacraft, and P. Fendley, *Boson Pairing and Unusual Criticality in a Generalized XY Model*, Phys. Rev. Lett. **107**, 240601 (2011).

[117] If the original system is finite with periodic boundary conditions, then we must sum over periodic and antiperiodic boundary conditions for $\phi_r$.

[118] C. D. Batista and Z. Nussinov, *Generalized Elitzur's Theorem and Dimensional Reductions*, Phys. Rev. B **72**, 045137 (2005).

[119] Z. Nussinov and G. Ortiz, *Sufficient Symmetry Conditions for Topological Quantum Order*, Proc. Natl. Acad. Sci. U.S.A. **106**, 16944 (2009).

[120] Formally, the relation between the models with $Z_2$ global and $Z_2$ 1-form symmetries is via gauging of these symmetries with flat gauge fields [35].

[121] M. B. Hastings and X. G. Wen, *Quasiadiabatic Continuation of Quantum States: The Stability of Topological Ground-State Degeneracy and Emergent Gauge Invariance*, Phys. Rev. B **72**, 045141 (2005).

[122] Reference [121] gives a rigorous result for the case when all excitations are gapped.

[123] The various string operators can be connected to the Fredenhagen Marcu order parameter [124] and the Huse-Liebler horseshoe [6,66], long used as diagnostics for deconfinement. The dressed string operators above obey the property of invariance under deformations. “Bare” string or Wilson line operators do not, and their correlators generically include a product of local contributions from along the length of the string. However, universal data can still be extracted by dividing an appropriate open expectation value for an open Wilson line by the expectation value of a closed Wilson line in order to cancel the UV contributions [6,124].

[124] K. Fredenhagen and M. Marcu, *Confinement Criterion for QCD with Dynamical Quarks*, Phys. Rev. Lett. **56**, 223 (1986).

[125] We defer a discussion of numerical implementation to Ref. [78]. To have a simple convention in mind, we can define the surface associated with a given connected cluster of links in $\partial M$ as the (possibly self-intersecting) surface traced out when the cluster is shrunk down onto its center of mass.

[126] This is only a sufficient condition; see, e.g., the comment on doubled strands in Sec. IX B below.

[127] This is because there are, in principle, two ways for the loops $\partial M$ to undergo a percolation transition. Heuristically, these correspond to proliferation either of single strands or of doubled strands (thin ribbons). The former results in a Higgs transition, but the latter does not. This can be made slightly more precise by extending the lattice gauge-Higgs theory to include a replica index $\alpha = 1, ..., n$ on the matter field $\tau$, with the limit $n \to 1$ recovering the initial problem. Condensation of $\tau_\alpha$ results both in percolation and in a Higgs transition, while condensation of the composite field $\tau_\alpha \tau_\beta$, which is gauge neutral but charged under replica symmetry, results in percolation but not a Higgs transition [128] (in this case, the percolation transition has no thermodynamic effect [91]). The “double-line” percolation transition does not necessarily obstruct the membrane-patching procedure since there is no large-scale ambiguity about patching the surface between two nearby lines.

[128] A. Nahum and J. Chalker, *Universal Statistics of Vortex Lines*, Phys. Rev. E **85**, 031141 (2012).

[129] The probability that two links separated by a distance $r \gg 1$ lie on the same cluster decays as $r^{-2x_{\text{conn}}}$.

[130] A loop is defined to span the system in (say) the x direction if it visits each of the L distinct planes of x-directed links.

[131] F. Winter, W. Janke, and A. M. J. Schakel, *Geometric Properties of the Three-Dimensional Ising and XY Models*, Phys. Rev. E **77**, 061108 (2008).

[132] M. Kompaniets and K. J. Wiese, *Fractal Dimension of Critical Curves in the O(n)-Symmetric $\phi^4$ Model and Crossover Exponent at 6-Loop Order: Loop-Erased Random Walks, Self-Avoiding Walks, Ising, XY, and Heisenberg Models*, Phys. Rev. E **101**, 012104 (2020).

[133] Recall that the loops in Eq. (5) are really clusters since any even number of occupied links can meet at a node.
[134] To see this, let $i, j$ be indices running over the distinct loops in $C$. Let $s_i$ be a species index, with $s_i = 1$ for an $e$ worldline and $s_i = -1$ for an $m$ worldline. Finally, let $n_{i,j} = 0, 1$ be the $\mathbb{Z}_2$ linking number of loops $i$ and $j$, which is straightforwardly defined since all loops contract to a point. The linking sign for $C_e$ and $C_m$ may be written $(-1)^{\tilde{X}(C_e, C_m)} = e^{(i\pi/2) \sum_{i<j} n_{i,j}(2-s_i-s_j)}$. Summing over the $s_i$ gives the result in the text. The sum vanishes unless $\sum_{j(\neq i)} n_{i,j}$ is even for every $i$.

[135] B. Nienhuis, *Exact Critical Point and Critical Exponents of $O(n)$ Models in Two Dimensions*, Phys. Rev. Lett. **49**, 1062 (1982).

[136] In the model (37), as it stands, the smallest loop that can be nontrivially linked by another loop is a square of side length 2.

[137] We can see this by summing over all possible spatial patterns of insertions of $\tilde{A}$ (obtained by expanding in this perturbation) for a given configuration of loops. Since a large loop may lie on either an even or an odd number of $\tilde{A}$ insertions, its linking number may be either even or odd. Further, summing over patterns of insertions on an asymptotically large loop, with a fixed parity for the number of insertions, yields a factor of $1/2$ that cancels the fugacity 2 in Eq. (39).

[138] B. Kang, S. A. Parameswaran, A. C. Potter, R. Vasseur, and S. Gazit, *Superuniversality from Disorder at Two-Dimensional Topological Phase Transitions*, Phys. Rev. B **102**, 224204 (2020).

[139] In the classical interpretation of the critical point (where disorder is uncorrelated in all three directions, rather than being translationally invariant in the imaginary time direction) disorder that breaks self-duality is relevant (even if it preserves self-duality on average), while since $x_S \approx 1.5$, disorder that preserves self-duality is close to being marginal [91].

[140] B. Svetitsky, *Symmetry Aspects of Finite-Temperature Confinement Transitions*, Phys. Rep. **132**, 1 (1986).

[141] Such configurations are suppressed by a factor of about $\exp(const/T)$ because of the line tension of the worldlines. This factor can be mapped to an exponentially weak magnetic field $h \sim \exp(const/T)$ in an effective 2D Ising model [140].

[142] L. Genovese, F. Gliozzi, A. Rago, and C. Torrero, *The Phase Diagram of the Three-Dimensional $Z_2$ Gauge Higgs System at Zero and Finite Temperature*, Nucl. Phys. B Proc. Suppl. **119**, 894 (2003).

[143] G. Y. Zhu and G. M. Zhang, *Gapless Coulomb State Emerging from a Self-Dual Topological Tensor-Network State*, Phys. Rev. Lett. **122**, 176401 (2019).

[144] S. Pujari, F. Alet, and K. Damle, *Transitions to Valence-Bond Solid Order in a Honeycomb Lattice Antiferromagnet*, Phys. Rev. B **91**, 104411 (2015).

[145] H. Shao, W. Guo, and A. W. Sandvik, *Monte Carlo Renormalization Flows in the Space of Relevant and Irrelevant Operators: Application to Three-Dimensional Clock Models*, Phys. Rev. Lett. **124**, 080602 (2020).

[146] S. M. Chester, W. Landry, J. Liu, D. Poland, D. Simmons-Duffin, N. Su, and A. Vichi, *Carving out OPE Space and Precise $O(2)$ Model Critical Exponents*, J. High Energy Phys. 06 (2020) 142.

[147] M. Hasenbusch and E. Vicari, *Anisotropic Perturbations in Three-Dimensional $O(N)$-Symmetric Vector Models*, Phys. Rev. B **84**, 125136 (2011).

[148] M. Hasenbusch, *Monte Carlo Study of an Improved Clock Model in Three Dimensions*, Phys. Rev. B **100**, 224517 (2019).

[149] Our result for the fractal dimension of an $e$ worldline is also consistent with the value for an $XY$ worldline, but it is unclear at present whether this is an independent exponent or whether $d_f = 3 - x_A$ [see the discussion below Eq. (36) and in Sec. X A]. In the $XY$ model, $d_f^{XY} = 3 - x_{q_2^2-q_3^2}$ by virtue of known scaling relations. Our data for the OPE coefficient $C_{AAS}$ in Fig. 19 are not sufficient to make a very useful comparison, but the known $XY$ OPE coefficient value $C_{q_2^2-q_3^2,q_2^2-q_3^2,\phi^2} = 1.25213(14)$ [88,139] cannot be ruled out.

[150] Heuristically, we may imagine a real-space RG transformation in the self-dual loop ensemble of Sec. II C that renormalizes the local interactions between large loops (worldlines) by integrating out the small ones. The minus sign associated with linking is the basic reason why this RG transformation is unlike that in a model with two mutually local particle species $\phi_x$ and $\phi_y$. A possible rationalization of the closeness of the present exponents to $XY$ exponents could be that the effect of this linking sign on the RG transformation of the local interactions is numerically small. Loosely speaking, the picture would be that linking events between “large” loops occur with a weight that is nonzero and universal but that happens to be small. This is motivated by the observation that at the self-dual critical point, the universal loop percolation probability is numerically small (Fig. 23).

[151] T. H. Hansson, V. Oganesyan, and S. L. Sondhi, *Superconductors are Topologically Ordered*, Ann. Phys. (Amsterdam) **313**, 497 (2004).

[152] M. C. Diamantini, P. Sodano, and C. A. Trugenberger, *Topological Order in Frustrated Josephson Junction Arrays*, Europhys. Lett. **83**, 21003 (2008).

[153] S.-P. Kou, M. Levin, and X.-G. Wen, *Mutual Chern-Simons Theory for $Z_2$ Topological Order*, Phys. Rev. B **78**, 155134 (2008).

[154] J. McGreevy, *TASI 2015 Lectures on Quantum Matter (with a View Toward Holographic Duality)*, in *New Frontiers in Fields and Strings (TASI 2015)-Proceedings of the 2015 Theoretical Advanced Study Institute in Elementary Particle Physics* (World Scientific, Singapore, 2016), p. 215.

[155] K. G. Wilson, *Confinement of Quarks*, Phys. Rev. D **10**, 2445 (1974).

[156] M. Billó, M. Caselle, D. Gaiotto, F. Gliozzi, M. Meineri, and R. Pellegrini, *Line Defects in the 3D Ising Model*, J. High Energy Phys. 07 (2013) 055.

[157] S. Edwards, *The Theory of Rubber Elasticity*, Br. Polym. J. **9**, 140 (1977).

[158] P. G. De Gennes, *Scaling Concepts in Polymer Physics* (Cornell University Press, Ithaca, 1979).

[159] J. des Cloizeaux, *Ring Polymers in Solution: Topological Effects*, J. Phys. (Paris), Lett. **42**, 433 (1981).

[160] A. R. Khokhlov and S. K. Nechaev, *Polymer Chain in an Array of Obstacles*, Phys. Lett. A **112**, 156 (1985).
[161] M. Rubinstein, *Dynamics of Ring Polymers in the Presence of Fixed Obstacles*, Phys. Rev. Lett. **57**, 3023 (1986).

[162] M. E. Cates and J. M. Deutsch, *Conjectures on the statistics of ring polymers*, J. Phys. (Paris) **47**, 2121 (1986).

[163] J. D. Halverson, J. Smrek, K. Kremer, and A. Y. Grosberg, *From a Melt of Rings to Chromosome Territories: The Role of Topological Constraints in Genome Folding*, Rep. Prog. Phys. **77**, 022601 (2014).

[164] M. V. Imakaev, K. M. Tchourine, S. K. Nechaev, and L. A. Mirny, *Effects of Topological Constraints on Globular Polymers*, Soft Matter **11**, 665 (2015).

[165] P. Serna, G. Bunin, and A. Nahum, *Topological Constraints in Directed Polymer Melts*, Phys. Rev. Lett. **115**, 228303 (2015).

[166] P. De Forcrand, *Progress on Lattice QCD Algorithms*, arXiv:hep-lat/9509082.

[167] M. Peardon, *Progress in Lattice Algorithms*, arXiv:hep-lat/0201003.

[168] S. Schaefer, *Algorithms for Lattice QCD: Progress and Challenges*, in *AIP Conference Proceedings* (American Institute of Physics, New York, 2011), Vol. 1343, pp. 93–98, https://doi.org/10.1063/1.3574948.

[169] K. I. Ishikawa, *Recent Algorithm and Machine Developments for Lattice QCD*, Proc. Sci. LATTICE 2008 (2009) 013.

[170] J. Liu, Y. Qi, Z. Y. Meng, and L. Fu, *Self-Learning Monte Carlo Method*, Phys. Rev. B **95**, 041101 (2017).

[171] L. Huang and L. Wang, *Accelerated Monte Carlo Simulations with Restricted Boltzmann Machines*, Phys. Rev. B **95**, 035105 (2017).

[172] H. Suwa, A. Sen, and A. W. Sandvik, *Level Spectroscopy in a Two-Dimensional Quantum Magnet: Linearly Distinguishing Spinons at the Deconfined Quantum Critical Point*, Phys. Rev. B **94**, 144416 (2016).

[173] L. Vanderstraeten, M. Mariën, F. Verstraete, and J. Haegelein, *Excitations and the Tangent Space of Projected Entangled-Pair States*, Phys. Rev. B **92**, 201111(R) (2015).

[174] L. Vanderstraeten, J. Haegeman, and F. Verstraete, *Simulating Excitation Spectra with Projected Entangled-Pair States*, Phys. Rev. B **99**, 165121 (2019).

[175] M. Dalmonte and S. Montangero, *Lattice Gauge Theory Simulations in the Quantum Information Era*, Contemp. Phys. **57**, 388 (2016).

[176] E. Zohar, A. Farace, B. Reznik, and J. I. Cirac, *Digital Lattice Gauge Theories*, Phys. Rev. A **95**, 023604 (2017).

[177] E. Zohar, A. Farace, B. Reznik, and J. I. Cirac, *Digital Quantum Simulation of $Z_2$ Lattice Gauge Theories with Dynamical Fermionic Matter*, Phys. Rev. Lett. **118**, 070501 (2017).

[178] C. Schweizer, F. Grusdt, M. Berngruber, L. Barbiero, E. Demler, N. Goldman, I. Bloch, and M. Aidelsburger, *Floquet Approach to 2 Lattice Gauge Theories with Ultracold Atoms in Optical Lattices*, Nat. Phys. **15**, 1168 (2019).

[179] L. Barbiero, C. Schweizer, M. Aidelsburger, E. Demler, N. Goldman, and F. Grusdt, *Coupling Ultracold Matter to Dynamical Gauge Fields in Optical Lattices: From Flux Attachment to 2 Lattice Gauge Theories*, Sci. Adv. **5**, eaav7444 (2019).

[180] M. C. Bañuls and K. Cichy, *Review on Novel Methods for Lattice Gauge Theories*, Rep. Prog. Phys. **83**, 024401 (2020).

[181] Z. Davoudi, M. Hafezi, C. Monroe, G. Pagano, A. Seif, and A. Shaw, *Towards Analog Quantum Simulations of Lattice Gauge Theories with Trapped Ions*, Phys. Rev. Research **2**, 023015 (2020).

[182] X. Cui, Y. Shi, and J. C. Yang, *Circuit-Based Digital Adiabatic Quantum Simulation and Pseudoquantum Simulation as New Approaches to Lattice Gauge Theory*, J. High Energy Phys. 08 (2020) 160.

[183] D. Paulson, L. Dellantonio, J. F. Haase, A. Celi, A. Kan, A. Jena, C. Kokail, R. van Bijnen, K. Jansen, P. Zoller, and C. A. Muschik, *Towards Simulating 2D Effects in Lattice Gauge Theories on a Quantum Computer*, PRX Quantum **2**, 030334 (2021).

[184] G. Magnifico, T. Felser, P. Silvi, and S. Montangero, *Lattice Quantum Electrodynamics in (3 + 1)-Dimensions at Finite Density with Tensor Networks*, Nat. Commun. **12**, 3600 (2021).

[185] D. Poland, S. Rychkov, and A. Vichi, *The Conformal Bootstrap: Theory, Numerical Techniques, and Applications*, Rev. Mod. Phys. **91**, 015002 (2019).

[186] S. El-Showk, M. F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin, and A. Vichi, *Solving the 3D Ising Model with the Conformal Bootstrap*, Phys. Rev. D **86**, 025022 (2012).

[187] F. Kos, D. Poland, and D. Simmons-Duffin, *Bootstrapping Mixed Correlators in the 3D Ising Model*, J. High Energy Phys. 11 (2014) 109.

[188] F. Kos, D. Poland, D. Simmons-Duffin, and A. Vichi, *Precision Islands in the Ising and $O(N)$ models*, J. High Energy Phys. 08 (2016) 036.

[189] A. M. Ferrenberg, J. Xu, and D. P. Landau, *Pushing the Limits of Monte Carlo Simulations for the Three-Dimensional Ising Model*, Phys. Rev. E **97**, 043301 (2018).

[190] To see conventional percolation exponents here, it is important that the geometrical objects we are considering are really clusters rather than strict loops: We have nodes where more than 2 occupied links connect at a site. If we adopted a definition where the geometrical objects were strictly looplike, we would obtain a different universality class for unoriented loops [128,191].

[191] A. Nahum, J. T. Chalker, P. Serna, M. Ortuno, and A. M. Somoza, *Phase Transitions in Three-Dimensional Loop Models and the $CP^{n-1} \sigma$ Model*, Phys. Rev. B **88**, 134411 (2013).

[192] X. Xu, J. Wang, J. P. Lv, and Y. Deng, *Simultaneous Analysis of Three-Dimensional Percolation Models*, Front. Phys. **9**, 113 (2014).