# Relations between $e$ and $\pi$: Nilakantha’s series and Stirling’s formula

V. Yu. Irkhin

###### Abstract

Approximate relations between $e$ and $\pi$ are reviewed, some new connections
being established. Nilakantha’s series expansion for $\pi$ is transformed to
accelerate its convergence. Its comparison with the standard inverse-factorial
expansion for $e$ is performed to demonstrate similarity in several first
terms. This comparison clarifies the origin of the approximate coincidence
$e+2\pi\approx 9$. Using Stirling’s series enables us to illustrate the
relations $\pi^{4}+\pi^{5}\approx e^{6}$ and $\pi^{9}/e^{8}\approx 10$. The
role of Archimede’s approximation $\pi=22/7$ is discussed.

The fundamental mathematical constants $e$ and $\pi$ are both transcendental,
and their simple combinations are probably transcendental too, except for
famous Euler’s identity $e^{i\pi}=-1$.

However, there exist very close numerical coincidences in some combinations of
$e$ and $\pi$, which seem to be accidental. Several examples can be found in
[1, 2]:

$\pi^{2}=0.9996(4e-1)\approx 4e-1$ (1) ${163}(\pi-e)=68.99966...$ (2)
$(\pi^{4}+\pi^{5})/e^{6}=0.999999956...\,,$ (3)
${\pi^{9}/e^{8}=9.9998...\approx 10},$ (4) $e^{\pi}-\pi=19.999...,$ (5)
${\pi^{2}{(\pi-e)^{3/2}}/e=0.9999869...,}$ (6)

For some of such approximate relations, complicated mathematical
justifications were found. For example, “Ramanujan’s constant”

${e^{\pi{\sqrt{163}}}=262\,537\,412\,640\,768\,743.999\,999\,999\,999\,25\ldots\approx
640\,320^{3}+744}$ (7)

is very close to integer value owing to that 163 is a Heegner number [3].

##  1 Nilakantha’s series

First we focus on the simplest relation

$e+2\pi=9.001...,$ (8)

or $2(\pi-3)\approx 3-e$. We demonstrate a reason of this coincidence by using
the series expansions. For $e$ we use the standard inverse-factorial series
which may be represented as

$e=\sum_{n=0}^{\infty}\frac{1}{n!}=3-\frac{1}{3}+\frac{1}{24}+\frac{1}{120}+\frac{1}{720}+...$
(9)

and is rapidly convergent.

As for $\pi$, the situation is more complicated. The simple Gregory-Leibniz
expansion

$\frac{\pi}{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+...,$ (10)

converges very slowly. We apply Nilakantha’s series which can be obtained from
(3) after acceleration of the convergence:

$\displaystyle\pi$ $\displaystyle=$ $\displaystyle
3+\sum_{n=1}^{\infty}(-1)^{n+1}\left(\frac{1}{n+1}+\frac{1}{n}-\frac{4}{2n+1}\right)$
(11) $\displaystyle=$ $\displaystyle
3+\sum_{n=1}^{\infty}(-1)^{n+1}\frac{1}{n(2n+1)(n+1)}$ $\displaystyle=$
$\displaystyle 3+\sum_{n=1}^{\infty}(-1)^{n+1}\frac{4}{(2n+1)^{3}-(2n+1)}.$

The series (11) was discovered in India already in the 15th century and
converges rather rapidly [1].

After grouping of pairs of opposite-sign terms, (11) can be transformed as

$\displaystyle 2\pi$ $\displaystyle=$ $\displaystyle
6+\frac{1}{3}-\sum_{n=1}^{\infty}\frac{3}{n(n+1)(4n+1)(4n+3)}$ (12)
$\displaystyle=$ $\displaystyle
6+\frac{1}{3}-\frac{3}{70}-\frac{1}{198}-\frac{1}{780}-...$

so that the convergence becomes still faster.

One can see that the expansions (9) and (12) demonstrate a close similarity of
several first terms, which justifies and clarifies the coincidence (8).

The structure of Nilakantha’s series also reproduces roughly the relation (1).
We have from (12)

$\displaystyle\frac{1}{4}(\pi^{2}+1)\approx\frac{1}{4}\left(3+\frac{1}{6}-\frac{1}{48}-\frac{1}{200}\right)^{2}+\frac{1}{4}$
$\displaystyle\approx 3-\frac{1}{3}+\frac{1}{12}-\frac{1}{32}-0.01\approx
3-\frac{1}{3}+\frac{1}{24}$ (13)

which corresponds to the expansion of $e$ (9).

The relation (8) enables us to present combinations of $e$ and $\pi$ in terms
of one constant and therefore can be useful to clarify other above relations
between $e$ and $\pi$. So, the relation (1) is transformed as

$\pi^{2}+8\pi=35$ (14)

which gives

$\pi=\sqrt{51}-4\approx 3.1414...$ (15)

Expanding, we obtain

$\pi=7\sqrt{1+2/49}-4\approx 22/7\approx 3.143,$ (16)

i.e., Archimede’s number.

On substituting (8) into (2) we have

$\pi=\frac{512}{163}=\frac{2^{9}}{163}=3.1411...,$ (17)

which corresponds to Stoschek’s approximation using powers of two and the
number 163 that is the largest Heegner number [2].

##  2 Stirling’s formula

Another way to relate $e$ and $\pi$ is using Stirling’s series [5] which gives
a good accuracy even for $n=1$:

$\displaystyle n!$ $\displaystyle\simeq{\sqrt{2\pi
n}}\left({\frac{n}{e}}\right)^{n}\left(1+{\frac{1}{12n}}+{\frac{1}{288n^{2}}}+\cdots\right)$
(18)

so that for any $n$

$e^{n}\simeq{\sqrt{2\pi
n}}\frac{n^{n}}{n!}\left(1+{\frac{1}{12n}}+{\frac{1}{288n^{2}}}+\cdots\right)$
(19)

Eq.(19) enables one to generate various approximations for $e$ and $\pi$. For
example, we obtain

$e=\frac{e^{n+1}}{e^{n}}=\left(1+\frac{1}{n}\right)^{n+1/2}\frac{1+\frac{1}{12(n+1)}+\cdots}{1+\frac{1}{12n}+\cdots}.$
(20)

For half-integer $n$ in (19), using the value of Euler’s Gamma function
$\Gamma(1/2)=\sqrt{\pi}$ we derive

$e^{n+1/2}=\sqrt{2}\frac{(2n+1)^{n+1}}{(2n+1)!!}\left(1+{\frac{1}{6(2n+1)}}+\cdots\right)$
(21)

which yields a good approximation even for $n=0$, $e=49/18\approx 2.72$.

To demonstrate the usefulness of the relation (19) we substitute $e$ from (19)
with $n=1$ into (3). Expanding up to $1/n^{2}$ we obtain

$\pi^{2}+\pi=8(1+1/2+1/8+\cdots)\approx 13.$ (22)

Then we derive again Archimede’s number:

$\pi=-1/2+(7/2)\sqrt{1+4/49}\approx 3+1/7$ (23)

Otherwise, we can substitute (1) into (22) to obtain

$4e+\pi=14.$ (24)

Solving the system of linear equations (8), (24) we derive

$\pi=3+1/7,\,e=3-2/7.$ (25)

Combining (8) and (24), we can generate approximate equations with decreasing
accuracy:

$3\pi-2e=4,\,3e-\pi=5,\,5e-4\pi=1...$

and also for any combinations $n\pi+me$ with integer $n$ and $m$ satisfying
$n-2m=7k$.

Another way to derive (22) from (3) is to take into account
$e^{3}=20.08\approx 20$ and the well-known relation $\pi^{3}\approx 31$ which
follows from the series

$\pi^{6}=960\sum_{n=0}^{\infty}\frac{1}{(2n+1)^{6}}\approx 960.$ (26)

Then we have

$\pi^{2}+\pi\approx 400/31\approx 13.$

Although the accuracy of (15), (23) is somewhat worse than that of original
relations (1), (3), they demonstrate close connection with ancient
investigations of the problem of calculating $\pi$.

Now we are also able to clarify the relation (4). Using (19) for both $n=1$
and $n=2$ we derive

$\displaystyle e^{8}$ $\displaystyle\simeq$ $\displaystyle
64\pi^{3}\left(1+{\frac{1}{12}}+\cdots\right)^{4}\left(1+{\frac{1}{24}}+\cdots\right)^{2}$
(27) $\displaystyle\approx$ $\displaystyle\frac{3}{2}64\pi^{3}=96\pi^{3}.$

Using (26) wee see the coincidence of (27) and (4).

The coincidence (5) can be related to the continued-fraction representation of
Gelfond’s constant, which also starts from Archimede’s approximation [6]:

$e^{\pi}=[23;7,9,1,1,591...]\approx 20+3+1/7+...$ (28)

In more detail, we can use $e^{3}\approx 20$ to write down

$e^{\pi}=e^{3}e^{\pi-3}\approx 20(1+1/7+0.01)\approx 20(1+\pi/20)=20+\pi.$
(29)

The relation (6) is most difficult. To obtain an estimation, we take into
account (8) and rewrite (6) as

$27\pi^{8}(\pi-3)^{3}/(\pi^{2}e)^{2}\approx 1.$ (30)

Further on, we use Archimede’s $\pi-3\approx 1/7$, $7^{3}\approx 350$,
$\pi^{2}e\approx 27$ (which follows from the geometric mean corresponding to
(8)), and

$\pi^{8}=9450\,\zeta(8)\approx 9450=350\cdot 27$

with $\zeta(x)$ the Riemann zeta function. Then Eq.(30) is satisfied.

Although the subject of the present work can seem somewhat artificial, we note
that combinations of $e$ and $\pi$ occur in physical problems, e.g., in the
exact solution of the Kondo model [7].

## References

  * [1] E. W. Weisstein, Almost Integer, MathWorld–A Wolfram Web Resource, https://mathworld.wolfram.com/AlmostInteger.html
  * [2] E. W. Weisstein, Pi Approximations, MathWorld–A Wolfram Web Resource, https://mathworld.wolfram.com/PiApproximations.html
  * [3] J. D. Barrow, The Constants of Nature. London: Jonathan Cape, 2002.
  * [4] R. Roy, The Discovery of the Series Formula for $\pi$ by Leibniz, Gregory and Nilakantha. Mathematics Magazine, Vol. 63, No. 5, pp. 291-306 (1990).
  * [5] Olver, F. W. J.; Olde Daalhuis, A. B.; Lozier, D. W.; Schneider, B. I.; Boisvert, R. F.; Clark, C. W.; Miller, B. R.; Saunders, B. V., 5.11 Gamma function properties: Asymptotic Expansions, NIST Digital Library of Mathematical Functions, Release 1.0.13 of 2016-09-16.
  * [6] The On-Line Encyclopedia of Integer Sequences, A058287.
  * [7] A.M. Tsvelick, P.B. Wiegmann, Exact results in the theory of magnetic alloys, Adv. Phys. 32, 453 (1983).


Pushing Blocks by Sweeping Lines

Hugo A. AkitayaUniversity of Massachusetts Lowell<EMAIL_ADDRESS>LöfflerDepartment of Information and Computing Sciences, Utrecht University<EMAIL_ADDRESS>VigliettaSchool of Information Science, Japan Advanced Institute of Science and Technology (JAIST), Japan<EMAIL_ADDRESS>
We investigate the reconfiguration of $n$ blocks, or “tokens”, in the square grid using line pushes. A line push is performed from one of the four cardinal directions and pushes all tokens that are maximum in that direction to the opposite direction by one unit. Tokens that are in the way of other tokens are displaced in the same direction, as well.

Similar models of manipulating objects using uniform external forces match the mechanics of existing games and puzzles, such as Mega Maze, 2048 and Labyrinth, and have also been investigated in the context of self-assembly, programmable matter and robotic motion planning.
The problem of obtaining a given shape from a starting configuration is know to be NP-complete.

We show that, for every $n$, there are sparse initial configurations of $n$ tokens (i.e., where no two tokens are in the same row or column) that can be compacted into any $a\times b$ box such that $ab=n$. However, only $1\times k$, $2\times k$ and $3\times 3$ boxes are obtainable from any arbitrary sparse configuration with a matching number of tokens.
We also study the problem of rearranging labeled tokens into a configuration of the same shape, but with permuted tokens.
For every initial configuration of the tokens, we provide a complete characterization of what other configurations can be obtained by means of line pushes.

§ INTRODUCTION

§ DEFINITIONS AND PRELIMINARIES

§ COMPACTION PUZZLES

§ PERMUTATION PUZZLES

§ CONCLUSIONS AND OPEN PROBLEMS


# Anderson localization of emergent quasiparticles: Spinon and vison interplay
at finite temperature in a $\mathbb{Z}_{2}$ gauge theory in three dimensions

Minho Kim  TCM group, Cavendish Laboratory, University of Cambridge. Cambridge
CB3 0HE, United Kingdom  Giuseppe De Tomasi  Department of Physics, University
of Illinois at Urbana-Champaign, Urbana, Illinois 61801-3080, USA  Claudio
Castelnovo  TCM group, Cavendish Laboratory, University of Cambridge.
Cambridge CB3 0HE, United Kingdom

###### Abstract

Fractional statistics of quasiparticle excitations often plays an important
role in the detection and characterization of topological systems. In this
paper, we investigate the case of a three-dimensional (3D) $\mathbb{Z}_{2}$
gauge theory, where the excitations take the form of bosonic spinon
quasiparticle and vison flux tubes, with mutual semionic statistics. We focus
on an experimentally relevant intermediate temperature regime, where sparse
spinons hop coherently on a dense quasistatic and stochastic vison background.
The effective Hamiltonian reduces to a random-sign bimodal tight-binding
model, where both the particles and the disorder are borne out of the same
underlying quantum spin liquid (QSL) degrees of freedom, and the coupling
between the two is purely driven by the mutual fractional statistics. We study
the localization properties and observe a mobility edge located close to the
band edge, whose transition belongs to the 3D Anderson model universality
class. Spinons allowed to propagate through the quasistatic vison background
appear to display quantum diffusive behavior. When the visons are allowed to
relax, in response to the presence of spinons in equilibrium, we observe the
formation of vison depletion regions slave to the support of the spinon
wavefunction. We discuss how this behavior can give rise to measurable effects
in the relaxation, response and transport properties of the system and how
these may be used as signatures of the mutual semionic statistics and as
precursors of the QSL phase arising in the system at lower temperatures.

##  I  Introduction

Quantum spin liquids [1, 2, 3, 4] (QSLs) are intriguing phases of matter
characterized by highly correlated spins that fail to order in a conventional
way down to temperatures much lower than the characteristic interaction energy
scale in the system. This leads to the appearance of exotic properties,
including emergent gauge symmetries, fractionalized excitations with anyonic
statistics, and topological order [5, 6, 7].

Such phenomenology is both of fundamental interest and has potential
applications in quantum information storage and processing [8, 9, 10]. Despite
the growing number of candidate spin liquid materials discovered in
laboratories in recent years, their experimental detection and
characterization remains a paramount challenge to date [4].

The fractional behavior of the quasiparticle excitations in QSLs is often
considered a promising handle to make them apparent in experiments [11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Authors of recent work [22, 24,
23] focused on a class of systems — encompassing spin ice [25, 26], Kitaev
materials [27], and valence bond systems [28] — where the QSL behavior is
borne out of the interplay of some large energy scale that acts projectively
on the Hilbert space — inducing on its own a classical spin liquid state [29]
— and smaller transverse (i.e., kinetic) terms [30, 31]. Considering two-
dimensional (2D) models, the presence of a finite temperature regime was
brought to light, where the anyonic statistics between the quasiparticles
leads to a rich phenomenology driven by Anderson localization [32] of emergent
tight binding particles, where the disorder is itself emergent from the same
spin degrees of freedom in an otherwise disorder-free system [33, 34, 35, 36,
37]. The resulting effects on relaxation and transport properties [22, 23] are
both of fundamental interest in their own right as well as of potential
experimental interest as precursor diagnostics tools and bellwethers for QSL
phases at lower temperatures.

In this paper, we revisit the intriguing finite-temperature juncture of
topological behavior and Anderson localization of emergent quasiparticles in
disorder-free systems by extending it to three dimensions (3D). This brings
about a few important changes, whose study is the key interest of this paper.
On the one hand, nontrivial mutual statistics in dimensions higher than 2D
requires at least one of the quasiparticle species to take the form of an
extended rather than pointlike object (e.g., a closed loop, in our system of
choice). On the other hand, 3D is the smallest dimension in which the system
has a metal-insulator Anderson transition as a function of energy, producing a
so-called mobility edge [38] that is absent in the 2D model considered earlier
[22, 23].

Specifically, we investigate a $\mathbb{Z}_{2}$ topological spin liquid system
on the cubic lattice [39, 40], which exhibits two species of fractionalized
excitations (spinons and visons) with mutual semionic statistics. The spinons
are pointlike (i.e., end points of open strings), and we assume their energy
cost to be the largest term in the Hamiltonian, with their transverse or
hopping term as the first subleading energy scale. As a first approximation,
the spinons can effectively be modeled as sparse tight-binding particles on
the lattice. The visons take instead the form of closed loops (i.e., edges of
open membranes embedded in 3D) and we assume their energy cost as well as
transverse or hopping terms to be negligibly small (effectively vanishing)
compared with the temperature. As such, the visons can effectively be modeled
as classical stochastic objects that are quasistatic on the characteristic
transport and equilibration timescales of the spinons. The visons therefore
act as a stochastic background that affects the behavior of the spinons
through mutual semionic statistics — a behavior that is equivalent to charges
moving on a lattice where plaquettes are randomly threaded by half quantum
($\pi$) fluxes [22]. Vice versa, the spinons endow the vison configuration
with an energy that can affect their correlations in an important way [23].

In the first instance, we study the localization properties of a spinon moving
across a stochastic and static vison background. This is equivalent to an
electron hopping on a cubic lattice with randomly distributed $\pi$ fluxes
threading its plaquettes — a model that has hitherto not been considered in
the literature. Whereas in 2D all states are localized away from the middle of
the band, we find clear evidence of a mobility edge in 3D, which occurs very
close to the band edge. We find that the Anderson transition of our model
belongs to the same universality class as the orthogonal 3D Anderson model
[41, 42], and the states inside the mobility edges are indeed extended.

We then proceed to study the behavior of spinons and visons, following the
path set by Refs. 22, 23. Despite the differences in the properties of the
eigenstates in 3D, the behavior we observe is very similar to the 2D case. We
find slightly weaker localization effects, as reflected by the fact that the
spinon propagation behavior in the $\pi$ flux background in 3D agrees well
with the diffusive like continuous random flux case — unlike in 2D where the
propagation showed a remarkable departure toward anomalous diffusion [22].
Notwithstanding, once the visons are allowed to move and relax stochastically
in the presence of the spinons, we do observe once again the formation of the
characteristic depletion regions [23] that are responsible for the predicted
anomalous relaxation and transport effects. This is once again in stark
contrast with the case where mutual statistics is replaced by interactions
between spinons and visons, where the formation of such depletion regions is
not observed.

Our results show that the rich phenomenology arising at the interface between
topological order and Anderson localization in disorder-free systems is not
limited to low-dimensional systems. Moreover, while the localization
properties are different, the notable effects on the relaxation and transport
properties underpinned by the formation of vison depletion regions around the
spinons survive. This is once again a behavior due to the mutual semionic
statistics and, if observed, can thus be taken as a signature of 3D QSL
behavior at finite temperature.

While we must wait for the discovery of candidate materials that realize
$\mathbb{Z}_{2}$ QSL phases in 3D, our results may be relevant to other
contexts, including frustrated magnetic pyrochlore oxides and resonant valence
bond systems. Moreover, the possibility of realizing $\mathbb{Z}_{2}$ spin
liquid Hamiltonians in our temperature regime with quantum annealers [43, 24,
44] and quantum simulators [45, 46] could provide a suitable arena where the
physics discussed here could be tested and explored further.

The outline of the paper is as follows. In Sec. II, we introduce the 3D model
and its Hamiltonian as well as the assumptions made for the temperature regime
of interest in this paper. In Sec. III, we discuss the localization properties
and verify the existence of a mobility edge in the system. We then proceed to
study the diffusion of spinons in a static random vison background in Sec. IV,
and we compare it with analytical results based on a Bethe lattice
approximation [23]. In Sec. V, we allow the visons to relax stochastically
while the spinons remain instantaneously in thermodynamic equilibrium. As in
2D, we observe the formation of mutual-semionic-statistics-induced vison
depletion regions around each spinon. We describe the effects that this
behavior has on the response and transport properties of the system,
particularly on the out-of-equilibrium dynamics when the temperature is varied
quickly enough, in Sec. VI. Finally, we present our conclusions and outlook in
Sec. VII.

##  II  Model

In this section we introduce the 3D toric code model [39, 40] (see also
Appendix D in Ref. 47) and the corresponding effective quasiparticle
Hamiltonian. The model is defined on an $L^{3}$ cubic lattice with periodic
boundary conditions. Spin-$\frac{1}{2}$ degrees of freedom $\sigma_{i}$ live
on the bonds of the lattice, labeled by the index $i=1,2,\cdots,3L^{3}$ (see
Fig. 1).

(a)  (b)

Figure 1:  Visualization of the lattice system considered in this paper, with
spin-$\frac{1}{2}$ degrees of freedom living on the bonds of a cubic lattice
(black solid dots). Star and plaquette operators at site $s$ and plaquettes
$p_{1},p_{2},p_{3}$, respectively, are shown for illustration purposes (left
panel). Acting on any given spin (e.g., the red dot labeled $\sigma$ in the
right panel) can affect the state of the two adjacent star operators (not
shown), and/or the state of the four adjacent plaquette operators
($B_{p_{1}}$, $B_{p_{2}}$, $B_{p_{3}}$, and $B_{p_{4}}$ in the figure).

We further label the elementary square plaquettes on the lattice with the
index $p=1,2,\cdots,3L^{3}$, and the lattice sites with $s=1,2,\cdots,L^{3}$.

We introduce the plaquette and star operators:

$B_{p}=\prod_{i\in p}\sigma_{i}^{z}\qquad{\rm and}\qquad A_{s}=\prod_{i\in
s}\sigma_{i}^{x}\,,$ (1)

where for convenience of notation $i\in p$ and $i\in s$ indicate the four
spins around a plaquette $p$ and the six spins adjacent to a site $s$,
respectively. Examples of spin and plaquette operators are shown in Fig. 1.

The Hamiltonian reads

$H=-\frac{\Delta_{s}}{2}\sum_{s}A_{s}-\frac{\Delta_{v}}{2}\sum_{p}B_{p}\,,$
(2)

where $\Delta_{s},\Delta_{v}>0$. All plaquette and star operators commute with
one another, and the ground state of the toric code Hamiltonian is obtained
when they all have eigenvalue $+1$. Negative eigenvalues of $A_{s}$ and
$B_{p}$ label the excited states. Both types of excitations are gapped; the
former are dubbed star defects or spinons (energy cost $\Delta_{s}$), whereas
the latter are dubbed plaquette defects or visons (energy cost $\Delta_{v}$).
Notice that both types of excitations are static according to Eq. (2). In the
following, it will be useful to introduce a vison occupation variable
$n_{p}=0,1$ [namely, $n_{p}=(1-B_{p})/2$].

Changing the sign of the $\sigma^{x}$ component of a spin on a given bond $i$
results in flipping the sign of the expectation values of the pair of star
operators on the $s,s^{\prime}$ sites adjacent to $i$. This amounts to the
creation/annihilation of a pair of spinons or to the hopping of a spinon
between the two sites $s$ and $s^{\prime}$, depending on the local initial
state of the spins. This directly relates to the well-known gauge symmetry in
the $\sigma^{x}$ degrees of freedom, whereby spinons act as sinks/sources of
the gauge field.

Similarly, changing the sign of the $\sigma^{z}$ component of a spin on a
given bond $i$ results in flipping the sign of the expectation values of the
four plaquette operators on the plaquettes adjacent to $i$ (see Fig. 1, right
panel). This leads to an equivalent gauge structure that is, however,
formulated in a redundant way, with an extensive number of conserved
quantities [40]. As a consequence of the cooperative flipping of the
eigenvalues of plaquette operators, visons must necessarily form edges of open
2D membranes embedded in the dual 3D cubic lattice (see Refs. 39, 40 for
details). This is indeed a key feature that allows spinons and visons to
exhibit nontrivial mutual statistics in 3D, as the operation of winding a
spinon around a vison remains well defined: as a spinon moves along a closed
loop on the cubic lattice, it either does or does not cross a given vison
membrane (an odd number of times); if it does, the final state acquires an
overall minus sign with respect to the starting one — spinons and visons obey
mutual semionic statistics. On their own, instead, both spinons and visons are
bosonic excitations.

In this paper, we are interested in a modification of the toric code
Hamiltonian in Eq. (2), where a longitudinal field $h\sum_{i}\sigma^{z}_{i}$,
$h<\Delta_{s}$, is added and the vison gap $\Delta_{v}$ is taken to be small
enough to be negligible (e.g., much smaller than some thermal energy scale
$0<T\ll\Delta_{s}$). The effects of the field are two fold: (i) spinon
excitations become itinerant and form a dispersive band, and (ii) the vison
gap gets perturbatively renormalized (which is a small effect and does not
conflict with the assumption above that it is negligible, at least at small
but finite temperatures). For a justification of these assumptions in relation
to experimentally relevant models of spin liquids, we refer the reader to the
introductory sections in Refs. 22, 23.

In this regime, the spinons behave as quantum itinerant quasiparticles that
can be effectively described in the sparse limit by a tight-binding
Hamiltonian [22]. The visons, on the other hand, remain static. Given the
mutual semionic statistics between spinons and visons [8], the presence of the
latter means that the spinon wavefunction in the effective tight-binding model
picks up a minus sign as the spinon hops around a closed loop on the lattice
traversing a vison membrane an odd number of times. This is equivalent to the
behavior of a unit charge moving around a loop threaded by half an elementary
magnetic flux unit ($\pi$ flux). The vison occupancy $n_{p}=0,1$ defines the
$\pi$ flux pattern through which spinons move. It is therefore convenient to
rewrite the effective Hamiltonian as

$H_{\textrm{eff}}=\Delta_{s}\sum_{s}b_{s}^{\dagger}b_{s}-\eta_{s}\sum_{\langle
s,s^{\prime}\rangle}\left[\exp(iA_{ss^{\prime}})b_{s}^{\dagger}b_{s^{\prime}}+{\rm
H.c.}\right]\,,$ (3)

where $(\nabla\times A)_{p}=\Phi_{p}=\pi\,n_{p}$ and
$A_{ss^{\prime}}=-A_{s^{\prime}s}$. Here, $b_{s},\;b_{s}^{\dagger}$ are the
spinon annihilation and creation operators at lattice site $s$, which obey
bosonic statistics. The second sum runs over pairs of nearest-neighbor lattice
sites, and we relabeled the field $h\to\eta_{s}$ for notational convenience in
the new language. One should further constrain the spinon density operator
$b_{s}^{\dagger}b_{s}$ to take on values $0,1$ only and include spinon
annihilation and creation terms. However, in this paper we will primarily
concern ourselves with the behavior of isolated spinons in the low spinon
density limit. Therefore, we will safely ignore the constraint as well as the
additional terms.

Generally speaking, one may expect small perturbations in the system to give
rise to hopping amplitudes for visons as well (the visons ought to move
subject to the gauge constraint discussed above, which is not equally amenable
to a pointlike particle hopping description, as is the case for the spinons).
In this paper we assume that such terms are small compared with the already
small vison gap $\Delta_{v}$ (otherwise, they would drive the system out of
the topological phase that underpins the physics of interest in this paper)
and hence also small compared with temperature. In this regime, coherent vison
dynamics is unlikely to survive on timescales of significance to the results
discussed in this paper; therefore, where appropriate, we shall model the
vison dynamics stochastically in a Markovian approximation (namely, using
Monte Carlo (MC) simulations that respect the gauge structure of the system).
Moreover, we are justified to take a Born-Oppenheimer perspective and
approximate the visons to be quasistatic on the timescale of
motion/equilibration of the spinons; namely, at short times, we can assume the
visons to be static (Sec. IV) with respect to spinon motion, whereas at long
times, we can assume the spinons to be instantaneously in equilibrium on the
timescales over which the vison stochastic dynamics occurs (Sec. V).

##  III  Energy eigenstate localization properties

As discussed in Sec. II, see, e.g., Eq. (3), visons introduce random $\pm$
phases in the spinon hopping amplitudes $\eta_{s}$. One can therefore draw a
parallel with random-hopping tight-binding particles in 3D lattices. However,
exploration of the case of bimodal hopping ($\pm\eta_{s}$) where the
distribution of signs is induced by an emergent $\mathbb{Z}_{2}$ gauge
structure is lacking. For this reason, we include here a brief study of the
localization properties of the energy eigenstates of our system while leaving
a more thorough characterization thereof to future work.

Anderson localization with off-diagonal disorder is somehow unique [48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60]. The presence of chiral
(sublattice) symmetry leads to a divergence of the density of states (DOS) at
the band center ($E\approx 0$). Consequently, the eigenstates around the
divergence might behave differently from all others. For instance, in one
dimension, all eigenstates away from $E\approx 0$ are exponentially localized
111The case of random phases is special in one dimension. Indeed, using gauge
transformation $b_{s}->\exp(i\phi_{s})b_{s}$ one can absorb the phase..
However, around $E\approx 0$, there exists a divergence in the DOS which
scales as $\rho(E)\sim\frac{1}{E\ln^{3}{E}}$ [48, 54], and the eigenstates
with $E\approx 0$ are only weakly localized [50, 62, 51] with $\psi(E\approx
0)\sim\exp(-\gamma\sqrt{r})$ 222Formally, the localization length is
infinite.. In higher dimensions, the question of localization at the band
center is still unresolved. At large enough disorder, a general argument
predicts $\psi(E\approx 0)\sim\exp({-\gamma\sqrt{\ln{r}}})$ [51]. However,
finite-sized numerical simulations found the existence of power law localized
states [56].

In addition to the localization behavior for states around $E\approx 0$, the
general belief is that off-diagonal systems are like the ones with onsite
disorder. Thus, the critical dimension is still $D=2$, and for larger
dimensions, one expects an Anderson transition as a function of disorder
strength and the existence of mobility edges within the extended phase. In our
case, the fluctuations of the magnitude of the disorder is fixed, i.e.,
$|\exp{iA_{s,s^{\prime}}}|=1$. This has two immediate consequences: first, we
do not have a tuning parameter for the magnitude of the disorder; and second,
the disorder is weak [53]. Therefore, we expect our system to be in the regime
where a mobility edge exists, and eigenstates close to the center of the
spectrum are delocalized, whereas states at the edge of the spectrum are
localized [53, 55, 59, 58, 57, 64, 65].

To probe the existence of mobility edges, we compute the generalized inverse
partition ratio IPRq [38], for several system sizes, focusing on $q=2$ and
$\frac{1}{2}$ and the associated multifractal exponents:

${\rm IPR}_{q}=\sum_{i}\lvert\psi_{i}\rvert^{2q}\qquad
D_{q}=\frac{1}{1-q}\frac{\ln({\rm IPR}_{q})}{\ln L^{3}}\,,$ (4)

where $L^{3}$ is the total number of sites on the cubic lattice, and
$\psi_{i}$ is the spinon eigenstate wave function at site $i$. In the
localized phase, both IPR${}_{2,\frac{1}{2}}$ saturate to a constant value as
a function of $L$, i.e., IPR${}_{2,\frac{1}{2}}\sim\mathcal{O}(L^{0})$, for
$L\gg 1$. Instead, in the extended phase, we expect IPR${}_{q}\sim
L^{3(1-q)}$, for $L\gg 1$. In the following, we shall denote disorder-averaged
quantities with an overline, e.g., $\overline{D_{q}}$. From Eq. (4), one
clearly sees that the multifractal exponent tends to $0$ when the sampled
states are localized and tends to $1$ when the sampled states are delocalized,
for both $q=\frac{1}{2}$ and $2$ in the thermodynamic limit $L\to\infty$.

The overall behavior of $\overline{D_{q}}$ across the energy spectrum is shown
in Fig. 2, for several $L$ as a function of $E$. The data were obtained using
exact diagonalization, averaging over vison configurations 333Here, $1000$
vison configurations for $L=15$, $500$ for $L=23$, and $100$ vison
configurations for larger system sizes. and a small energy shell around the
target energy 444The width of the energy shell has been chosen to be small
enough for the change in the density of states to be negligible..

We observe a general tendency toward delocalization in the middle of the
spectrum and toward localization at the edge (akin to what was found for
random bond disorder [53, 55, 59, 57]), providing evidence of the existence of
mobility edges. Upon closer inspection (see Fig. 2), our results suggest the
presence of a mobility edge around $E_{c}\sim-4.35$ (in units of the spinon
hopping). Indeed, for $|E|>|E_{c}|$, $\overline{D_{q}}$ tends to zero with
increasing $L$, and for $|E|<|E_{c}|$, we have $\overline{D_{q}}\rightarrow
1$. To better pin down the value of $E_{c}$, we perform finite-scaling
analysis on $\overline{D_{q}}$ by collapsing the curves of $\overline{D_{q}}$
with different $L$. As shown in the insets in Fig. 2, a good collapse is
obtained using $E_{c}=-4.33$ for $q=2$ and $E_{c}=-4.37$ for $q=\frac{1}{2}$,
with the critical exponent being $\nu=1.6$ for both cases, in good agreement
with the critical exponent of the 3D Anderson model with onsite disorder
$\nu_{\mathrm{AT}}\approx 1.58$ [41, 42], belonging to the orthogonal
universality class. This is expected, given that our Hamiltonian is real and
symmetric and therefore has time-reversal symmetry. In contrast, when one
assumes random phases as off-diagonal disorder, time-reversal symmetry is
broken, and it is known that the localization transition belongs to the
unitary class [68]. We note that this region, being close to the band edge,
has a small DOS which is rapidly increasing as energy increases, as shown in
Appendix A.

Figure 2:  Disorder-averaged multifractal coefficients $\overline{D_{q}}$ for
$q=2$ (top panel) and $q=\frac{1}{2}$ (bottom panel) for various system sizes
from $L=10$ to $75$, near the edge of the spectrum (in units of the spinon
hopping). Our results suggest the existence of a mobility edge near
$E_{c}\sim-4.35$. The insets show the scaling with system size for
$\overline{D_{q}}$. The curves for various system sizes collapse onto the same
curve when scaled with the empirical critical exponent $\nu=1.6$.

To further support the existence of mobility edges, we analyze the energy
spectrum of the effective Hamiltonian $H_{\textrm{eff}}$ in Eq. (3). We focus
on the energy correlation between close-by energies, quantified by the so-
called $r$-gap ratio [69, 70]:

$r_{n}=\frac{\mathrm{min}(\delta_{n},\delta_{n-1})}{\mathrm{max}(\delta_{n},\delta_{n-1})},\qquad\delta_{n}=E_{n}-E_{n-1}\,.$
(5)

It is known that the mean value for $r$ should be $r_{\text{Poisson}}=2\log
2-1=0.386$ for localized states and $r_{\text{GOE}}=0.5295$ for Gaussian
orthogonal ensembles (GOEs, extended states) [69, 70]. As shown in Fig. 3, at
low energies, $r\approx r_{\text{Poisson}}$, and the system is localized,
while for energies closer to the band center, $r_{\text{GOE}}\approx 0.52$,
and the system is extended. Furthermore, in Fig. 3 (bottom panel), we perform
the same finite-sized scaling analysis for $r$, as in Fig. 2. The analysis of
the energy spectrum confirms $E_{c}\approx-4.37$ and the critical exponent
$\nu\approx\nu_{\mathrm{AT}}$.

Figure 3: $r$ values window-averaged with window width $\Delta E=0.1$ (top
panel), and the same data scaled using the empirical critical exponent
$\nu=1.6$. The $r$ values transition from the Poisson distribution value
(localized states, lower dotted line) to the Wigner value (extended states,
upper dotted line), with a critical point at $E/\eta_{s}\sim-4.37$ (vertical
dotted line). The results for various system sizes show a good collapse onto a
single curve when scaled by system size near the mobility edge, analogously to
Fig. 2, suggesting the existence of a localization-delocalization transition.

Finally, for the sake of completeness, we inspect the probability distribution
of $D_{2}$ and of the energy level spacing
$s=\frac{E_{n+1}-E_{n}}{\overline{E_{n+1}-E_{n}}}$, for two energy windows,
one in the localized phase and another in the extended phase (see Fig. 4). In
the localized phase (Fig. 4, left panels), the probability distribution of
$P(\log{D_{2}})$ shifts uniformly to lower values $\overline{D_{q}}\ll 1$ with
increasing $L$, ruling out the existence of delocalized states.
Correspondingly, the probability distribution of the level spacing is
Poissonian, $P(s)\sim\exp(-s)$. In the extended phase (Fig. 4, right panels),
both $P(\log{D_{2}})$ and $P(s)$ show full delocalization.

Figure 4:  The probability distribution of the multifractal coefficient
$D_{2}$ for localized ($E/\eta_{s}\in[-4.6,-4.5]$, left), and extended
($E/\eta_{s}\in[-4.3,-4.2]$, right) states, alongside the level statistics
$P(s),s=(E_{n+1}-E_{n})/\overline{(E_{n+1}-E_{n})}$ in the same energy
windows. In the localized regime, we see a Poissonian distribution of the
level statistics, while in the extended regime, we observe a Gaussian
orthogonal ensemble (GOE) distribution (Wigner surmise). The $r$ value for
each lattice size in the given energy window is also marked, showing once
again that the states are localized (left panel) and extended (right panel).
The dashed lines are the Poissonian distribution
$P_{\text{Pois.}}(s)=\exp(-s)$ (left panel) and Wigner surmise
$P_{\mathrm{GOE}}(s)=\frac{\pi s}{2}\exp(-\pi s^{2}/4)$ (right panel).

One should note that, at finite temperature, thermal occupation of excited
states and coupling to the bath can cause dephasing, which destroys
localization. This is driven by variable-range hopping [71], and the resulting
conduction current (in 3D) can be expressed as
$I=I_{s}\exp\left[-\left(T_{\mathrm{M}}/T\right)^{1/4}\right]$, where
$T_{\mathrm{M}}$ is the so-called Mott temperature. Then $I_{s}/I$ gives us an
estimate of the timescale up to which the results discussed above may be
expected to hold: $t_{\rm
max}\sim\exp\left[\left(T_{\mathrm{M}}/T\right)^{1/4}\right]$. In a localized
state, $T_{\mathrm{M}}\propto[\xi_{\mathrm{loc}}\rho(E)]^{-1}$ [71], where
$\rho(E)\sim 10^{-2}\eta_{s}^{-1}$ is the DOS in the localized phase and
$\xi_{\mathrm{loc}}\sim 1$ is the localization length. Therefore, in
comparison with the temperatures of interest $T/\eta_{s}<10^{-2}$, the Mott
temperature is relatively large $T_{M}/\eta_{s}\gtrsim 10^{2}$, and
localization survives up to long times $t_{\rm max}\gtrsim 10^{4}$.

##  IV  Spinon propagation

In this section, we look at the case of a single spinon moving in a static
random vison background (with probability $\frac{1}{2}$ of vison occupancy
$n_{p}=0,1$, or equivalently, density $\langle n_{p}\rangle=\frac{1}{2}$),
which corresponds to the model introduced in Sec. II where the energy cost of
visons is set to zero (i.e., an infinite temperature regime for the visons)
555We remark that our choice of random vison background respects the emergent
gauge structure of the system, discussed in Sec. II.. This regime is relevant
to the short-time, short-distance regime where the visons are treated
quasistatically with respect to the motion of the spinons.

From the time dependence of the spinon wave function, we compute the root mean
square displacement:

$\langle\mathbf{r}^{2}(t)\rangle=\sum_{i}d^{2}(i,i_{0})|\langle
i|\exp(-iH_{\textrm{eff}}t)|i_{0}\rangle|^{2},$ (6)

where $d(i,i_{0})$ is the Euclidean distance between the initial site $i_{0}$
and site $i$. To compute $\langle\mathbf{r}^{2}(t)\rangle$, we employ fourth-
order Trotterized [73] time evolution on a $L^{3}=149^{3}$ cubic lattice, with
open boundary conditions. The initial state is chosen to be a spinon wave
function localized on a single lattice site at the center of the lattice.

The results are shown in Fig. 5, where we see a crossover from short-range,
short-time ballistic propagation $\langle\mathbf{r}^{2}(t)\rangle\sim t^{2}$
to diffusivelike behavior $\langle\mathbf{r}^{2}(t)\rangle\sim t$. This agrees
with the previous results that our system belongs to the 3D universality
class, which is characterized by a diffusive behavior within the extended
phase [74, 75, 76, 77, 78].

Figure 5:  Numerical simulations of the spinon mean square displacement
$\langle\mathbf{r}^{2}(t)\rangle$ (blue line) compared with the analytical
solution on the Bethe lattice (red line). We used fourth-order Suzuki-Trotter
evolution on a $149^{3}$ lattice, averaged over $128$ different vison
background realizations. The inset shows the same data in linear scale, to
show more clearly the departure of the numerics from quantum diffusive
behavior.

The numerical results are contrasted with a self-retracing approximation that
allows for an analytical solution on the Bethe lattice, as discussed in Ref.
22. Details of the Bethe lattice calculation are given in Appendix B. Here, we
report and discuss only the main results, namely, the equation for the $2k$th
moment $\langle\mathbf{r}^{2k}(t)\rangle$:

$\begin{split}\left\langle\mathbf{r}^{2k}(t)\right\rangle&=\oiint\frac{dw_{1}}{2\pi
i}\frac{dw_{2}}{2\pi
i}\frac{\exp[i\eta_{s}t\left(w_{1}-w_{2}\right)]}{w_{1}w_{2}}T\left(\frac{1}{w_{1}}\right)T\left(\frac{1}{w_{2}}\right)\\\
&\times\mathcal{R}_{2k}\left[S\left(\frac{1}{w_{1}}\right)S\left(\frac{1}{w_{2}}\right)\right]\,,\end{split}$
(7)

where $T(x)$ and $S(x)/x$ are generating functions for self-retracing paths
from the root node (depth $0$) and nonroot node (depth $\geq 1$) for a Bethe
lattice of coordination number $z$, and $\mathcal{R}_{2k}(x)$ is the generator
for the $2k$th moment. The structure of $\mathcal{R}_{2k}(x)$ alongside $T(x)$
and $S(x)$ (see Appendix B) implies that the pole signatures are the same as
in 2D, with a line of double poles at $w_{1}=u+i0^{+},w_{2}=u-i0^{-}$ for
$u\in\mathbb{R}$ and $|u|<2\sqrt{z-1}$. While there are subtle differences due
to the additional lattice dimension (2D $\to$ 3D), the resulting identical
pole signatures give the same power law dependence as in 2D; hence, we expect
$\langle\mathbf{r}^{2}\rangle\sim t^{2}$ at short times (ballistic
propagation) and $\langle\mathbf{r}^{2}\rangle\sim t$ at long times (quantum
diffusion), with the crossover happening at $\eta_{s}t\sim 1$. This is indeed
observed in Fig. 5. In the self-retracing approximation, the 3D cubic lattice
asymptotically exhibits the same diffusion coefficient as the 2D triangular
lattice:

$D_{6}^{\mathrm{diff}}=\frac{3}{\pi}\left[3\sqrt{5}-2\ln\left(\frac{3+\sqrt{5}}{3-\sqrt{5}}\right)\right]\simeq
2.72968\,.$ (8)

In this paper, we use the superscript diff to unequivocally distinguish the
diffusion constant $D_{k}^{\mathrm{diff}}$, with $k$ being the coordination
number of the lattice, from the multifractal exponent $D_{q}$. The full
solution of the self-retracing approximation (Fig. 5) was obtained
numerically, using the $\lvert w_{1}\rvert=\lvert w_{2}\rvert=6$ torus (we are
performing a double-loop integral; hence, the integration region is a torus,
not a loop). This integration region encircles all the double-pole
singularities. A numerically exact integral is achievable only up to
$\eta_{s}t\sim O(1)$. Beyond that, one can demonstrate that the equations lead
to asymptotic linear-in-time behavior, with the diffusion coefficient given
above, Eq. (8) [which is in excellent agreement with the behavior in the time
window accessible by exact numerical integration of Eq. (7)].

We observe good agreement between numerics and the Bethe lattice result up to
$\eta_{s}t\sim 10^{1}$. We note that the subsequent departure from quantum
diffusion, which is fitted by an anomalous exponent $t^{0.98}$ and unlikely to
be significant, is much less pronounced than in 2D [22]. This is in
qualitative agreement with the results presented in Sec. III, namely, the fact
that the 3D system is generally less localized than its 2D counterpart.

One should bear in mind that, in general, the effective tight-binding picture
derives from perturbation theory; hence, there is a time limit to the validity
of our model. This is given by $t\lesssim\Delta_{s}/\eta_{s}^{2}$ when next-
nearest neighbor hopping effects are expected to become relevant [22]. After
this timescale, the chirality of the system need not be obeyed, and the long-
timescale behavior could deviate from typical expectations.

##  V  Spinons and visons

Next, we look at the long time regime, where we adopt as in 2D [23] a Born-
Oppenheimer perspective, with the spinons instantaneously in equilibrium on
the timescales over which visons relax (via a stochastic MC process). Given
the system sizes we can access, we limit ourselves again to the case of a
single spinon in Eq. (3).

In our simulations, we use a cubic lattice of size $L^{3}=16^{3}$ with
periodic boundary conditions. We initialize the system in a random vison
background of density $\langle n_{p}\rangle=\frac{1}{2}$ (subject to the gauge
conditions discussed in Sec. II). We then run MC updates using a Metropolis
algorithm as follows: (i) We pick one bond $i$ at random. There are four
plaquettes which contain $i$. We propose an update which flips the vison
number of those four plaquette $n_{p}\to 1-n_{p}$. (ii) A new effective spinon
Hamiltonian $H^{\prime}$, following Eq. (3), is constructed by changing
$A_{ss^{\prime}}$ to $A_{ss^{\prime}}+\pi$ for the bond $i\equiv ss^{\prime}$
that has been selected at point (i). (iii) We accept the proposed change
according to the Metropolis probability
$\textrm{min}[1,\textrm{tr}\>\exp(-\beta H^{\prime})/\exp(-\beta H)]$, where
$H$ is the spinon Hamiltonian before the update, and $\beta=1/T$ is the
inverse temperature (working in energy units where the Boltzmann constant is
$k_{\mathrm{B}}=1$).

We are interested in a regime where the vison energy cost $\Delta_{v}$ does
not play a significant role and we set it to zero for convenience [23].
Moreover, at the temperatures of interest in this paper ($T/\eta_{s}\leq
0.01$), we find that it is sufficient to consider only the six lowest
eigenstates when computing the traces above, and the contribution of higher
excited states can be safely neglected. This simplification greatly enhances
the performance of our code.

Despite the presence of a mobility edge and delocalized states in the spectrum
(see Sec. III), we clearly observe the formation of a vison depletion region
around the spinon wave function and a rapid decay of the latter outside the
depleted region. Correspondingly, most of the results appear to be like the 2D
case [23]. At low temperature, we observe the appearance of a spherical region
where the vison density vanishes, and the spinon wave function is finite
almost exclusively inside this region. Outside, the spinon density vanishes,
and the vison density rapidly approaches the random uncorrelated value
$\frac{1}{2}$. This is illustrated in Fig. 6, where a specific instance is
shown of a vison configuration and corresponding spinon density.

(a)

(b)

(c)

Figure 6:  Vison configuration (blue squares) overlapped with the spinon
probability density (yellow spheres, radius proportional to probability) at
$T_{1}/\eta_{s}=1.0\times 10^{-4}$ (top left panel). The depletion region has
been shifted to the center of the system for visualization convenience, taking
advantage of the periodic boundary conditions. The top right panel shows the
vison density for each layer of the $L^{3}=16^{3}$ lattice, as a function of
layer index for each of the three coordinate axes; this illustrates
quantitatively the presence of a spherical vison depletion region. The bottom
panels are the connected vison correlators $C_{\rho}(r)$, Eq. (9), at the same
temperature: left panel from Monte Carlo simulations; right panel from a
phenomenological model which assumes a vison-depleted sphere surrounded by a
randomly half-filled vison background.

(The spinon density was approximated here for simplicity by the square modulus
of the ground state spinon wave function since it carries most of the
Boltzmann weight; indeed, we find that $\Delta_{01}/\eta_{s}~{}O(10^{-1})\gg
T_{1}/\eta_{s}$, where $\Delta_{01}$ is the spinon energy gap above the ground
state.) We further confirm this behavior statistically by computing the vison
connected correlation function:

$C_{\rho}(s,s^{\prime})=\frac{1}{3}\sum_{\alpha=x,y,z}\left[\langle
n_{p_{\alpha}}n_{p^{{}^{\prime}}_{\alpha}}\rangle-\langle
n_{p_{\alpha}}\rangle\langle n_{p^{{}^{\prime}}_{\alpha}}\rangle\right]\,,$
(9)

also shown in Fig. 6. Here, $p_{\alpha}$ labels the plaquette adjacent to site
$s$ that lies perpendicular to the direction $\alpha$; similarly,
$p^{{}^{\prime}}_{\alpha}$ labels the plaquette adjacent to site $s^{\prime}$
that lies perpendicular to the direction $\alpha$ (where we have uniquely
labeled all plaquettes on the cubic lattice by associating one plaquette for
each of the three directions to one given site of the lattice). The angular
brackets represent the thermodynamic equilibrium average computed in our MC
simulations.

###  V.1  Effective free energy

As for the 2D case [23], the behavior of the system can be understood using a
simple effective free energy model, where the spinon is confined to a
spherical region of radius $\xi$ devoid of visons, and the visons outside the
region take on a random configuration of density $\langle
n_{p}\rangle=\frac{1}{2}$. In this limit, the spinon energy levels are known
analytically [we consider for simplicity only the ground state
$E(\xi)=\pi^{2}\eta_{s}/\xi^{2}$], and the visons contribute only entropically
to the free energy of the system. To account in part for the decay of the
spinon wave function into the region where visons are present (due to a short
but finite localization length), an additional temperature-independent
penetration length scale $\xi_{0}$ is introduced: $E(\xi)\to E(\xi+\xi_{0})$.
The effective free energy $F(\xi)$ can be written as

$F(\xi)=\frac{\pi^{2}\eta_{s}}{(\xi+\xi_{0})^{2}}+\frac{4}{3}\pi{\xi}^{3}T\mathrm{ln}\,[1+\exp(-\beta{{{\Delta}}}_{\text{v}})]\,.$
(10)

Minimization with respect to $\xi$ gives the saddle-point equilibrium value
$\xi_{*}$, dependent on $T$ and $\xi_{0}$:

$\left[1+\left(\frac{\xi_{0}}{\xi_{*}}\right)^{3}\right]\xi_{*}^{5}=\frac{\pi}{2\ln
2}\frac{\eta_{s}}{T}\,,$ (11)

where we set for simplicity $\Delta_{v}=0$. This implicit equation can be
solved numerically, and in the $\xi_{*}\gg\xi_{0}$ regime, it gives the
scaling behavior $\xi_{*}\sim T^{-1/5}$, in contrast to the 2D scaling
$\xi\sim T^{-1/4}$ [23].

To account for thermal fluctuations, one can compute instead the equilibrium
value of $\xi$ using the formula:

$\langle\xi\rangle=\frac{1}{Z}\int_{0}^{R}{d}\xi\ \xi\exp[-\beta
F(\xi)]\,,\quad Z=\int_{0}^{R}{d}\xi\ \exp[-\beta F(\xi)]\,.$ (12)

Here, $R$ is a cutoff that captures the finite-sized effect in the MC
simulations. At the temperatures of interest in this paper, we observe that
the Boltzmann distribution for $\xi$ is in fact highly peaked, and the saddle-
point value $\xi_{*}$ is a very good approximation for $\langle\xi\rangle$.

To determine $\xi_{0}$, we compare the analytical result
$E(\xi)=\pi^{2}\eta_{s}/(\xi+\xi_{0})^{2}$ with simulations of a toy system
where the visons are artificially set to random configurations of probability:

$p({n}_{p})=\left\\{\begin{array}[]{ll}0,&\,\text{if}\,|{{\bf{r}}}_{p}|\,<\,\xi\,,\\\
\frac{1}{2},&\,\text{otherwise}\,.\end{array}\right.$ (13)

The ground state energy of the spinon averaged over $100$ independent vison
configurations is shown in the inset of Fig. 7, demonstrating good agreement
for $\xi_{0}=1.814(0)$.

Equipped with the best fit value of $\xi_{0}$, we can then compare the
equilibrium behavior predicted by the effective model in Eqs. (10 and (12)
with the MC simulations of the spinon and vison system introduced at the start
of this section. This is done in Fig. 7 for the vison density and the radius
of the vison depletion region.

Figure 7:  Dependence of the average vison density $\langle n_{p}\rangle$ per
plaquette (blue circles, left vertical scale) and of the radius of the
depleted region $\sqrt[3]{\xi^{3}}$ (red triangles, right vertical scale) as a
function of temperature. The dashed black line is the analytical result from
the effective model in the main text, Eq. (10). The vertical dotted lines
indicate the temperatures at which finite-sized instabilities occur in the
shape of the depleted region, from spherical to cylindrical $T_{\mathrm{sc}}$
and from cylindrical to planar slab $T_{\mathrm{csl}}$ as discussed in
Appendix C. $T_{1}$ is the temperature used in Fig. 6, as a reference. The
simulations were performed on an $L^{3}=16^{3}$ lattice with periodic boundary
conditions, averaged over $50$ independent runs. The inset shows fitting the
ground state energy of a spinon in an artificial vison background generated
according to Eq. (13) on a $L^{3}=20^{3}$ lattice, as a function of $\xi$ and
averaged over $500$ realizations (squares), to the energy
$E(\xi)=\pi^{2}\eta_{s}/(\xi+\xi_{0})^{2}$ discussed in the main text (solid
line) to obtain $\xi_{0}=1.814(0)$. The vertical axis has been shifted and
rescaled for convenience.

We show $22$ temperature points from $\textrm{log}_{10}T/\eta_{s}=-4.625$
(where the vison density reaches the $1/{\rm volume}$ limit for our finite
system size) to $\textrm{log}_{10}T/\eta_{s}=-2$, each averaged over $50$
independent cooling histories. The effective free energy model appears to
capture the behavior of the system decently well.

As temperature is lowered, the size of the vison depletion region increases.
In a finite-sized system, it eventually percolates, and in doing so, we
observe that it deforms first into a cylinder and then into a planar slab.
This is a finite-sized effect that is largely outside the scope of this paper;
however, for completeness, we present some results about it in Appendix C. The
temperatures at which these two instabilities occur, according to our
simulations, are $T_{\mathrm{sc}}/\eta_{s}=6.1\times 10^{-5}$ (sphere to
cylinder) and $T_{\mathrm{csl}}/\eta_{s}=4.5\times 10^{-5}$ (cylinder to
planar slab) and are indicated for reference by vertical dotted lines in Fig.
7.

###  V.2  Mutual statistics vs interactions

The behavior discussed so far is entirely driven by the mutual statistics
between the noninteracting quasiparticles. The spinon localization that
underpins it, however, could in principle be driven also by possible spinon-
vison interactions. Therefore, it is important to contrast the two cases
explicitly.

For concreteness, we consider the case of contact interactions of strength
$U=\eta_{s}/4$ added to the spinon tight-binding model above:

$H_{\textrm{int}}=H_{\textrm{eff}}+U\sum_{s}\left(\frac{1}{12}\sum_{s\in
p}n_{p}\right)b_{s}^{\dagger}b_{s}\,.$ (14)

Here, $H_{\textrm{eff}}$ is the effective Hamiltonian in Eq. (3), and $s\in p$
indicates that the sum runs over the plaquettes $p$ which are adjacent to the
spinon lattice site $s$.

We then proceed to compare the behavior of the new system with
($A_{ss^{\prime}}\neq 0$, $\pi$ flux) and without ($A_{ss^{\prime}}=0$, zero
flux) mutual statistics between the spinons and the visons. In the latter
case, $H_{\rm int}$ describes a 3D tight-binding system with uniform hopping
strength $\eta_{s}$ and diagonal disordered potential, which is generally
expected to contribute to the spinon localization.

In the regime of interest to us ($U<\eta_{s}$, as we expect any interactions —
if at all present — to be weak perturbations of the original model),
localization due to diagonal disorder occurs over long length scales. This is
in stark contrast with the very pronounced localization effects due to
semionic statistics (equivalently, $\pi$ fluxes). Correspondingly, in our
simulations, we do not observe a well-defined vison depletion region in the
purely interacting case; the vison density is suppressed due to the presence
of spinons in the system, but this occurs in a uniform way that is in fact
well captured by a simple modeling with a tunable uniform chemical potential.
Figure 8 illustrates the dependence of the vison density as a function of
temperature for the two cases at hand.

(a)

(b)

Figure 8:  Vison density as a function of temperature in the presence of
spinon-vison contact interactions of strength $U=\eta_{s}/4$, with (red
triangles) and without (blue circles) a mutual statistical angle ($\phi=\pi$
vs $\phi=0$). The behavior when there are no mutual statistics between spinons
and visons is well captured by a trivial change in the overall chemical
potential for the visons ($\mu=2.134\times 10^{-5}\eta_{s}$ gives a good fit),
in stark contrast with the dependence observed in the original model.
Correspondingly, the vison density decreases in an essentially uniform way
across the system, instead of forming a spherical depleted region. This is
illustrated by the behavior of the connected correlator $C_{\rho}(r)$ in the
bottom panels, for the case with (left) and without (right) mutual statistics.
The correlators are computed at $T_{2}/\eta_{s}=3.16\times 10^{-4}$. Each data
point is averaged over $32$ separate histories.

Once again, as it was the case in 2D [22, 23], we find that the effects of
mutual statistics are robust and more pronounced than (reasonable) spinon-
vison interactions, suggesting that evidence in candidate $\mathbb{Z}_{2}$
spin liquid systems of spinon localization in vison-depleted regions would be
a signature of anyonic statistics, fractionalization and of a topologically
ordered QSL.

##  VI  Response out-of-equilibrium

Finally, we briefly discuss some of the effects that the phenomenology
encountered in our system has on its behavior out of equilibrium. Following
the 2D case in Ref. 23, we consider sweeps where the temperature is varied at
a finite rate, starting from an equilibrium state at $T_{b}$ where the spinon
density $\rho_{s}\sim\exp(-\beta\Delta_{s})$ is much smaller than $1/\xi^{3}$
(as schematically illustrated in Fig. 9).

Figure 9:  Schematic illustration of the spinon density out of equilibrium.
The system is initially at $T_{b}$ (see main text). Upon lowering the
temperature, the spinon density $\rho_{s}$ is unable to relax due to the slow
diffusion of the spinons dressed by their vison depletion regions. The latter,
however, grow in size as temperature is lowered and eventually come into
contact with one another at $T_{c}$, when the correlation radius $\xi^{*}$
becomes comparable with the spinon separation distance $\sim\rho_{s}^{-1/3}$.
This produces a kinematic locking of the spinon density to the
$1/(\xi^{*})^{3}\sim T^{3/5}$ curve. When the thermal sweep is reversed at
$T_{d}$ and the temperature is increased, once again the spinon density enters
a plateau as annihilation events stop, until thermal equilibrium is restored
$\rho_{\rm eq}\sim\exp(-\Delta_{s}/T)$.

Upon lowering $T$, the system ought to lower its spinon density to remain in
equilibrium. Given the nature of the initial state, spinon annihilation
processes (which can only occur in pairs) require them to travel cooperatively
with their own vison depletion region for large distances before they meet
another spinon. Like the 2D case, this can be expected to be a slow process,
and if $T$ is lowered at a sufficiently fast rate, the system then falls
strongly out of equilibrium at some constant spinon density.

While the spinon motion is expected to be slow (cooperative), the size of the
vison depletion region can easily adjust to remain in thermal equilibrium
(driven by vison creation/annihilation events at the edge of the region). As
temperature is lowered, one therefore expects the depletion regions to grow
around each spinon — thus ultimately and necessarily curtailing the constant
spinon density regime. Once the depletion regions, and therefore the
corresponding spinon wave functions, become large enough to overlap, spinon
annihilation processes resume. In the effective free energy picture, this
occurs when

$\rho_{0}\,\xi_{*}^{3}\sim
1\>\leftrightarrow\>T\sim\eta_{s}\,\rho_{0}^{5/3}\,,$ (15)

where $\rho_{0}$ is the initial spinon density. The right-hand side of the
equation is obtained using the scaling relation $\xi_{*}\sim T^{-1/5}$ derived
in the large $\xi_{*}$ limit in Sec. V.1.

As we continue to lower the temperature, the spinon density enters a new
kinematically locked regime where it decreases at the rate set by the
condition that the typical depletion region diameter equals the mean
separation between spinons (see Fig. 9). Using again the left-hand side of the
expression in Eq. (15), one finds the kinematically locked spinon density:

$\rho_{\rm kl}\sim\xi_{*}^{-3}\sim T^{3/5}\,.$ (16)

This differs quantitatively from the 2D case [23], where it was found that
$\rho_{\rm kl}(T)\sim T^{1/2}$.

When the thermal sweep direction is reversed, the vison depletion regions
contract in size as temperature rises, and spinon annihilation stops. This
generates a new plateau in the spinon density at some $\rho_{s}$ much greater
than the thermal equilibrium value. This plateau lasts until the temperature
becomes high enough to restore equilibrium $\rho_{s}\sim\exp(-\Delta_{s}/T)$,
at which point the spinon density resumes increasing following the adiabatic
curve, completing the hysteresis loop (as illustrated in Fig. 9).

##  VII  Discussion and conclusions

In this paper, we studied the behavior of a 3D $\mathbb{Z}_{2}$ lattice gauge
theory (toric code) [39, 40], in the finite-temperature regime where the
spinon excitations are sparse and can hop coherently across the system,
whereas the vison excitations are quasistatic and stochastic. The latter
therefore act as $\mathbb{Z}_{2}$ ($\pi$) fluxes for the former via their
mutual semionic statistics. This model exhibits an instance of Anderson
localization [32] driven by anyonic statistics, where both the noninteracting
tight-binding particles as well as the disorder are borne out of the same spin
degrees of freedom in an otherwise disorder-free topological QSL.

We investigated the localization properties of the system and found that it
exhibits a mobility edge, which is pushed toward the end of the spectrum. We
determined the corresponding transition energy using various standard
techniques such as level statistics [79, 57], $r$-statistics [69, 70], and
IPR/multifractal coefficients [38]. The vicinity of the mobility edge to the
tail end of the spectrum leaves only a small proportion of states in the
localized region at the extreme ends of the spectrum. By using finite-sized
scaling, we provide evidence that the system belongs to the 3D Anderson
universality class [41, 42], and therefore, states inside the mobility edges
are fully extended and delocalized.

Furthermore, we inspect the out-of-equilibrium dynamics. As expected, we found
that, at long times, the system shows diffusion as the 3D Anderson model with
onsite disorder [74, 75, 76, 77, 78]. This result should be compared with the
2D case, where a strong subdiffusive-like behavior was observed [22]. Finally,
we support our numerical evidence with an analytical computation using a Bethe
lattice approximation.

Interestingly, once the visons are allowed to relax stochastically in the
presence of a spinon, the behavior observed in 3D remains remarkably close to
the one encountered in 2D [22], despite the presence of a mobility edge and
delocalized states. We observe that visons are expelled from the high spinon
density area, forming a vison-depleted region determined by the balance
between spinon kinetic energy and vison configurational entropy. The effective
energy stipulates that the size of the patches have a typical radius that
scales with temperature as $T^{-1/5}$. The close resemblance between 2D and 3D
behavior despite a mobility edge is probably due to the fact that low-energy
states close to the lower edge of the band dominate in importance in the
thermodynamic behavior that leads to the appearance of the vison depletion
region –– hence, the localized behavior dominates. We note, however, a longer
penetration depth of the spinons into the random vison background compared
with 2D, hinting at weaker localization in 3D.

As in the 2D case, such nontrivial interactions between spinons and visons are
driven by mutual statistics. We have indeed shown that the spinon-vison
coupling by itself is not enough to induce clearly visible vison-depleted
regions. In this paper, we ignored possible interaction terms between the same
quasiparticle species, which may however be present in real materials. While
the topological spin liquid behavior is protected against sufficiently weak
interactions by the presence of a finite gap, the localization properties are
subject to change in what would then be a many-body localized
system.Nonetheless, the qualitative picture of the physics we observe, up to
the intermediate length and timescales considered in this paper, is not likely
to change. Indeed, it is reasonable to expect that vison-vison interactions
only alter the entropic term of the effective free energy quantitatively but
not qualitatively, whereas spinons, in the low-temperature low-density regime,
are separated by length scales much longer than any reasonable characteristic
interaction range, and their behavior should remain unchanged up to reasonably
long timescales.

We also briefly discussed the nonequilibrium behavior of the system when the
temperature is varied in time at some finite rate. We expect a similar
hysteretic behavior as in 2D [23]. When the system is cooled, the spinon
density $\rho_{s}$ falls off equilibrium and forms a plateau, before it
eventually becomes kinematically locked to $\rho_{s}\sim T^{3/5}$ when the
vison-depleted regions become large enough to overlap (i.e., their size is
approximately equal to the characteristic spinon separation length).

Our results demonstrate how the rich phenomenology at the interface between
topological order and Anderson localization in disorder-free systems arises
also in higher dimensional systems. Whereas the localization properties are
different, the notable effects on the relaxation and transport properties
underpinned by the formation of vison depletion regions around the spinons
survive. This behavior is driven by the mutual semionic statistics and can
thus be taken as a signature of 3D QSL behavior at finite temperature. The
hysteretic behavior discussed in Sec. VI could be probed, for instance, with
techniques that access spinon density, which is expected, for example, to
directly affect the magnetic susceptibility [23]. The presence of vison
depletion regions bound to the spinons is also likely to have distinctive
repercussions on transport properties where spinons or visons contribute
(e.g., thermal transport). An interesting future direction could be to model
more extensively the dynamical interplay of spinons and visons, including
spinon annihilation events when vison depletion regions overlap, and possibly
simulating the out-of-equilibrium behavior directly. Another direction could
be to extend our analysis to other classes of Hamiltonians such as
$\mathbb{Z}_{n}$ gauge theories where the effective fluxes break time-reversal
symmetry, and to general 3D string-net models [9, 39]. As we await the
discovery of candidate materials that realize $\mathbb{Z}_{2}$ QSL phases in
3D, our results may be relevant to other contexts, including frustrated
magnetic pyrochlore oxides and resonant valence bond systems. Some of these
systems exhibit further gapless excitations, and it would be interesting to
see how the behavior discussed in this paper is affected by their presence.
Moreover, the possibility of realizing $\mathbb{Z}_{2}$ spin liquid
Hamiltonians in our temperature regime with quantum annealers [24, 44] and
quantum simulators [45, 46] (albeit typically limited to 2D) could provide a
suitable arena where the physics discussed here could be tested and explored
further.

###### Acknowledgements.

We are very grateful to O. Hart for the generous guidance in the early stages
of this project and for several useful discussions thereafter. This paper was
supported in part by the Engineering and Physical Science Research Council
(EPSRC) Grants No. EP/P034616/1, No. EP/T028580/1, and No. EP/V062654/1 (CC).
GDT acknowledges support from the EPiQS Program of the Gordon and Betty Moore
Foundation. MK developed and performed the calculations and numerical
simulations.

##  Appendix A  Further properties of localized and extended states

In this appendix, we present details of the study of the DOS and the scaling
properties of the IPR for the localized and extended states, and we further
discuss the behavior of the probability distributions and level statistics in
regions around the mobility edge to complement our discussion in Sec. III.

In the main text, we determined the mobility edge to be around
$E/\eta_{s}\sim-4.35$, which is close to the band edge, as shown in Fig. 10.
This leads to a suppression of the DOS in the localized region of the spectrum
and a relatively quickly varying DOS near the mobility edge. The smaller DOS
also leads to larger Mott temperature, making finite-temperature dephasing
effects relevant only at long timescales, as discussed in the main text.

Figure 10: Density of states (DOS) for various system sizes. The mobility
edge, being close to the band edge, is on the region where the DOS profile is
not flat. We note the generally higher DOS for larger system sizes closer to
the band tail, which is compensated by the relative depression near the band
center. The dashed vertical line indicates the mobility-edge and is a guide
for eye.

In Fig. 11, we show directly the scaling of IPRq, contrasting states in the
localized vs delocalized portions of the energy spectrum.

Figure 11:  Scaling with system size of the IPRq for localized (top panels)
and delocalized (bottom panels) eigenstates, for $q=2$ (left panels) and
$q=\frac{1}{2}$ (right panels). For localized states, $\ln{\rm IPR}_{q}$ are
sampled from states in the energy window $E/\eta_{s}\in[-4.65,-4.55]$ in units
of the spinon hopping. The extended states are sampled in the window
$E/\eta_{s}\in[0.95,1.05]$.

The localized states are consistent with IPRq taking a constant value (the
fitted slope being close to zero), whereas the delocalized states are in very
good agreement with the expected scaling, $\ln\langle{\rm IPR}_{q}\rangle\sim
3(1-q)\ln L$.

Near the mobility edge, the level statistics and the probability distributions
of IPR2 and $D_{2}$ gradually change from localized to extended behavior as we
sample higher energy states. This is illustrated in Fig. 12 for a range
energies from localized and extended.

Figure 12:  Probability distributions for the multifractal coefficient $D_{2}$
(left panels) and the IPR2 (insets), alongside the level statistics (right
panels), for a fixed energy window of width $0.1\eta_{s}$ placed in different
regions of the spectrum, spanning localized and extended states. From top to
bottom, the values of the energy in the middle of each window are
$E/\eta_{s}=-4.5$, $-4.4$, $-4.3$, respectively. We note a semi-Poissonian
distribution appearing just before the mobility edge is reached, due to the
fact that the probability distributions of the IPR do not exactly align with
each other at energies lower but close to the mobility edge (e.g., near
$E/\eta_{s}\sim-4.4$).

As the energy increases from the lowest end of the spectrum, the level
statistics transitions from a Poissonian distribution to a Wigner surmise, via
a semi-Poissonian regime [79, 57]. The transition starts from around
$E/\eta_{s}\sim-4.5$, which indicates the existence of level repulsion even in
regions slightly below the mobility edge. The emergence of the semi-Poissonian
is due to the fact that the states are not fully localized and it is the
typical behavior around the Anderson transition [80].

Correspondingly, around this energy, the tail ends of the probability
distributions of the IPR for different system sizes start to deviate from each
other. The behavior of the multifractal coefficient distribution is like IPR2,
given that the shape of the distribution curves starts to transition at
similar energies. These results are consistent with the $r$ statistics (see
also Fig. 3).

Figure 13 shows the probability distribution results for $D_{\frac{1}{2}}$ and
IPR${}_{\frac{1}{2}}$, which exhibit the same behavior as $q=2$ above.

Figure 13:  Probability distributions for the multifractal coefficient and the
IPR for $q=\frac{1}{2}$, for a range of energies between $E/\eta_{s}=-4.55$ to
$-4.25$, with energy spacing $0.1\eta_{s}$ and energy window width of
$0.1\eta_{s}$. The energies increase from left to right and from top to
bottom. We see an overall similar pattern with $D_{2}$ and IPR2.

##  Appendix B  Bethe lattice calculation

Following the Bethe lattice calculation in Ref. 22, we derive similar results
for our model defined on a 3D cubic lattice. We start from the generating
functions for closed and open walks on a Bethe lattice of coordination number
$z=6$ (or branching ratio $z-1$), which are given as

$\displaystyle T^{(z)}(x)$ $\displaystyle=$
$\displaystyle\frac{2(z-1)}{z-2+z\sqrt{1-4(z-1)x^{2}}},$ (17) $\displaystyle
S^{(z)}(x)$ $\displaystyle=$
$\displaystyle\frac{1-\sqrt{1-4(z-1)x^{2}}}{2(z-1)x}\,.$ (18)

We can also directly use the analytical form for the $2k$th moment of the
displacement:

$\displaystyle\left\langle\mathbf{r}^{2k}(t)\right\rangle$ $\displaystyle=$
$\displaystyle\oiint\frac{dw_{1}}{2\pi i}\frac{dw_{2}}{2\pi
i}\frac{\exp[iht\left(w_{1}-w_{2}\right)]}{w_{1}w_{2}}T\left(\frac{1}{w_{1}}\right)T\left(\frac{1}{w_{2}}\right)$
(19) $\displaystyle\times$
$\displaystyle\mathcal{R}_{2k}\left[S\left(\frac{1}{w_{1}}\right)S\left(\frac{1}{w_{2}}\right)\right]\,.$

The nonreversing constraints, however, need to be modified. Unlike the 2D
$z=6$ case, the constraints now have three variables, as 3D lattice walks
require three independent vectors. We introduce the generating variables
$x,\delta,\epsilon$, and $\zeta$ which count the length of the walk and the
number of steps taken in the lattice directions
$\mathbf{e}_{1},\mathbf{e}_{2}$, and $\mathbf{e}_{3}$, respectively. At each
lattice point, there are six possible walks:
$\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}$ and their reverses, each noted
by $\delta,\epsilon,\zeta,\delta^{-1},\epsilon^{-1}$, and $\zeta^{-1}$,
respectively, in generating variables. In addition to the first step, however,
the immediate reversing of the previous step is forbidden. We can denote this
in a matrix, given as

$N=x\left(\begin{array}[]{cccccc}\delta&\epsilon&\zeta&\zeta^{-1}&\epsilon^{-1}&0\\\
\delta&\epsilon&\zeta&\zeta^{-1}&0&\delta^{-1}\\\
\delta&\epsilon&\zeta&0&\epsilon^{-1}&\delta^{-1}\\\
\delta&\epsilon&0&\zeta^{-1}&\epsilon^{-1}&\delta^{-1}\\\
\delta&0&\zeta&\zeta^{-1}&\epsilon^{-1}&\delta^{-1}\\\
0&\epsilon&\zeta&\zeta^{-1}&\epsilon^{-1}&\delta^{-1}\end{array}\right),$ (20)

with the initial condition
$N_{0}=x\operatorname{diag}\left(\delta,\epsilon,\zeta,\zeta^{-1},\epsilon^{-1},\delta^{-1}\right)$.
We note that, for each step, the length of the path is advanced by one, and
the matrix multiplication ensures that all possible combinations of steps are
accounted for. The zero entries represent the nonreversing constraint, while
$N_{0}$ is unconstrained.

The matrix elements $N_{0}N^{l}$ therefore give the paths of length $l$ that
are consistent with the nonreversing constraint. The generation function
$\mathcal{N}(x;\delta,\epsilon,\zeta)$ for nonreversing paths can then be
defined as the sum over all matrix elements and all possible path lengths $l$,
including also the trivial walk of zero length. This can be written as

$\begin{split}\mathcal{N}(x;\delta,\epsilon,\zeta)&=1+\sum_{i,j}\sum_{\ell=1}^{\infty}\left[N_{0}N^{\ell-1}\right]_{ij}\\\
&=1+\sum_{i,j}\left[N_{0}\left(\mathbbm{1}_{z}-N\right)^{-1}\right]_{ij}\\\
&=\frac{1-x^{2}}{1+5x^{2}-x\left(\delta+\delta^{-1}+\epsilon+\epsilon^{-1}+\zeta+\zeta^{-1}\right)}\,.\end{split}$
(21)

The generator for the $2k$th moment of the displacement therefore is

$\begin{split}\mathcal{R}_{2k}(x)&\equiv\sum_{s\in\mathcal{L}_{z}}\left(s_{x}^{2}+s_{y}^{2}+s_{z}^{2}\right)^{k}C_{s}(x)\\\
&=\left.\left\\{\left[\left(\delta\partial_{\delta}\right)^{2}+\left(\epsilon\partial_{\epsilon}\right)^{2}+(\zeta\partial_{\zeta})^{2}\right]^{k}\mathcal{N}(x;\delta,\epsilon,\zeta)\right\\}\right|_{\delta=\epsilon=\zeta=1}\,.\end{split}$
(22)

For $k=1$, the generator takes the simpler form:

$\mathcal{R}_{2}(x)=\frac{6x(1+x)}{(1-5x)^{2}(1-x)}\,.$ (23)

When we substitute the relevant generators into Eq. (19), we notice that there
is a line of poles due to
$S\left(u+i0^{+}\right)S\left(u-i0^{+}\right)=(z-1)^{-1}$ for $u\in\mathbb{R}$
and $|u|<2\sqrt{z-1}$. The pole signature and positions are the same as in 2D;
hence, we can directly use

$2D_{z}^{\mathrm{diff}}=\frac{1}{2\pi}\mathcal{F}(z)\lim_{w\rightarrow(z-1)^{-1}}\left[1-w(z-1)\right]^{2}\mathcal{R}_{2}^{(z)}(w)\,,$
(24)

where

$\begin{split}\mathcal{F}(z)&\equiv(z-1)\int_{-2\sqrt{z-1}}^{2\sqrt{z-1}}~{}du\frac{4(z-1)-u^{2}}{z^{2}-u^{2}}\\\
&=(z-1)\left[4\sqrt{z-1}-z\left(\frac{z-2}{z}\right)^{2}\ln\left(\frac{z+2\sqrt{z-1}}{z-2\sqrt{z-1}}\right)\right],\end{split}$
(25)

which is the quantity we obtain after integrating over the relevant residue
line.

The final asymptotic diffusion coefficient is given by

$D_{6}^{\mathrm{diff}}=\frac{3}{\pi}\left[3\sqrt{5}-2\ln\left(\frac{3+\sqrt{5}}{3-\sqrt{5}}\right)\right]\simeq
2.72968\,,$ (26)

which is coincidentally identical to the case of a triangular lattice in 2D.

##  Appendix C  Low temperature finite-size instabilities of the vison
depletion region

At the lowest temperatures, the vison depletion region discussed in Sec. V
becomes comparable in size with the simulation cell. The system displays
finite-sized instabilities which alter the shape of the depletion region from
a sphere to a cylinder and eventually to a planar slab (see Fig. 14) as the
temperature is lowered.

(a)  (b)

(c)

Figure 14:  Snapshots of the spinon probability density (yellow spheres) and
the corresponding vison background (in blue) from our Monte Carlo simulations.
The spinon wavefunction and depletion region are generally spherical (top left
panel). At low temperatures, due to finite-sized effects, we observe an
instability first to a cylindrical shape (top right panel) and finally to a
planar slab (bottom panel) before the visons are completely expelled from the
system. The temperatures for these visualizations are $T/\eta_{s}=1.0\times
10^{-4}$, $5.2\times 10^{-5}$, and $4.0\times 10^{-5}$ for the sphere,
cylinder, and slab, respectively.

This is observed in the numerics and can be understood by effective free
energy considerations like those of Eq. (10). Here, we only provide an
analytical calculation of the transition temperatures at saddle-point
approximation and in the limit $\xi\gg\xi_{0}$ and $\Delta_{v}=0$, for
convenience. In the case of a spherical region, we minimize Eq. (10) with
respect to $\xi$ to obtain

$\displaystyle\xi_{*}^{s}$ $\displaystyle=$
$\displaystyle\left(\frac{\pi\eta_{s}}{2T\ln 2}\right)^{1/5},$ (27)
$\displaystyle F(\xi_{*}^{s})$ $\displaystyle=$
$\displaystyle\frac{5}{3}\pi^{2}\left(\frac{2T\ln
2}{\pi\eta_{s}}\right)^{2/5}\,.$ (28)

In the case of a cylindrical region, Eq. (10) needs to be modified to

$F_{c}(\xi)=\frac{j_{0}^{2}\eta_{s}}{(\xi+\xi_{0})^{2}}+\pi{\xi}^{2}LT\mathrm{ln}\,[1+\exp(-\beta{{{\Delta}}}_{\text{v}})]\,,$
(29)

shown for completeness including $\xi_{0}$ and $\Delta_{v}$. Minimizing again
with respect to $\xi$ gives

$\displaystyle\xi_{*}^{c}$ $\displaystyle=$
$\displaystyle\left(\frac{j_{0}^{2}\eta_{s}}{\pi TL\ln 2}\right)^{1/4},$ (30)
$\displaystyle F_{c}(\xi_{*}^{c})$ $\displaystyle=$ $\displaystyle
2j_{0}^{2}\eta_{s}\left(\frac{\pi LT\ln 2}{j_{0}^{2}\eta_{s}}\right)^{1/2}\,.$
(31)

Finally, in the case of a planar slab depletion region,

$F_{sl}(\xi)=\frac{\pi^{2}\eta_{s}}{4(\xi+\xi_{0})^{2}}+2\xi
L^{2}T\mathrm{ln}\,[1+\exp(-\beta{{{\Delta}}}_{\text{v}})]\,,$ (32)

shown again for completeness including $\xi_{0}$ and $\Delta_{v}$. We minimize
with respect to $\xi$ to obtain

$\displaystyle\xi_{*}^{sl}$ $\displaystyle=$
$\displaystyle\left(\frac{\pi^{2}\eta_{s}}{4TL^{2}\ln 2}\right)^{1/3},$
$\displaystyle F_{sl}(\xi_{*}^{sl})$ $\displaystyle=$
$\displaystyle\frac{3}{4}\pi^{2}\eta_{s}\left(\frac{4L^{2}T\ln
2}{\pi^{2}\eta_{s}}\right)^{2/3}\,.$ (33)

The stable equilibrium geometry of the vison depletion region is determined by
the lowest free energy, and transitions between them occur at the crossing
temperature where Eq. (28) equals Eq. (31):

$T_{\mathrm{sc}}=\left(\frac{5\pi}{6j_{0}}\right)^{10}\frac{16\pi}{\ln
2}\frac{\eta_{s}}{L^{5}}\simeq 1.62\times 10^{-4}\eta_{s}\,,$ (34)

and where Eq. (31) equals Eq. (33):

$T_{\mathrm{csl}}=\left(\frac{2j_{0}}{3\pi^{1/3}}\right)^{6}\frac{16\pi}{\ln
2}\frac{\eta_{s}}{L^{5}}\simeq 1.19\times 10^{-4}\eta_{s}\,,$ (35)

where indeed $T_{\mathrm{sc}}>T_{\mathrm{csl}}$. Both instabilities take place
within the physical range of parameters, namely, well before the vison density
falls below the $1/{\rm Volume}$ limit.

The same calculations can be carried out including a finite $\xi_{0}$ fitting
parameter as in the main text. The equations are less transparent, and for
convenience, we limit ourselves to illustrate the results graphically.
Following the same procedure as in Fig. 7, we fit the $\xi_{0}$ values for the
cylindrical and planar slab geometries in Fig. 15.

Figure 15:  Ground state energy of a spinon in an artificial vison background
generated with depleted cylindrical and planar slab regions spanning the
entire system, as a function of the linear size $\xi$ of the vison depleted
region. The data are averaged over $500$ realizations of systems of size
$L^{3}=20^{3}$. The error bars are shown (and generally smaller than the
symbol size). The dashed lines are the best fits to the energy of an infinite
three-dimensional (3D) potential well of corresponding shape but with linear
size $\xi+\xi_{0}$ (the radius of the cylinder and the thickness of the slab);
this gives $\xi_{0}=1.835(3)$ and $=1.869(2)$ for the cylinder and the slab,
respectively. The vertical axis has been shifted and rescaled for convenience.

With the respective fit values $\xi_{0}=1.814(0),\,1.835(3),\,1.869(2)$
(sphere/cylinder/slab), the improved saddle-point free energies are plotted in
Fig. 16 as a function of temperature, showing the finite-sized instabilities
at the crossing points $T/\eta_{s}=1.50\times 10^{-4}$ (sphere to cylinder)
and $T/\eta_{s}=7.67\times 10^{-5}$ (cylinder to planar slab).

Figure 16:  Effective free energy plots for temperatures between
$T/\eta_{s}=10^{-5}$ and $10^{-3}$ for the three geometries of the vison
depletion region: spherical, cylindrical, and planar slab. Here, we use the
improved free energies that account for the penetration depth $\xi_{0}$
obtained from the fits in Figs. 7 and 15 for system size $L^{3}=16^{3}$. The
free energy diagram predicts a transition from sphere to cylinder at
$T_{\mathrm{sc}}/\eta_{s}=1.50\times 10^{-4}$ and a transition from cylinder
to strip at $T_{\mathrm{csl}}/\eta_{s}=7.67\times 10^{-5}$.

The saddle-point approximation discussed in this appendix paints a picture
that is consistent with the one empirically observed in our MC simulations,
which show a transition from a sphere to a cylinder around
$T_{\mathrm{sc}}/\eta_{s}\simeq 6.1\times 10^{-5}$ (see also Fig. 7) and a
transition from cylinder to planar slab around
$T_{\mathrm{csl}}/\eta_{s}\simeq 4.5\times 10^{-5}$.

## References

  * [1] L. Savary and L. Balents.  Quantum spin liquids: a review.  Rep. Prog. Phys., 80(1):016502, nov 2017.
  * [2] Y. Zhou, K. Kanoda, and T.-K. Ng.  Quantum spin liquid states.  Rev. Mod. Phys., 89:025003, Apr 2017.
  * [3] J. Knolle and R. Moessner.  A field guide to spin liquids.  Annu. Rev. Condens. Matter Phys., 451, 04 2019.
  * [4] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil.  Quantum spin liquids.  Science, 367(6475):eaay0668, 2020.
  * [5] X.-G. Wen.  Quantum Field Theory of Many-Body systems: From the Origin of Sound to an Origin of Light and Electrons.  (Oxford University Press, Oxford, 2007).
  * [6] F. Wilczek.  New Kinds of Quantum Statistics, pages 61–69.  (Birkhäuser Basel, Basel, 2009).
  * [7] R. Moessner and J. E. Moore.  Topological Phases of Matter.  (Cambridge University Press, Cambridge, 2021).
  * [8] A. Kitaev.  Fault-tolerant quantum computation by anyons.  Ann. Phys., 303(1):2–30, 2003.
  * [9] M. A. Levin and X.-G. Wen.  String-net condensation: A physical mechanism for topological phases.  Phys. Rev. B, 71:045110, Jan 2005.
  * [10] 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–1159, Sep 2008.
  * [11] T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee.  Fractionalized excitations in the spin-liquid state of a kagome-lattice antiferromagnet.  Nature, 492(7429):406–410, Dec 2012.
  * [12] J. Nasu, J. Knolle, D. L. Kovrizhin, Y. Motome, and R. Moessner.  Fermionic response from fractionalization in an insulating two-dimensional magnet.  Nat. Phys., 12(10):912–915, Oct 2016.
  * [13] J. Yoshitake, J. Nasu, and Y. Motome.  Fractional Spin Fluctuations as a Precursor of Quantum Spin Liquids: Majorana Dynamical Mean-Field Study for the Kitaev Model.  Phys. Rev. Lett., 117:157203, Oct 2016.
  * [14] Y. Shen, Y.-D. Li, H. Wo, Y. Li, S. Shen, B. Pan, Q. Wang, H. C. Walker, P. Steffens, M. Boehm, Y. Hao, D. L. Quintero-Castro, L. W. Harriger, M. D. Frontzek, L. Hao, S. Meng, Q. Zhang, G. Chen, and J. Zhao.  Evidence for a spinon fermi surface in a triangular-lattice quantum-spin-liquid candidate.  Nature, 540(7634):559–562, Dec 2016.
  * [15] J. A. Paddison, M. Daum, Z. Dun, G. Ehlers, Y. Liu, M. Stone, H. Zhou, and M. Mourigal.  Continuous excitations of the triangular-lattice quantum spin liquid YbMgGaO4.  Nat. Phys., 13(2):117–122, Feb 2017.
  * [16] S. C. Morampudi, A. M. Turner, F. Pollmann, and F. Wilczek.  Statistics of Fractionalized Excitations through Threshold Spectroscopy.  Phys. Rev. Lett., 118:227201, May 2017.
  * [17] J. Yoshitake, J. Nasu, Y. Kato, and Y. Motome.  Majorana dynamical mean-field study of spin dynamics at finite temperatures in the honeycomb Kitaev model.  Phys. Rev. B, 96:024438, Jul 2017.
  * [18] J. Yoshitake, J. Nasu, and Y. Motome.  Temperature evolution of spin dynamics in two- and three-dimensional Kitaev models: Influence of fluctuating $\mathbb{Z}_{2}$ flux.  Phys. Rev. B, 96:064433, Aug 2017.
  * [19] S.-H. Do, S.-Y. Park, J. Yoshitake, J. Nasu, Y. Motome, Y. Kwon, D. T. Adroja, D. J. Voneshen, K. Kim, T.-H. Jang, J.-H. Park, K.-Y. Choi, and S. Ji.  Majorana fermions in the Kitaev quantum spin system $\alpha$-RuCl3.  Nat. Phys., 13(11):1079–1084, Nov 2017.
  * [20] S. Chatterjee, J. F. Rodriguez-Nieva, and E. Demler.  Diagnosing phases of magnetic insulators via noise magnetometry with spin qubits.  Phys. Rev. B, 99:104425, Mar 2019.
  * [21] G. B. Halász, S. Kourtis, J. Knolle, and N. B. Perkins.  Observing spin fractionalization in the Kitaev spin liquid via temperature evolution of indirect resonant inelastic x-ray scattering.  Phys. Rev. B, 99:184417, May 2019.
  * [22] O. Hart, Y. Wan, and C. Castelnovo.  Coherent propagation of quasiparticles in topological spin liquids at finite temperature.  Phys. Rev. B, 101:064428, Feb 2020.
  * [23] O. Hart, Y. Wan, and C. Castelnovo.  Correlation holes and slow dynamics induced by fractional statistics in gapped quantum spin liquids.  Nat. Commun., 12(1):1459, Mar 2021.
  * [24] C. Chamon, D. Green, and Z.-C. Yang.  Constructing Quantum Spin Liquids Using Combinatorial Gauge Symmetry.  Phys. Rev. Lett., 125:067203, Aug 2020.
  * [25] M. J. P. Gingras and P. A. McClarty.  Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets.  Rep. Prog. Phys., 77(5):056501, may 2014.
  * [26] M. Udagawa and L. Jaubert, editors.  Spin Ice.  (Springer, Cham, 2021).
  * [27] M. Hermanns, I. Kimchi, and J. Knolle.  Physics of the kitaev model: Fractionalization, dynamic correlations, and material connections.  Annu. Rev. Condens. Matter Phys., 9(1):17–33, 2018.
  * [28] R. Moessner and K. S. Raman.  Quantum Dimer Models, pages 437–479.  Springer, Berlin, Heidelberg, 2011.
  * [29] C. Castelnovo and C. Chamon.  Topological order and topological entropy in classical systems.  Phys. Rev. B, 76:174416, Nov 2007.
  * [30] L. Balents, M. P. A. Fisher, and S. M. Girvin.  Fractionalization in an easy-axis kagome antiferromagnet.  Phys. Rev. B, 65:224412, May 2002.
  * [31] M. Hermele, M. P. A. Fisher, and L. Balents.  Pyrochlore photons: The ${U}(1)$ spin liquid in a $s=\frac{1}{2}$ three-dimensional frustrated magnet.  Phys. Rev. B, 69:064404, Feb 2004.
  * [32] P. W. Anderson.  Absence of diffusion in certain random lattices.  Phys. Rev., 109:1492–1505, Mar 1958.
  * [33] M. Schiulaz, A. Silva, and M. Müller.  Dynamics in many-body localized quantum systems without disorder.  Phys. Rev. B, 91:184202, May 2015.
  * [34] N. Y. Yao, C. R. Laumann, J. I. Cirac, M. D. Lukin, and J. E. Moore.  Quasi-Many-Body Localization in Translation-Invariant Systems.  Phys. Rev. Lett., 117:240601, Dec 2016.
  * [35] A. Smith, J. Knolle, D. L. Kovrizhin, and R. Moessner.  Disorder-free localization.  Phys. Rev. Lett., 118:266601, Jun 2017.
  * [36] H. Yarloo, A. Langari, and A. Vaezi.  Anyonic self-induced disorder in a stabilizer code: Quasi many-body localization in a translational invariant model.  Phys. Rev. B, 97:054304, Feb 2018.
  * [37] A. Smith.  Disorder-free Localization.  Doctoral Thesis, University of Cambridge (Springer Nature Switzerland AG, Cham, 2019).
  * [38] F. Evers and A. D. Mirlin.  Anderson transitions.  Rev. Mod. Phys., 80:1355–1417, Oct 2008.
  * [39] A. Hamma, P. Zanardi, and X.-G. Wen.  String and membrane condensation on three-dimensional lattices.  Phys. Rev. B, 72:035307, Jul 2005.
  * [40] C. Castelnovo and C. Chamon.  Topological order in a three-dimensional toric code at finite temperature.  Phys. Rev. B, 78:155120, Oct 2008.
  * [41] K. Slevin and T. Ohtsuki.  Corrections to Scaling at the Anderson Transition.  Phys. Rev. Lett., 82:382–385, Jan 1999.
  * [42] A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. Römer.  Critical Parameters from a Generalized Multifractal Analysis at the Anderson Transition.  Phys. Rev. Lett., 105:046403, Jul 2010.
  * [43] K. Boothby, P. Bunyk, J. Raymond, and A. Roy.  Next-generation topology of D-Wave quantum processors.  arXiv:2003.00133, (2020).
  * [44] S. Zhou, D. Green, E. D. Dahl, and C. Chamon.  Experimental realization of classical $\mathbb{Z}_{2}$ spin liquids in a programmable quantum device.  Phys. Rev. B, 104:L081107, Aug 2021.
  * [45] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin.  Probing topological spin liquids on a programmable quantum simulator.  Science, 374(6572):1242–1247, 2021.
  * [46] L. Lumia, P. Torta, G. B. Mbeng, G. E. Santoro, E. Ercolessi, M. Burrello, and M. M. Wauters.  Two-dimensional $\mathbb{Z}_{2}$ lattice gauge theory on a near-term quantum simulator: Variational quantum optimization, confinement, and topological order.  PRX Quantum, 3:020320, Apr 2022.
  * [47] C. Castelnovo, C. Chamon, and D. Sherrington.  Quantum mechanical and information theoretic view on classical glass transitions.  Phys. Rev. B, 81:184303, May 2010.
  * [48] F. J. Dyson.  The dynamics of a disordered linear chain.  Phys. Rev., 92:1331–1338, Dec 1953.
  * [49] M. Weissmann and N. V. Cohan.  Density of states of a one-dimensional system with off-diagonal disorder.  J. Phys. C: Solid State Phys., 8(9):L145–L147, may 1975.
  * [50] L. Fleishman and D. C. Licciardello.  Fluctuations and localization in one dimension.  J. Phys. C: Solid State Phys., 10(6):L125–L126, mar 1977.
  * [51] T. A. L. Ziman.  Localization with off-diagonal disorder: A qualitative theory.  Phys. Rev. B, 26:7066–7069, Dec 1982.
  * [52] C. M. Soukoulis, I. Webman, G. S. Grest, and E. N. Economou.  Study of electronic states with off-diagonal disorder in two dimensions.  Phys. Rev. B, 26:1838–1841, Aug 1982.
  * [53] P. Cain, R. Römer, and M. Schreiber.  Phase diagram of the three-dimensional Anderson model of localization with random hopping.  arXiv:cond-mat/9908255, 1999.
  * [54] P. W. Brouwer, C. Mudry, and A. Furusaki.  Density of States in Coupled Chains with Off-Diagonal Disorder.  Phys. Rev. Lett., 84:2913–2916, Mar 2000.
  * [55] P. Biswas, P. Cain, R. Römer, and M. Schreiber.  Off-diagonal disorder in the Anderson model of localization.  Phys. Stat. Sol. (b), 218(1):205–209, 2000.
  * [56] S.-J. Xiong and S. N. Evangelou.  Power-law localization in two and three dimensions with off-diagonal disorder.  Phys. Rev. B, 64:113107, Aug 2001.
  * [57] S. N. Evangelou and D. E. Katsanos.  Spectral statistics in chiral-orthogonal disordered systems.  J. Phys. A: Math. Gen., 36(12):3237–3254, mar 2003.
  * [58] S. N. Taraskin and S. R. Elliott.  Disorder-induced zero-energy spectral singularity for random matrices with correlations.  Phys. Rev. B, 65:052201, Jan 2002.
  * [59] B. K. Nikolić.  Statistical properties of eigenstates in three-dimensional mesoscopic systems with off-diagonal or diagonal disorder.  Phys. Rev. B, 64:014203, Jun 2001.
  * [60] A. Eilmes, R. A. Römer, and M. Schreiber.  Exponents of the localization lengths in the bipartite Anderson model with off-diagonal disorder.  Phys. B: Condens. Matter., 296(1):46–51, Feb 2001.
  * [61] The case of random phases is special in one dimension. Indeed, using gauge transformation $b_{s}->\mathop{exp}\nolimits(i\phi_{s})b_{s}$ one can absorb the phase.
  * [62] C. M. Soukoulis and E. N. Economou.  Off-diagonal disorder in one-dimensional systems.  Phys. Rev. B, 24:5698–5702, Nov 1981.
  * [63] Formally, the localization length is infinite.
  * [64] S. Nishino, S. Katsuno, and M. Goda.  Disorder-induced quantum bond percolation.  J. Phys.: Condens. Matter., 21(3):035602, dec 2009.
  * [65] A. Tadjine and C. Delerue.  Anderson localization induced by gauge-invariant bond-sign disorder in square PbSe nanocrystal lattices.  Phys. Rev. B, 98:125412, Sep 2018.
  * [66] Here, $1000$ vison configurations for $L=15$, $500$ for $L=23$, and $100$ vison configurations for larger system sizes.
  * [67] The width of the energy shell has been chosen to be small enough for the change in the density of states to be negligible.
  * [68] T. Kawarabayashi, B. Kramer, and T. Ohtsuki.  Anderson transitions in three-dimensional disordered systems with randomly varying magnetic flux.  Phys. Rev. B, 57:11842–11845, May 1998.
  * [69] V. Oganesyan and D. A. Huse.  Localization of interacting fermions at high temperature.  Phys. Rev. B, 75:155111, Apr 2007.
  * [70] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux.  Distribution of the Ratio of Consecutive Level Spacings in Random Matrix Ensembles.  Phys. Rev. Lett., 110:084101, Feb 2013.
  * [71] N. F. Mott.  Conduction in non-crystalline materials.  Philos. Mag., 19(160):835–852, 1969.
  * [72] We remark that our choice of random vison background respects the emergent gauge structure of the system, discussed in Sec. II.
  * [73] H. De Raedt.  Product formula algorithms for solving the time dependent Schrödinger equation.  Comput. Phys. Rep., 7(1):1–72, 1987.
  * [74] P. Prelovšek.  Numerical simulation of diffusion in a three-dimensional disordered lattice.  Solid State Commun., 31(3):179–182, 1979.
  * [75] P. Prelovšek.  Numerical Study of the Conductivity in the Vicinity of Mobility Edges.  Phys. Rev. Lett., 40:1596–1599, Jun 1978.
  * [76] T. Ohtsuki and T. Kawarabayashi.  Anomalous diffusion at the Anderson transitions.  J. Phys. Soc. Jpn, 66:314–317, 1997.
  * [77] P. Sierant, D. Delande, and J. Zakrzewski.  Thouless Time Analysis of Anderson and Many-Body Localization Transitions.  Phys. Rev. Lett., 124:186601, May 2020.
  * [78] P. Prelovšek and J. Herbrych.  Diffusion in the Anderson model in higher dimensions.  Phys. Rev. B, 103:L241107, Jun 2021.
  * [79] D. Braun, G. Montambaux, and M. Pascaud.  Boundary Conditions at the Mobility Edge.  Phys. Rev. Lett., 81:1062–1065, Aug 1998.
  * [80] B. I. Shklovskii, B. Shapiro, B. R. Sears, P. Lambrianides, and H. B. Shore.  Statistics of spectra of disordered systems near the metal-insulator transition.  Phys. Rev. B, 47:11487–11490, May 1993.