Vison crystal in quantum spin ice on the breathing pyrochlore lattice

Recent excitement in the quantum spin ice community has come from the experimental discovery of pseudospin-$1/2$ breathing pyrochlores, including Ba$_3$Yb$_2$Zn$_5$O$_{11}$, in which inversion symmetry is broken by the `up' and `down' tetrahedra taking different physical sizes. We show here that the often-neglected $J_{z\pm}$ coupling between Kramers ions, in combination with the breathing nature of the lattice, can produce an imaginary ring flip term. This can lead to an unconventional '$U(1)_{\pi/2}$ phase', corresponding to a maximally dense packing of visons on the lattice. Coherent dynamics persists in all phases, together with its emergent QED description, in a manner reminiscent of fragmentation in spinon crystals. We characterize the enlarged QSI phase diagram and its excitations, showing that the imaginary ring flip acts both as a chemical potential for visons and as an effective three-photon vertex akin to strong light-matter coupling. The novel coupling causes a structured high-energy continuum to emerge above the photon dispersion, which is naturally interpreted as three photon up-conversion in a nonlinear optical crystal.


I. INTRODUCTION
Quantum spin ice (QSI) is one of the best known three dimensional models realising quantum spin liquid phases.In addition to the gapped spinon excitations of classical spin ice, QSI is characterized by the presence of an emergent U (1) gauge field bearing a striking resemblance to compact scalar QED [1], complete with emergent analogues of Dirac monopoles (visons) and gapless gauge photons.Here we adopt the 'physical' naming convention in which the spinons are said to carry magnetic charge in the emergent QED, while visons carry electric charge.This convention is adopted to reflect the well-established (real) magnetic charge of the spinons in spin ice materials [2,3], and the more recent suggestion that visons may experience Lorentz force in the presence of an externally applied magnetic field [4,5].Correspondingly, we shall use the language of Refs. 2 and 6, in which the S z operators correspond to the magnetic field, while canonically conjugate rotor variables a, [S z , a] = −i, are understood as the U (1) vector potential of an emergent electric field [7].
In this work, we study a realistic model of breathing quantum spin ice, showing that within the spin liquid phase the breathing anisotropy and generic nearestneighbor couplings can conspire to create an all-in, allout (AIAO) bias in the electric field.If tuned beyond a critical strength in parameter space, we demonstrate that such a bias field induces a phase transition and the * als217@cam.ac.uk emergence of a vison crystal, or equivalently, a U (1) π/2 phase.
The AIAO bias echoes an existing model in the magnetic sector, which we summarize here as a useful analogy for our generalization to the electric sector.In certain rare-earth pyrochlore iridates (RE 2 Ir 2 O 7 , RE ∈ {Ho, Dy, ...}), classical spin ice is interpenetrated by a pyrochlore lattice of AIAO-ordered iridium atoms that generate a staggered, AIAO magnetic field h aligned with the local z axis on every pyrochlore site j [18][19][20][21][22][23], where η t is +1 on the A tetrahedral sublattice and -1 on the B sublattice.Both inversion and time reversal symmetry are explicitly broken, prejudicing 'in' (−) and 'out' (+) spinons to lie on different sublattices.The staggered field h drives a strongly first order monopole condensation transition [24] from the divergence free 2-in, 2out (2I2O) phase to a 3-in, 1-out (3I1O) spin-fragmented phase, ultimately reaching an AIAO-ordered phase if h is further increased [20,[25][26][27].
We show that an analogue of the AIAO magnetic bias field is generally present in the electric sector of a breathing spin- 1  2 quantum spin ice.The role of the iridium AIAO magnetic field is played, subtly, by the J z± S z i S ± j term (acting within conserved spinon subspaces), assigning different phases to spinon hopping processes S ± j depending on the state S z i of the spin at site i.When these terms are used to transport a virtual spinon around a closed loop in the standard fashion of a ring flip process, a nonvanishing phase correction emerges in the effective Hamiltonian.The main consequences of this novel contribution to the Hamiltonian may be summarized as follows: 1. Nearest-neighbor microscopic spin interactions on the Kramers breathing pyrochlore generically generate an AIAO electric bias field.
2. Numerics reveal that this electric bias field selects between four possible phases: the known U (1) 0 and U (1) π phases, and two unconventional, fragmented vison crystal states which we refer to as the U (1) ±π/2 phases.
3. These phases are separated by lines of liquid-gaslike first order phase transitions, terminating at a critical end point and becoming all continuously connected via the disordered high-temperature phase.
4. The bias field is decoupled from the emergent photon dynamics to quadratic order.
5. The leading order correction to the photon propagator comes at cubic order in the emergent gauge field a, representing a three-photon vertex correction.
6. Visons interact via an emergent Coulomb law, with the bias field acting (within a large-S description) only as a chemical potential for visons, prejudicing them to lay on a certain sublattice set by the original breathing anisotropy.
The paper is structured as follows.We introduce the model in Sec.II.The microscopic perturbation theory responsible for the AIAO offset is outlined in Sec.III.Sec.IV presents certain symmetries of the minimal model, which we employ to explore the parameter space more efficiently.We then study in Sec.V the phase diagram of this AIAO-biased quantum spin ice using semiclassical simulations [6,28].We find a sharp, liquid-gas like phase transition between the known U (1) 0 QSI phase and an unconventional π/2-flux phase.This transition appears to terminate at a critical endpoint, which we investigate using finite-size scaling.Sec.VI presents static and dynamic structure factors within the U (1) π/2 phase.In Sec.VII we obtain a functional form fo the vison-vison interaction potential in the new regime.Finally, we draw our conclusions in Sec.VIII.

II. BREATHING PYROCHLORE QUANTUM SPIN ICE
We begin by choosing a notation with which to index the sites of the breathing pyrochlore lattice.The spins of the model, residing on the pyrochlore sites, can be conveniently thought of as living on the links of a parent diamond lattice, as illustrated in Fig. 1 [29].The pyrochlore tetrahedra are centered on the diamond sites, and the smallest closed loops of links form (buckled) hexagonal plaquettes.Four such plaquettes can be chosen to enclose three-dimensional volumes, which we will refer to as voids [30].It will be most convenient to use the language of pyrochlore sites and tetrahedra for the perturbative derivation of the effective Hamiltonian in Sec.III, while Gray tetrahedra share corners which form the pyrochlore sites i (red sphere, right), which are in one-to-one correspondence with the links (red line, left) of the diamond lattice.Each diamond plaquette p (green even-sides hexagon, left) is similarly in unique correspondence with the hexagonal plaquettes of the pyrochlore lattice (green alternating-sided hexagon, right).Four pairwise adjacent plaquettes enclose a void v (purple region, bottom left).
the effective theory of the subsequent sections is more elegantly described in the diamond lattice language.
We obtain the Hamiltonian of breathing quantum spin ice by modifying the standard interaction terms in nonbreathing, Kramers pyrochlore spin ice (J zz , J ± , J z± , J ±± ) to take different strengths on the two inequivalent tetrahedra, which we name A and B [8,[31][32][33][34][35][36][37]: where, if µ(i) = 1, 2, 3, 4 denotes the pyrochlore sublattice of link i, then σ is a dummy index for the diamond (i.e., tetrahedral) sublattices A, B and ⟨ij⟩ σ denotes a pyrochlore bond belonging to a tetrahedron on sublattice σ.The sum over i in the fourth line runs over all pyrochlore sites, with the subsequent sum ⟨j → i⟩ σ running over the six (three A, three B) nearest neighbors of i.Note that each spin component is understood with respect to the local (pyrochlore sublattice dependent) axes tabulated in Appendix A. Our attention is restricted to the case of Kramers doublets; the additional conditions satisfied by non-Kramers and dipolar-octupolar doublets are given in Table III for completeness.We suppose that both J zz,A and J zz,B are much stronger than the remaining couplings.Without loss of generality let the A-sublattice coupling be stronger, J zz,A ≳ J zz,B > 0. If we neglect all other terms, the Hamiltonian becomes and the ground states belong to the canonical 2I2O ice ensemble.The energy of a spinon excitation becomes sublattice dependent, rendering the B-sublattice spinons (irrespective of their charge) less energetically expensive.The spinon gap is therefore controlled by J B zz .In this respect, we note that the known breathing spin- 1  2 pyrochlore material to date, Ba 3 Yb 2 Zn 5 O 11 , suffers from an exceedingly small value of J zz,B , leaving the A tetrahedra paramagnetically decoupled down to ∼ 0.1 K [38].

III. PERTURBATION THEORY
We next revisit the perturbation theory of Hermele, Fisher and Balents [1], starting from the general symmetry-allowed nearest-neighbor Hamiltonian in Eq. ( 1) compatible with the effective spin-1/2 C 3v symmetry.
In the standard fashion, we consider We construct an effective Hamiltonian that operates solely within the 2I2O spin ice ensemble using the standard expansion where P projects onto the 2I2O ice manifold, and the interaction term is We also remind in passing that second order perturbation theory is known to only generate farther range Ising (namely, S z ) interactions that -when sufficiently small -do not disrupt the QSL phase outright, but merely renormalize the emergent speed of light [39,40].In what follows, we suppose that the system remains continuously connected to the QSL phase, with the Ising interactions simply weighting the state distribution.We similarly neglect any four-S z and longer range Ising terms (e.g., those dubbed γ and I 2 in Table I) that occur at higher orders.
Selection of processes (i.e., ring flip loops and retracing paths) appearing up to fourth order (see also Table I), together with their associated (unordered) operator sets.Black and white symbols denote spins that are flipped in the process, respectively up-to-down and down-to-up.Circles denote pyrochlore sites, i.e. spins.Spins that are flipped twice (returning to their original state) are half-filled in black and white.Gray spins are not flipped by the process.Curved purple lines indicate the six nearest neighbors of a ladder operator S ± , over which the sums over i,j run.Loops a) and b) illustrate the two classes of terms arising from H 2 ± H 2 z± that are responsible for generating complex ring flips, as explained in the main text.Loop c) is a fourth-order 'teardrop' correction to the conventional ring flip term.Self-retracing loop d) generates only an Ising (namely, S z ) interaction.Note that the smaller (i.e., A) tetrahedra are taken to face out of the plane of the page, as in right panel of Fig. 1.
The processes generated at third and fourth order have sufficiently many vertices to create a virtual spinon pair, transport the spinons around a non-trivial (i.e., nonself-retracing) loop on the diamond lattice, and annihilate them.The prototypical example of such a process is the well-known ring flip term from H 3 ± , which generates g (3) and the corresponding operator (where 1, 2, . . ., 6 label the six sites around the hexagon ).Another example at higher order comes from H 4 ± , which generates 'teardrop' ring flips (see Fig. 2d), as well as some additional trivial energy shifts.
An illustrative selection of the resulting terms is given in Fig. 2, and detailed in Table I.The coupling H ±± can be shown (see Appendix B) to generate only multi-S z interactions up to fifth order, and we therefore neglect it for the purpose of our discussion.Importantly, in addition to the canonical inversioneven g terms, we obtain an imaginary ring flip ig ′ O − O † from the eight-spin ring exchange H 2 ± H 2 z± .This term arises from the sum over all strings of the form , where 1-6 refer to the pyrochlore site labeling of Figs. 2 and 3. When both i and j lie on the ring (i.e.i, j ∈ {1, 2, ..., 6}), as they do in Fig. 3a,b, repeated applications of the identity S z S ± ≡ ±S ± /2 replace both factors of S z by scalars.The effective operator is then proportional to O , with a possibly complex coefficient.Such complex ring flips are generally forbidden by inversion symmetry, but are allowed on the breathing pyrochlore where this symmetry is explicitly broken.We illustrate here an example of a non-reciprocal process responsible for these imaginary ring flips.Consider the process shown in Fig. 3a In the term (5), the rightmost factor of H 0 measures two spinons, one on each of the A and B tetrahedral sublattices, giving a factor of 2/(J A zz + J B zz ).The corresponding propagator in (6) instead measures two spinons on the B sublattice, giving a factor of 1/J B zz .The terms therefore do not cancel.
There is an important subtlety associated with orienting the ring flip operators O on the breathing pyrochlore lattice.The loops in Fig. 2 are oriented with the smaller A tetrahedra facing out of the plane of the page, with ladder operators arranged such that a counterclockwise traversal of the ring raises (lowers) the spins encountered when moving from an A (B) tetrahedron to a B (A) tetrahedron.The opposite convention would generate instead O † .For consistency, it follows that g must be even under exchange of tetrahedron sublattices, while g ′ must be odd.
For completeness, we mention that these H 2 ± H 2 z± -type ring processes also generate ring flips of the form S z S z O where the dangling S z factors are not part of the ring (i, j / ∈ {1, 2, . . ., 6} in the notation above, see Fig. 3c).We interpret these as ring flips modulated by S z S z correlations in the local environment, which we elaborate on in Appendix B. We will ultimately drop such eightspin operators, restricting our attention to the 'pure' ring flips [6,41,42].Note that for any C ∈ C \ {0}, terms of the form CS z O+C * S z O † violate time reversal symmetry and do not appear in the effective Hamiltonian; it follows that 1-on, 1-off complex loops, as the case illustrated in Fig. 3d, give no overall contribution.
If we further parameterize the couplings by introducing the anisotropy parameter κ, where (1 + κ)/(1 − κ) = J zz,A /J zz,B = J z±,A /J z±,B = J ±,A /J ±,B , it is possible to give an asymptotic expression for g, g ′ in the vicinity of κ = 0. Denoting the average couplings with an overline (e.g., J zz = (J zz,A + J zz,B )/2), we have As required by the definition of O, it is true to all orders that g ′ is an odd function of κ, whereas g is even.
Take note that the fourth-order contribution to g has a far larger combinatorial prefactor than the third-order term (216 as opposed to 12).This is suggestive of a slowly converging perturbative series in cases where J ± /J zz is moderately large.Indeed, recent gauge mean field theory calculations [43] suggest that the U (1) π phase remains stable up to J ± ≈ J zz .
IV. THE g + ig ′ MODEL In the rest of the discussion, we restrict our attention to a minimal model containing g and g ′ as the sole parameters: The sign of g ′ reflects a choice of labeling for the A and B tetrahedra.The earlier convention J ±,A ≥ J ±,B renders, without loss of generality, g ′ non-negative.It will be later convenient to express g + ig ′ = ρ g e iϕg (recall that g and g ′ are real).

A. Large-S limit and semiclassical theory
We apply the canonical Villain expansion [44], introducing an operator a canonically conjugate to S z (with [a, S z ] = i).As the a variables are quantum rotors, they are compact; our work, following standard convention, understands a to be an operator on the Hilbert space L 2 (U (1)) of square integrable functions on the periodic [0, 2π) interval, which will ultimately be represented as a c-number in U (1) = [0, 2π) in the semiclassical theory.Ladder operators for spin S are then replaced by the expression where s = S + 1/2 (see also Refs. 1,6,and 41).We expand O p to leading order in 1/s to obtain a modified form of the canonical U (1) lattice gauge theory, where we introduced the lattice curl ∇ × p a = j∈p (−) j a j , with an identical sign convention to that which was used in the definition of O, chosen to respect the orientation of the plaquette O p operators.In the second line above we reparameterized for convenience the coupling constants as g = ρ g cos(ϕ g ) and g ′ = ρ g sin(ϕ g ).With this definition, it becomes clear that a is a compact U (1) gauge field -∇ × a is the discrete line integral of a around a closed loop, i.e. a gauge flux.
Notice the subtle importance of working first on the perturbation theory of the spin-1/2 system, to produce an effective Hamiltonian that contains terms with the same number of spin operators, and then taking the large-S limit.While we have combined third and fourth order terms, which nominally have 6 and 8 spin operators respectively, we argued in Sec.III that the interesting and novel contributions from the latter arise when 2 of the 8 spin operators (namely, the S z terms) have a known expectation value and can therefore be included in the numerical coefficients.This reduces the 8 spin operators of interest to 6 spin operators, thence Eq. (9).Turning then to large-S theory produces a well-defined Hamiltonian where all the contributions are of order s6 .Of course, the fact that some terms come from third order perturbation theory, and some from fourth order, remains reflected in the parametric form of the Hamiltonian coefficients.
The Hamiltonian in Eq. ( 11) is cast in terms of a U (1) lattice gauge field a.It will be useful to the later discussion to introduce the unitary C, which rotates all spins about the local x axis.
This operation reverses all S z moments, inverting the sign of all magnetic monopoles.Further, its action on S ± , suggests that Ca = −a, reversing the electric field strength.Inspired by the observation that C reverses both electric and magnetic field lines, we name this operation charge conjugation.It follows immediately from this discussion that C maps (g, g ′ ) → (g, −g ′ ).Further, it may be straightforwardly verified that, in the effective Hamiltonian, C is equivalent to the formal inversion operation in which A and B sublattice labels are interchanged, in a QSI analogue of QED's charge-parity symmetry.
Note that H zz and H ± are invariant under C.It is therefore a symmetry of all pure-XXZ spin models, and indeed all dipolar-octupolar spin ice Hamiltonians.This highlights breathing Kramers QSI as a promising system for generating complex ring flips [45]: both inversion and charge conjugation symmetry need to be broken explicitly.

B. Definition of the periodic electric field
As a is a rotor variable belonging to U (1), ∇ × p a is ambiguous up to factors of 2π.Indeed, any summations of a must be understood as an addition of phases, intrinsically modulo 2π.This ambiguity makes it somewhat subtle to define the presence (or absence) of an electric monopole, i.e., a vison.
There is no ambiguity if one works directly with the Wilson loops exp(i∇ × p a) (deriving indeed from the product of spin operators around an elementary hexagonal plaquette).Thinking of a as a lattice gauge field, it may be readily observed that this quantity is the total phase acquired by the transport of a virtual U (1) charge around the plaquette p, the Wilson loop which contains all the gauge-invariant information of the lattice gauge theory.
For each void v of the lattice (see Fig. 1), we introduce the lattice divergence of ∇ × p a to be where η v = +1 on the A-sublattice voids and . Since all a variables are defined modulo 2π, we are forced to consider Ṽ modulo 1. Due to the geometric cancellation illustrated in the top panel of Fig. 4, we have Ṽ ≡ 0 (mod1).
Nonetheless, it is clearly useful to be able to measure the sense in which the divergence is zero: there is a physical difference, for example, between a void with all plaquette fluxes close to zero and a void with predominantly π/2 fluxes.In order to construct a vison charge operator that is sensitive to this difference, we have to first choose a chart for the lattice U (1)-valued field ∇ × a, i.e., a partial mapping to an open subset of the real axis; we shall refer to ∇ × a in a given chart as e p , to indicate this subtle but important change.The convention used One may regard the link variables a as either U (1) rotors or real numbers.Bottom panel: Summation of phases involved in V (v) (Eq.( 19)) as viewed in the (−π, π) chart.The topology of U (1) is depicted using a circle, with the four red arrows representing the four ep fluxes associated with the void v.As discussed in the main text, the path must return to its starting point to remain consistent with ∇ •v ∇ ×p a = 0(mod1).The removal of ±π phases forbids any red arrow to enter the midpoint of the circle -this is in effect a topological singularity about which a winding number can be defined.
in the majority of the literature [28] (implicitly or otherwise) is to take e p ∈ (−π, π), which for weak fields is indistinguishable from standard QED.Regarding e p ∈ R as ∇ × p a in a particular chart, the vison charge operator can then be defined as which in general takes any value in Z.This step makes V sensitive to the journey taken by the phase as it loops back to zero, identifying the vison charge as a winding number (bottom panel of Fig. 4) that, in the chart (−π, π), may only take values −1, 0, 1.
For consistency with the existing QSI literature, we will also work in the chart-fixed '0-flux convention' in which e ∈ (−π, π), leaving further comments on chart dependence to Appendix C. As a manifestly chartdependent quantity, V should be interpreted as a vison marker, rather than the exact vison charge operator itself.Where appropriate, we will also take steps (see Appendix C) to illustrate that any electric sector behavior we remark on is also observed in terms of purely chartindependent variables.

C. Quadratic Hamiltonian
In the special cases (g, g ′ ) = (±1, 0) and (g, g ′ ) = (0, ±1), the ground state of the semiclassical Hamiltonian in Eq. ( 11) can be read off as ∇ × p a = 0, π, ±π/2, ∀ p, which are the four uniform states consistent with the constraint ∇ • v ∇ × p a = 0 (mod 2π).We name these states U (1) 0 , U (1) π and U (1) ±π/2 , respectively.Using these four states as variational Ansätze, one finds that is an upper bound on the ground-state energy.These uniform-flux states are robust against perturbations.Indeed, if we write the Hamiltonian (11) in terms of e p (upon fixing a chart) and we consider small fluctu- it becomes clear that in the vison-free ∇ • v e p = 0 sector (as appropriate for the low energy U (1) 0 state), the g ′ term contributes only at higher than quadratic order.Note that the divergence term, is a multiple of the large-S vison charge operator V (v).The sum appearing in Eq. ( 18) runs over A voids only, and we can thus identify the linear-order term as a staggered electric potential that favors positive visons on the A voids (and therefore negative ones on the B voids, by global charge neutrality).Due to the vison quantization, there is no linear response to a small g ′ perturbation of the U (1) 0 phase, up until g ′ reaches a value comparable to the vison gap.Similarly, one can show that the U (1) ±π/2 phases have vanishing linear response to g perturbations.
To measure the excess of positive visons on the A sublattice, we define the order parameter as the average value of the vison charge operator restricted to voids lying on the A sublattice: which has vanishing expectation value at all temperatures when g ′ = 0, and tends to 1 at low temperatures in the vison crystal phase.This quantity essentially measures the electric AIAO polarization (recall that by global vison charge neutrality, V B = −V A ), playing the role of the staggered magnetization for a Néel antiferromagnet.

D. Enhanced spectral periodicity and duality in the minimal model
Examples of gauge field configurations for the U (1) lattice gauge theory on the diamond lattice, generating uniform π/2 and π flux phases are given in Fig. 5. Though the fluxes themselves are commensurate with the periodicity of the breathing lattice, the enlargement of the unit cell retains a physical significance in the periodicity of the spinon continuum [46,47].In essence, the Aharonov-Bohm phase acquired from spinon transport fractionalizes the lattice translational symmetry into a projective space group representation, which only 'squares out' to a true representation when an enlarged unit cell is considered.The zero-flux phase can be described by tiling a primitive unit cell with a rr ′ = 0 over the whole lattice, while the smallest link variable configuration of a π flux phase requires two primitive unit cells.It can be verified explicitly that the smallest possible unit cell describing a U (1) π/2 phase consists of four primitive unit cells, distinguishing the state from the U (1) 0 and U (1) π phases.
Fig. 5 may also be interpreted as a duality of the effective Hamiltonian.This is achieved by recognising that the Abelian group structure of the lattice gauge theory extends to the spin model.If a(j) is a lattice gauge field, then the unitary operation on # Charge conjugation, combined with the π/2 rotation duality, makes the physics of the g, g ′ phase diagram D 4 symmetric.The g, g ′ model can therefore be fully characterized by the half-quadrant g > 0, 0 < g ′ < g.In the special cases ϕ g ∈ π 4 Z 8 , these dualities are elevated to Z 2 symmetries of the Hamiltonian, which we further distinguish into the 'axis' symmetry lines g ′ = 0, g = 0 and the 'diagonal' symmetry lines g = ±g ′ .While the axis symmetries remain intact down to zero temperature, we will show that the diagonal symmetries are spontaneously broken below a critical temperature.

V. LARGE-S PHASE DIAGRAM
We are finally in the position to study the phase diagram of the model in the large-S limit.We use a semiclassical Monte Carlo approach developed by Szabó and Castelnovo [28], based on the large-S expansion [6], adapted to our effective quantum Hamiltonian in Eq. ( 9).We perform classical Monte Carlo on the low-energy effective Hamiltonian (Eq.9), capturing the gauge-sector dynamics of the quantum theory in a large-S sense [49].The ring flip operator O p = S + S − S + S − S + S − is regarded as a complex number of norm ≤ 1.The presence of a vison in void v is measured by V (v) = p∈v arg O p /2π, using the chart arg : C → [−π, π).We initialize the system in a random state, then lower the temperature from T = 10ρ g to T = 10 −5 ρ g in logarithmically spaced steps.At every temperature, we run a number of annealing steps known from benchmarks to equilibrate the system (typically ranging from 128 at high temperature to 512 at low temperature).At the end of a full cooling cycle, we then re-heat the system up to the initial temperature, following the same temperature steps and comparing order parameters throughout.The absence of any hysteresis is used as indication that equilibrium has been likely attained in the simulations.
A sharp phase transition is observed in Fig. 6 at g = g ′ at low temperatures, separating the conventional U (1) 0 phase from a phase in which all A sites are occupied by visons.In the latter, all electric fluxes take the value π/2, hence the name U (1) π/2 phase, previously identified using a projective symmetry group approach in Ref. 47.
The transition from 0-flux to π/2-flux phases may be understood as the chemical forcing g ′ overwhelming the intrinsic vison chemical potential µ g ′ =0 = 7.872367608(68)g [28].Neglecting any cohesive energy from interactions, the phase transition would occur at g ′ /g = µ g ′ =0 /g ≃ 7.872 (ϕ g ≃ 0.460π).However, we know from the (unitary) duality operation R a π/2 C, defined in Eqs. ( 21) and (13), that the phase-boundary structure of the global phase diagram must be symmetric under g ↔ g ′ .It follows that the critical g ′ /g is at most 1 (since g/g ′ is also a critical point).Our numerics clearly show g ′ /g = 1 to be the critical point, emphasizing the substantial renormalization of short-ranged µ by interactions.
We verified that the zero-temperature system energy saturates the trivial bound in Eq. ( 16) from Sec. IV C to within Monte Carlo error, suggesting that we have successfully identified all relevant low temperature phases.Further, it suggests that the smaller of g, g ′ plays a trivial role at zero temperature, merely shifting the ordinary QSI vacuum.
Fixing the overall energy scale ρ g , we identify a characteristic temperature T pm associated with the crossover to the disordered phase.Guided by the quadratic theory of Sec.IV C, we are able to establish a phenomenological model that captures this behaviour.Recall that expanding the Hamiltonian about the 0-flux vacuum, Eq. ( 18), yields g as the tension of the electric field (i.e., the permittivity of free space) and g ′ as a staggered chemical potential for visons.Applying the duality R a π/2 interchanges g and g ′ , rendering g ′ the effective field tension for excitations above the π/2-flux state, which is indeed the preferred ground-state flux configuration when g ′ > |g|.It may be seen by a straightforward generalization that the quantities T = max(|g|, |g ′ |) and U s = min(|g|, |g ′ |) act respectively as field tension and chemical potential for all values of ϕ g .We see that T pm ∼ √ 2T c (1.3T − 0.3U s ) captures well the general trend of the crossover.Though the dominant contribution to the crossover's shape is the weakening of field tension, including a small U s contribution appears to improve the visual agreement.
Raising T above T pm brings about a vison liquid state, featuring a random disordered electric field.We identify this with spinon-free classical spin ice.The transition to the true paramagnetic phase [50] is not observable in our approximation scheme, as we work within the no-spinon ensemble of states.

A. Return to microscopic coupling constants
It is useful to now reconsider the phase diagram in terms of the original microscopic variables of Sec.II, J zz , J ± , and J z± .A particular slice of parameter space is shown in Fig. 7, using the uniform scaling parametrization of Sec.III.Close to the critical point κ = 1, Jz± / Jzz − ∼ 0.2, it is possible to use J z± /J zz and κ−1 κ+1 as proxies for g, g ′ , see Eqs. (7) and (8).We thereby obtain a clear physical interpretation of the D 4 duality structure derived in Sec.IV D -the four white lines (two dotted, two solid) are lines of emergent symmetry generated by the conjugate-shift operations R n a π/2 C, n ∈ Z 4 .Admittedly, the value of J z± needed for g ′ to overcome g is very large compared with J ⊥ .There is therefore a risk of the Ising interactions I 2 causing the system to order before any appreciable g ′ can be generated.We stress, however, that the calculation here is truncated at fourth order.We expect very substantial renormalization of g and g ′ at higher orders.

B. Critical end point
The duality arguments of Sec.IV D are consistent with a continuous connection between the U (1) 0 and U (1) π/2 phases via the phase-incoherent spin ice paramagnetic phase -there is no symmetry breaking across the transition.The line of first order phase transitions at g ′ /g = 1 therefore ends at a critical end point.
We probed the critical exponents of the end point using finite size scaling analysis.We used 3D periodic systems of cubic side lengths L = 4−12 unit cells.At each system size, we initialized 256 independent systems to random initial states, then performed a temperature sweep from 10ρ g down to 0.3ρ g , then back up to 0.45ρ g .Annealing times were chosen to remove any hysteresis in the observables.In the case of specific heat, we computed an ad hoc cost function [51] designed to estimate the quality of the scaling collapse over a grid of T c , α and ν values.From this tableaux, we estimated both the optimal fitting parameters and their uncertainties, including T c , based on the width of the optimal fit basin.We repeated this analysis for the staggered polarization data ⟨V A ⟩, with the additional fitting parameter ⟨V A ⟩ c , obtaining also values for β and γ.
Our results (shown in Fig. 8) are suggestive of 3D Ising criticality, finding its critical exponents β, γ and ν to lie within our error bars.Recall from Sec. IV D that on the line ϕ g = π/4, the operator R a π/2 C is precisely a Z 2 symmetry that is spontaneously broken across the g ′ = g line -consistently with Ising criticality at the end point.The story is, however, not quite so neat: the exponent α This peculiarity is perhaps less surprising against the background of Debye-Hückel plasma criticality, which has been variously reported as Ising, Gaussian, crossover-like and tricritical [54][55][56], depending sensitively on the form of the short-range cutoff [57].between the U (1) 0 and U (1) π phases above T c .However, the critical value of the vison order parameter ⟨V A ⟩ c ≃ 0.38, and perhaps even more surprisingly the high-temperature phase-incoherent state sees an even smaller value, ⟨V A ⟩ ≃ 0.
To understand this, let us firstly note that systems with ϕ g = x and ϕ g = π/2 − x are dual to one another.Indeed, the system energy seen in Fig. 6b is symmetric about the coexistence line ϕ g = π/4 to within simulation accuracy, and the expectation values ⟨cos e⟩ and ⟨sin e⟩ in Fig. 6c,d are related to one another by mirroring about the coexistence line.
Numerically, the order parameter V A is computed by V = p arg(O p ), where O p is a complex number formed by a product of S x ± iS y terms using the chart arg(z) ∈ (−π, π].As already discussed in Sec.IV B, this is a choice of chart convention that privileges zero, as at high tem- perature, the system samples a uniform distribution on (−π, π].This average contains no physical information, only a marker of which chart was used for e p .As the temperature is raised, ⟨V A ⟩ 'crosses over' from a useful indicator of the number of topological defects to a completely unphysical constant.This emphasizes the ill-definedness of visons in the high density regime [58].For further discussion of the precise meaning of the vison in QSI lattice gauge theory, see Appendix C.

VI. ELECTRODYNAMICS AT FINITE g ′
The dynamical structure factor of our system (see Fig. 9) continues to feature the known large-S photon dispersion in the presence of g ′ .We interpret these plots as probes of the photon spectrum alone.The visons possess no kinetic properties in the large-S limit and only hop via instanton processes; we are therefore unable to meaningfully probe their real-time evolution [28].For this reason, our simulations show a gapless photon dis-persion even when temperatures are high enough to have a substantial vison population [28], as opposed to the Debye-screened mass one would expect from a charged plasma [41,42].
The ⟨S z S z ⟩ correlations are invariant under the R a π/2 duality transformation, and therefore the behavior in the U (1) π/2 phase is identical to that shown in Fig. 9 (with the roles of g and g ′ reversed).This statement was verified numerically as a benchmark.
At g ′ = 0 we see that the familiar large-S photon dispersion is in excellent agreement with the analytic result at low temperature.Less obviously, the vast majority of the spectral weight remains identical to the g ′ = 0 large-S photon curve even when g ′ /g is of order 1.It is only at higher temperatures, where the photon states become highly populated, that the e 3 p vertex -see Eq. ( 23) -begins to have a measurable influence.In this regime we observe a continuum above the photon line that terminates abruptly at twice the renormalized maximum dispersion energy.There is some structure discernible in the continuum, which reflects the density of two-photon states.Excitingly, these anomalous features are visible in the experimentally accessible S z correlation function, and could in principle be used to detect a system sufficiently close to the phase transition between the U (1) 0 and U (1) π/2 phases, at intermediate temperatures.
The leading order contribution due to the g ′ term in the U (1) 0 phase is cubic, which we suggest is analogous to strong light-matter coupling in a nonlinear crystal.From this perspective, the high-energy features are naturally interpreted as a two-photon continuum arising from three photon up-conversion [59], which we compute analytically in the large-S limit for comparison hereafter.

A. Large-S calculation of two-photon continuum
We show here how the continuum seen in the ⟨S z S z ⟩ dynamical structure factor in Fig. 9 can be explained within perturbative large-S field theory.We perform the calculation for the U (1) 0 phase, noting that the calculation would be identical in the U (1) π/2 or U (1) π phases, up to duality relations.
We decompose the large-S Hamiltonian into quadratic and interacting parts in the standard fashion: The key insight is that the vison contribution V (v) is a non-dynamical, topological term that does not enter the perturbation theory.
To pick out the most relevant terms from this expansion, we follow Ref. 6 and perform naive RG scaling on where b > 1.In particular, we note that a pair of contracted three-photon vertices scales as b −1−d s−1 , which is the same as a Hartree-Fock vertex.It follows that we need only consider the (momentum-space) diagrams a), b) and c) shown in Fig. 10: FIG. 10.Leading-order loop corrections in the large-S field theory.
The Hartree-Fock contribution a) is known [6] to be a renormalization of the photon continuum.The contribution b) can be shown to vanish exactly because the central propagator carries zero momentum.Remarkably, only the leading-order non-trivial term c) is found to give a contribution from g ′ to the spectral weight (see Appendix G for details).Looking at Fig. 11, we see excellent agreement between analytical and numerical results.The leading-order calculation somewhat overestimates the maximum energy of the continuum, and in particular the spectral weight at high energy, we expect however that such discrepancies will be resolved by dressing the propagator at higher orders in perturbation theory.

B. Static Structure Factor
Perhaps the most direct measurement of vison crystal formation is the appearance of Bragg peaks in the ⟨ee⟩ correlations, shown in Fig. 12. Starting from ϕ g = 0, screening of the Coulomb interaction by the Debye plasma brings about a Lorentzian blurring of the ⟨ee⟩ pinch points.This crossover is controlled by two factors: the lowering of the vison gap upon increasing g ′ , and the simultaneous increase of the ratio T /g (recall that we are keeping the ratio T /ρ g fixed at 1/2).Both effects independently lead to an increase in the vison density.Interestingly, we observe a blurring of the pinch points even if g ′ is increased at constant T /g (likely due to charge disorder upon approaching coexistence at the transition line).
On either side of the coexistence line ϕ g = π/4, a Bragg peak emerges at the Γ point due to the statistically significant excess of net vison charge on the A sites (as already mentioned in preceding sections, e(k = 0) is indeed proportional to the vison order parameter V A ).
At the coexistence point ϕ g = π/4, we observe rod correlations associated with domain walls that are frozen in as they become metastable at low temperature.Rods extend in the [110] and (to a lesser extent) in the [111] directions, reflecting an anisotropic energy cost to domain wall orientation.

VII. VISON INTERACTIONS
Finally, we study the energetics of the visons in our system -namely their chemical potential and interactions -within the large-S approximation.In order to do so, we create spin configurations where we arrange visons in a polarization-free octupole pattern and measure the total field energy at vanishing temperature, which we then decompose into a chemical potential part and an interaction part by changing system size.
We place positive and negative visons respectively on the void sublattice where their charge is preferred by the staggered chemical potential g ′ , making them metastable.It can be shown that a zinc blende superlattice (see Fig. 13), corresponding to the original diamond lattice rescaled by a factor of 4m + 1, m ∈ Z, places the rescaled A(B) sites coincident with the A(B) sites of the original lattice.We exploit this to generate a sequence of increasingly dilute artificial vison supercrystals, from which we may analyse the scaling of the crystal cohesive energy (thus being able to differentiate between chemical potential and interaction effects).
We do not have access to a closed form for the vison creation operator, although a good approximation may in principle be obtained from the Green's function of the lattice Laplacian [1].Here we adopt instead a different procedure to create the desired artificial vison arrangement: VC0 We find a polarization-free arrangement of electric Dirac strings that introduces the desired vison structure.
VC1 For every pair of connected visons, we simultaneously rotate the fluxes of all plaquettes crossed by the Dirac string until branch cuts are crossed at the beginning and end of the string.This places the system in an excited state corresponding to a superposition of the desired vison (super)crystal structure and some distribution of photons.
VC2 We anneal the system from an initial temperature T hot > ρ g down to T cold ≪ ρ g , allowing the energy of the metastable configuration to be read off once the photons have thermalized.

A. Numerical results
The scaling of the system energy with the linear dimension L of the system is show in Figs. 14 and 15.
Inspired by the quadratic theory in Eq. ( 18), we model the visons as a Coulombic plasma with cohesive energy where Z v,w = ±1 is the sign of the charge of the visons labelled by v, w, r v,w are their positions and N vison is the number of visons [60].It follows immediately from trivial rescaling of the Hamiltonian that, at zero temperature, q 2 e (g, 0) ∝ g and µ(g, 0) ∝ g.The weak-field expansion in Eq. ( 18) suggests that q 2 e is independent of g ′ to leading order, with corrections arising from the polarized region in the immediate vicinity of a vison.
The numerical scaling of the vison energy with system size reveals that this intuition is essentially correct.The Coulomb functional form is an excellent fit to the system energy in Fig. 14a at all length scales, though the short-range energy for L = 5 is marginally overestimated.Large values of the electric field are more heavily penalized by the quadratic theory than they are by the original cosine potential; near-neighbor vison dipoles are therefore more favorable in our numerics than the Coulomb estimate would predict.Indeed, we believe that this excess of vison dipoles is responsible for the homogeneous background seen in the static e correlators of Fig. 12.
The dependence of the zero-temperature vison chemical potential on g ′ , shown in Fig. 14c, is very robustly linear, suggesting that the chemical potentials associated with g and g ′ are essentially decoupled.The gradient is slightly less than the quadratic estimate (−1), reflecting the influence of nonlinearities in sin(e).0.000π 0.500π Fig. 14b reveals the best-fit vison charge q 2 e to be very close to the quadratic estimate at all values of g ′ .The apparent dependence of q 2 e on g ′ is due to the failure of the simple Coulomb model to properly capture the higher-order moments of the vison-vison interaction.We see that as short-ranged points are excluded from the fit, the dependence disappears.Fig. 15 revisits this issue, subtracting off the robustly known chemical potential component to isolate the interaction part of the vison energy.We clearly see L −2 scaling consistent with the coupling of local dipole moments to the inhomogeneous field in the immediate vicinity of a charge.Though the fit data points at larger radius appear to cross over to Coulomb 1/L scaling, this is an artefact of the Gaussian error model used for the linear regression.Perhaps more significantly, the L −2 scaling remains consistent with all fits to within 1σ at all lengths.We can neither confirm nor rule out some Coulomb contribution arising from g ′ (we leave such costly numerical exploration for future work).
Combining this measurement of the vison charge with the photon dispersion measurement from Sec. VI, we are now in principle able to compute a large-S approximation of the fine-structure constant.There is, however, a subtlety arising from the electric-magnetic duality: In the presence of both gapped magnetic and electric monopoles, there are a priori two dimensionless numbers of interest, the 'electric fine structure constant' α e = q 2 e /ℏc and the 'magnetic fine structure constant' α m = q 2 m /ℏc, where q m is the charge of a spinon.Dirac's quantization condition imposes [61] a reciprocity relationship on the two fine structure constants, This condition forces at least one of α m and α e to exceed the confining threshold α c ≃ 0.21 [61][62][63][64], seemingly suggesting that either spinons, visons or both are confined.Some care should be taken before jumping to this conclusion, however.Though it is known that matter-free quantum U (1) lattice gauge theory has a strong-weak duality mapping electric and magnetic fields to one another [1,[65][66][67], it is a decidedly different question to ask whether this duality survives in the presence of dynamical matter.Previous numerics have shown coexisting, Coulombic spinons and visons even with α e ≫ 1 [61].
All this is to say that we should not naively interpret α e as a tuning parameter for vison confinement.Substitution of the large-S speed of light ℏc S=∞ ≃ 0.15ga 0 [6,41] yields α m ≃ 21.This number is independent of g and g ′ , and well above α c .Nonetheless, the Coulomb form of the vison potential is characteristic of deconfined visons.
The Dirac quantization condition Eq. ( 25) allows us to calculate the value of α m = 1/(4α e ).At S = ∞, it is of order ≃ 0.01, interestingly close to the value obtained from recent exact diagonalization (α m,ED = 0.07 [39]).Accounting for Hartree-Fock corrections to the speed of light [6] renormalizes α e to ≃ 7.7; in surprisingly reasonable agreement with 1/(4α m,ED ) ≃ 3.6.e /(La0) + µ to the data points with L ≥ 17.This model marginally overestimates the cohesive at L = 5 (see the data points outlined by a solid black line box).Note that qe may be interpreted as the elementary vison charge in Gaussian CGS units.The inset shows the diamond superlattice arrangement of visons, see also Fig. 13 for details.The lower panels show the dependence of b) the Coulomb constant q 2 e and c) the chemical potential µ on g ′ .Panel b) shows the dependence of the q 2 e fits on the short-range interaction cutoff, suggesting convergence to the quadratic estimate q 2 e = πg (horizontal dashed black line).Panel c) shows the dependence of the vison chemical potential µ on g ′ , where a linear best fit µ(g, g ′ )/g = −1.1154(2)g′ /g + 7.873(1) appears to be in excellent agreement with the simulations.These results are independent of the cutoff choice.For all plots, the absence of error bars indicates uncertainties smaller than the marker size (for reference, the standard error on the µ fits is of the order of 10 −4 g).

VIII. CONCLUSIONS
In this work, we studied the behaviour of (Kramers) quantum spin ice on a breathing pyrochlore lattice, by means of perturbative arguments and large-S approximations.We demonstrated that the broken inversion symmetry allows for a new term (ig ′ (O − O † )) that would Uv − ∂µ ∂g g − ∂µ ∂g ′ g ′ with g ′ .b) Log-log plot of the gradient of the interaction part of the vison energy with respect to g ′ against linear system size L, exhibiting a potential crossover from L −2 to L −1 scaling.The error bars reflect 1σ standard error from linear regression.
have otherwise been forbidden.As shown in Fig. 7, the anisotropy itself needs not be particularly large for g ′ to be comparable with g.We therefore suggest that the U (1) π/2 phase could be most easily realized in a 'minimally breathing' pyrochlore compound -the effect requires breaking inversion symmetry without decoupling the A tetrahedra.
The g ′ interaction effectively frustrates the emergent electric field, as there is no path in state space that can smoothly interpolate between the zero-flux state and the π/2-flux state.Instead, the fluxes jump abruptly between the four values allowed by the vison quantization condition.At finite temperature, the phases are separated by first-order liquid-gas like phase transition lines, terminating at critical end points which we tentatively placed in the 3D Ising universality class based on a finitesize scaling analysis.
When g > |g ′ |, we found that the interactions between visons are well approximated by an essentially g ′independent Coulomb law, in addition to a small correction term arising from deviations from the Coulomb law at short range.The more significant effect of g ′ is the emergence of a staggered chemical potential for the visons.It follows from duality that the excitations above the U (1) π/2 vacuum are similar -one may readily view them as vison holes interacting with each other via a Coulomb potential with the roles of g and g ′ interchanged.
The behaviour we have uncovered can be seen as the (quantum) electric analogue of spin fragmentation in classical spin ice [20,26,68,69], occurring for instance in iridate pyrochlores, where the interpenetrating AIAO magnetic bias field causes a monopole crystal to condense while the spins retain a classical U (1) fluctuating component.Similarly in our system, the U (1) π/2 vison crystal phase remains a spin liquid dual to the original U (1) 0 phase, undergoing fragmented Coulomb phase dynamics.Unlike the magnetic sector analogue, the electric phenomenology is not due to an extrinsic (Ir magnetism) contribution but it is inherent to the quantum spin ice system and its interactions; moreover, in our case there is a microscopic one-to-one mapping between the effective Hamiltonians of both the U (1) 0 and U (1) π/2 phases.
We further studied the dynamical structure factor of the system, revealing, to leading order, a photon dispersion independent of g ′ .We interpret this as the longwavelength photons decoupling from a deep-UV fluctuating background, just as real photons do in ionic crystals.More interestingly, we found a correction at order e 3 that directly influences the photon self-energy, modifying the experimentally accessible two-point S z correlator.This correction is in essence a QSI analogue of second harmonic generation in a nonlinear crystal.In the presence of a g ′ e 3 vertex, two photons may be up-converted to a single, higher-energy photon lying well above the original photon dispersion.
Realization of non-negligible g ′ requires an admittedly fine-tuned energy hierarchy: J zz,B /J zz,A must not be so small as to decouple the tetrahedra, yet also not so close to 1 that inversion symmetry suppresses g ′ .Nonetheless, the apparent ability of applied external electric fields to change the sign of the effective value of g [70,71] suggests that in the presence of non-zero J z± , an electric field tuning g between U (1) 0 and U (1) π could drive the system into an intermediate U (1) π/2 phase.
In our work, we have identified a number of distinctive phenomena that could be used as possible experimental probes to test the effects of g ′ in candidate materials.To maximize the chances of hitting a 'sweet spot' in parameter space, we suggest that a candidate material should be subjected to an applied [100] electric field, which is known to effectively tune the g parameter [70].We propose that an experiment sweeping g from negative to positive should search for three phenomena: a) two divergences in thermodynamic observables as the system transits across the phases U (1) 0 → U (1) π/2 → U (1) π ; b) the appearance of a structured high-energy continuum in ⟨S z S z ⟩ correlations (see Figs. 9 and 11), which can be accessed, e.g., using inelastic neutron scattering; and c) rod-like correlations from domain wall formation in ⟨ee⟩ at criticality (see Fig. 12).Though these 12-spin correlators are for all practical purposes inaccessible, we suspect that some vestige thereof may yet appear in the ⟨S x S x ⟩ and ⟨S y S y ⟩ channels.
This work highlights the richness of emergent QED in frustrated magnets with broken inversion symmetry, representing an unconventional lever with which the strongcoupling QED of QSI may be tuned towards yet new phases and phenomena.
It is conventional to write the QSI Hamiltonian with respect to the local reference frame for the spins.The convention adopted here is the same as in Ref. 28: S j = α=x,y,z S α j e µ(j),α , where µ(j) denotes the pyrochlore sublattice of site j and Further, we introduce the four inequivalent vectors connecting the centre of a tetrahedron to the four nearest-neighbor pyrochlore sites at its corners.This section presents the finer details involved in obtaining the quoted values for g and g ′ in perturbation theory.Many of the more verbose calculations involved were performed using our noncommutative algebra package commutation [74].
For completeness, we include an enumeration of the symmetry properties of the three allowed crystal field doublets in psuedospin- 1  2 pyrochlores in Table III.It fol- lows that complex ring flips are not possible in dipolaroctupolar QSI: a basis exists in which all coupling constants are real.

Geometric language for perturbation theory
The task before us is to simplify terms in the series expansion of H ef f (9).As the calculation is fairly technical, it is useful to develop a graphical language to express terms in the perturbative expansion.To this end, we establish a rigorous correspondence between terms in the perturbation expansion and links (paths) on the diamond lattice: Spinon transport between second neighbor (i.e., same sublattice) diamond sites.We associate this term with the path across links i and j highlighted by a yellow arrow in Fig. 16.
J z± S z i S + j : Spinon transport between nearest-neighbor (i.e., different sublattice) diamond sites, with a S z i dependent Z 2 phase.These are associated with a diamond link j highlighted a purple arrow in Fig. 16.
J ±± S + i S + j : Double-hop of two A site spinons to the same B site (or vice versa), illustrated for example by the green arrows in Fig. 18.
These associations are illustrated in the legend of Fig. 16.

Algorithm for generating terms
We generate terms in the effective Hamiltonian as follows: 1. Fix an order N in perturbation theory.This sets the number of diamond lattice paths available for use.

Complex ring flips
We identify two distinct channels that contribute to the fourth-order complex ring flips, depending on the placement of the two J z± S + terms on the hexagon, which we label ortho and para (see Fig. 16c,d).For the purposes of this section, we define A := J zz,A /2, B := J zz,B /2. .
The term ξ2B may be obtained from interchanging A and B in ξ 2A .

Other terms
In addition to the real and complex ring flips (g+ig ′ )O p discussed in the main text, perturbation theory generates eight-spin operators of the form S z S z O p coupling the ring flip to dangling S z moments nearest neighbor to the sites visited by the ring flip, parameterized by four independent complex couplings ξ1 , ξ2A , ξ2B , and ξ3 as illustrated in Fig. 17 where the sum runs over all pairings of spins i and j that are normal to, but not lying on, the ring.There are four symmetry-inequivalent pairings, here labeled 1, 2A, 2B, 3, corresponding to nearest, second-nearest and third-nearest neighbour normal spins.Note that this expression could have been equivalently written in terms of the radial spins (green dots in Fig. 17b).The equivalence stems from the fact that we work in the 2I2O manifold: if S z n and S z r are, respectively, normal and radial spins belonging to a tetrahedron consisting of spins n, r, 1, 2, the 2I2O constraint may then be used to replace S z r by −S z n − S z 1 − S z 2 .If sites 1 and 2 coincide with a ring-flip, the factors of S z may be replaced by eigenvalues as described in Sec.III.

Relevance of H±±
In principle, the spinon sublattice mixing from H ±± may also contribute to g ′ .However, as we mentioned in the main text, such contributions are higher order than those from J z± and have been neglected in our work.Effective operators within the ice ensemble must annihilate all virtual spinons they create.H ±± makes this task particularly difficult: It creates doubly-charged spinons, necessitating additional terms to return to the spin ice ensemble.The lowest-order term capable of annihilating these spinons by nontrivial loop is of the form H 2 z± H ±± H 2 ± , sketched in Fig. 18.All such terms are higher order in 1/J zz than the ring flips considered in the main text, so they may be justifiably dropped.This presents an issue if we wish to treat these vortex states as 'matter', decoupled from the gauge field -it seems that we can fluidly move between them by a trivial change of coordinates.
There is no such problem in a continuum gauge field.Starting from the closed surface bounding a void v of the diamond lattice, we can arbitrarily subdivide each plaquette into infinitesimally small sub-plaquettes, ultimately converging to the integral which is Dirac quantized in the usual way.In the winding-number language, we subdivide our straight line phases into arbitrarily small segments, ultimately reaching a point where each small plaquette carries an arbitrarily small field.
A partial resolution of this chart dependence on the lattice comes from carefully constructing a coarse-grained continuum limit.Suppose for now that the system contains only zero fluxes.Let us then introduce a defect on some void v using a half-infinite Dirac string, and fix some ϵ > 0 to be a small parameter.There exists a large volume Ω containing the point defect, such that all elementary plaquettes making up its surface ∂Ω satisfy |e| < ϵ, which may be verified by noting that at large separations e ∼ 1/r 2 .This procedure resolves the ambiguity associated with assigning directions to the four phases of Fig. 4: Any large phases are simply subdivided by taking larger Gaussian surfaces until all fields are small.In the limit ϵ → 0 + , the vison charge approaches an integrated field curvature: We have defined n p = ±1 to be orientation consistent with the outward normal of the volume Ω.We recognize this integral to be the usual Chern number [58] of the Dirac monopole.The argument presented here does not generally apply in the presence of multiple visons, which ultimately place a lower bound on ϵ depending on their separation.Indeed, the visons become increasingly illdefined as their mean separation approaches the lattice scale.
Appendix D: Chart-independent static structure factors The Z 2 duality transform R a π/2 C interchanges the phases U (1) 0 ↔ U (1) π/2 , suggesting that the physics of the model should be the same for ϕ g = θ and ϕ g = π/2 − θ.It is however not immediate to see this symmetry in the static electric-sector correlators ⟨ee⟩, illustrated in Fig. 12, as a consequence of having chosen a specific chart, namely the (−π, π) one (see Sec. IV B).
Fig. 20 shows the static structure factors as seen through manifestly chart invariant plaquette correlations, ⟨sin e sin e⟩ and ⟨cos e cos e⟩.Where ϕ g < π/2, we see ⟨sin e sin e⟩ correlations with very similar features to ⟨ee⟩ (Fig. 12).Similarly, when ϕ g > π/2, we see that ⟨cos e cos e⟩ ≃ ⟨(e + π/2)(e + π/2)⟩ reveals (shifted) ⟨ee⟩ correlations.At ϕ g = 0, 0.125π, the ⟨cos e cos e⟩ correlations reveal the Bragg peaks at (000) and (111) expected of a static order on the pyrochlore lattice.Both the (000) and (111) peaks are surrounded by Gaussian diffuse features corresponding to thermal fluctuations.Though the energy scale ρ g = g 2 + (g ′ ) 2 is held constant in this plot, the strength of the fluctuations is controlled by T /g rather than T /ρ g , and the diffuse features are therefore more intense at ϕ g = 0.125.Of course, an identical discussion applies to the sin e correlations at ϕ g = 0.375π, 0.5π.
Yet another view of the phase transition comes from Fig. 21, showing cuts along the line h = 0. Photon correlations are known to vanish exactly along this line; this observable probes only effects due to the visons.The Debye screening of the photons manifests as a Lorentzian centred on (002), of characteristic width set by the Debye length (itself a function of vison density) [28].As one may expect, the vison density is maximal at the phase transition where the field tension T = max{|g|, |g ′ |} is smallest relative to the temperature.Thermally excited, isolated vison dipoles manifest as a uniform background in k-space: indeed, this background is largest at the degeneracy point ϕ g = π/4.
At low temperature (Fig. 21a) we see that traces at ϕ g and π/2 − ϕ g overlap for most of the Brillouin zone, differing only by the presence of a Gaussian diffuse peak at (000).This feature is due to fluctuations in the Bragg peak at (000) in the U (1) π/2 phase.Note that the duality transform R a π/2 C maps e → π/2 − e, preserving all non-singular correlations.At higher temperature (see panel 21b), the chart boundary for e distorts this picture: when a fluctuation in e is large enough to cross the chart boundary at π, the fluctuation manifests in the correlation function as a uniform background, associated with a nearest-neighbour vison dipole [28].This issue is more pronounced in the U (1) π /2 phase than the U (1) 0 one due to the closer proximity of the chart boundary to the mean value of e.

Appendix E: Semiclassical simulation
Our Monte Carlo algorithm is identical to that established by Szabó and Castelnovo [28].The method is based on classical Monte Carlo simulations of the large-S path integral formulation of the U (1) gauge theory, implemented by a Metropolis algorithm consisting of the following stages: MC1 Apply a random gauge transformation to every tetrahedron.
MC2 Attempt to rotate spins about the local z axis, accepting the move with probability ∝ exp(−β∆E) MC3 'Ring flip' step: attempt to change the z components of the six spins around a plaquette by some angle δ in an alternating ring flip pattern.For a hexagon p of spins labelled by 1, 2, 3, 4, 5, 6, this is the transformation The move is either rejected if any S z is moved outside of [−1, 1], or accepted with probability ∝ exp(−β∆E).This step mimics the infinitesimal, large-S action of the O p operator.
Depending on the observable we are attempting to measure, we employ the following two basic processes.a. Annealing.Starting from 1/β > ρ g , we run a fixed number N step of Metropolis steps and then lower the temperature.We use N temp temperatures, which are logarithmically spaced between T hot and T cold .These two variables must be tuned carefully.To simulate a realistic annealing process, one should choose N step to day law implies that ∂ t b, ∇ × e, and J m must all transform the same way under parity and time reversal.This in turn guarantees that e and b have opposite properties under these transformations, i.e., precisely one of the two is a vector and precisely one is even under time reversal.
In conventional electromagnetism, we take e to be a T -even vector, automatically making b a TRS-odd pseudovector.The spin ice literature is divided on this assignment -the convention established by Hermele et al. [1] calls the coarse-grained S z moments, a TRS odd vector, an electric field; the convention established by Castelnovo et al. [2] calls this the magnetic field.Neither convention is ultimately equivalent to Maxwell theory.
The units of (emergent) electomagnetism we work in are Gaussian.Explicitly, the units of electric charge q e and magnetic charge q m are defined such that the visonvison and spinon-spinon interaction energies are as the quasiparticle separation R → ∞.
Appendix G: Large-S perturbation theory in g ′ /g Starting from a system in the U (1) 0 phase, we develop a quantitative understanding of the perturbative effect of g ′ on the system.We will show that the sin e term acts as an anharmonic medium for the photons analogous to a nonlinear crystal.The Hamiltonian of Eq. ( 11) is first expanded about the 0-flux vacuum in falling powers of 1/s = (S + 1/2) −1 , where it is noted that e ∼ S z /s ∼ s−1/2 , and then separated into quadratic and interacting parts: By taking the expansion, we are supposing that e p is in some sense small, and it is therefore safe to assume in this section that e p and ∇ × p a are essentially synonymous.
(Strictly speaking, a chart is being chosen for a in addition to e p .)As noted in the main text (Sec.VI A), V (v) in the large-S limit is a non-dynamical topological term decoupled from the gauge-field dynamics, and it may be set to zero for the purposes of this calculation.The term proportional to gz in Eq. (G1), where z = 6 counts the number of spins on a plaquette, arises from the Villain expansion of the ring-flip [1,6,42].
It can be shown that the next-leading, mixed-operator corrections from g ′ are of the form g ′ (S z /s) 2 e p .In the spirit of presenting a minimal effective theory that captures the high-energy features from semiclassical numerics, such terms have been dropped from H I .These terms essentially amount to a renormalization of the threephoton scattering vertex e 3 already discussed in the main text (Sec.VI A).The partition function of the system can finally be written as: ⟨S z (q, τ )S z (−q, 0)⟩ = 1 Z Da S z (q, τ )S z (−q, 0) 1 + S con-tains only powers of a 3 and higher.Due to the vanishing of odd-powered terms in the action, the leading order in the series is then S I .The Hartree-Fock corrections arising from S (4) I [a] were calculated in Ref. 6.They correct the energy of the existing photon modes and renormalise the speed of light, but do not introduce any additional modes.We shall therefore neglect them hereafter.
It is necessary at this point to move from indexing by links and plaquettes to indexing by tetrahedron and pyrochlore sublattice.We index a spin uniquely by the position R t of its nearest A-type tetrahedron and its pyrochlore sublattice λ ∈ {0, 1, 2, 3}.Similarly, a plaquette may be uniquely indexed by its nearest-neighbour A-void location R v and its plaquette sublattice µ ∈ {0, 1, 2, 3}.It may be readily verified that the spin (R t , λ) is located at the position R t +(1−χ)b λ /2, and the centre of the plaquette (R v , µ) lies at R v −(1−χ/3)b µ /2, where χ ∈ [0, 1) is the breathing anisotropy of the original pyrochlore lattice.For notational simplicity we will set χ = 0, as it has only a marginal effect on the derivation to follow.
Following the standard QED approach [6,41], we introduce bosonic operators associated with each site c Rt,λ , [c Rt,λ , c † Rs,ρ ] = δ ts δ λρ and define their associated Fourier modes to be c λ (k) = 2/N t t∈A exp(−ik • R t )c Rt,λ .Here N t is the total number of tetrahedra, which is equal to the total number of voids.Inspired by QED, we make the definitions where Π = S z /s is now understood as the canonical momentum conjugate to a, and ϵ k is to be determined.The Fourier transforms of e and Π may then be written in terms of the a variables, where After making this substitution, the Wick rotated quadratic action Eq.(G4) reads from which we see that the anomalous terms may be canceled by diagonalizing z −1/2 Z and setting ϵ k to its unique nonvanishing (doubly degnerate) eigenvalue.This fixes ϵ k , recovering the standard photon dispersion.The novel term introduced by the g ′ coupling arises from the third order term in the expansion of sin e,

(G14)
Hidden in the structure of the Γ terms is a symmetryenforced Kronecker delta between µ and ρ.We plot (the trace of) this this term in Fig. 11, with an arbitrary intensity scale.

FIG. 1 .
FIG.1.Illustration of the pyrochlore and diamond lattices.Gray tetrahedra share corners which form the pyrochlore sites i (red sphere, right), which are in one-to-one correspondence with the links (red line, left) of the diamond lattice.Each diamond plaquette p (green even-sides hexagon, left) is similarly in unique correspondence with the hexagonal plaquettes of the pyrochlore lattice (green alternating-sided hexagon, right).Four pairwise adjacent plaquettes enclose a void v (purple region, bottom left).

5 FIG. 3 .
FIG. 3. Four of the individual terms contributing to Fig. 2a.The curved purple lines denote a factor of ζijS + i S z j or ζ * ij S − i S z j .As indicated in the main text, a), b) lead to complex ring flips whereas c) leads to real ring flips.The term d) cancels with another term in the expansion and does not contribute.

FIG. 4 .
FIG.4.Top panel: Exact vanishing of ∇ •v ∇ ×p a on a particular lattice void, composed of the four plaquettes p1,2,3,4.One may regard the link variables a as either U (1) rotors or real numbers.Bottom panel: Summation of phases involved in V (v) (Eq.(19)) as viewed in the (−π, π) chart.The topology of U (1) is depicted using a circle, with the four red arrows representing the four ep fluxes associated with the void v.As discussed in the main text, the path must return to its starting point to remain consistent with ∇ •v ∇ ×p a = 0(mod1).The removal of ±π phases forbids any red arrow to enter the midpoint of the circle -this is in effect a topological singularity about which a winding number can be defined.
represents a transformation acting on the spin model which rotates the phase of the ring flip operator on any plaquette p by O p → e i∇×pa O p .The π-flux configuration a π (j) of Fig. 5 (right panel) results in R aπ [O] = −O, reproducing the well known duality [46, 48] between the U (1) 0 and U (1) π phases.The left panel of the same figure is a 'square root' of the U (1) π ↔ U (1) 0 transformation, generating R a π/2 [O p ] = iO p for all 16 plaquettes in the unit cell.

FIG. 6 .
FIG. 6. Phase diagram of the g + ig ′ = ρg exp(iϕg) model from semiclassical numerics, showing a) the average expectation value of the large-S average vison charge VA on A tetrahedra, b) the energy of the system relative to the 0-flux state, c) and d) the expectation values ⟨sin e⟩ = p sin ep/Np and ⟨cos e⟩ = p cos ep/Np, where Np denotes the total number of plaquettes.The black, dotted curve Tpm/Tc = √ 2 [1.3 max(| cos(ϕg)|, | sin(ϕg)|) − 0.3 min(| cos(ϕg)|, | sin(ϕg)|)], is given as a guide to the eye, outlining the crossover between paramagnetic and U (1) behaviour.The resolution of ϕg is ∼ π/100 in the region |ϕg − π 4 | < 0.1, and π/19 elsewhere, with data points indicated by the x-axis ticks.
overestimates the Ising value by a factor of 3, well beyond what can be explained by our fitting uncertainty.It shall remain an intriguing open question for future work.

CLFIG. 8 .
FIG. 8. Finite-size scaling collapse at the critical end point of the g = g ′ line, in terms of τ = (T − Tc)/Tc, for systems consisting of L 3 cubic unit cells.a) Specific heat capacity, var[E].b) Reduced thermally averaged staggered vison order parameter, ⟨VA⟩ − ⟨VA⟩c, where ⟨...⟩c denotes the thermal average at the critical point.c) Susceptibility of the order parameter, var[VA].We performed a fit for Tc and the relevant scaling exponent for each of the three quantities separately, with results given directly in each panel.For reference, the critical exponents of the 3D Ising model are ν = 0.629971(4), α = 0.11008(1), β = 0.326419(3), γ = 1.237075(10) [52, 53].

2 S
FIG.11.Two-photon continuum in the g, g ′ model.a) Contribution from the g ′2 term to the self energy (arb.intensity units), calculated analytically as discussed in the main text and Appendix G. b) Monte Carlo (MC) numerical results from Fig.9, for g ′ /g = 0.8.Both panels correspond to T = 0.1g.

FIG. 14 .
FIG.14.Scaling of vison interaction energy with system size within the U (1)0 phase.Panel a) shows the fits of the Coulomb functional form Uv = MZnS q 2 e /(La0) + µ to the data points with L ≥ 17.This model marginally overestimates the cohesive at L = 5 (see the data points outlined by a solid black line box).Note that qe may be interpreted as the elementary vison charge in Gaussian CGS units.The inset shows the diamond superlattice arrangement of visons, see also Fig.13for details.The lower panels show the dependence of b) the Coulomb constant q 2 e and c) the chemical potential µ on g ′ .Panel b) shows the dependence of the q 2 e fits on the short-range interaction cutoff, suggesting convergence to the quadratic estimate q 2 e = πg (horizontal dashed black line).Panel c) shows the dependence of the vison chemical potential µ on g ′ , where a linear best fit µ(g, g ′ )/g = −1.1154(2)g′ /g + 7.873(1) appears to be in excellent agreement with the simulations.These results are independent of the cutoff choice.For all plots, the absence of error bars indicates uncertainties smaller than the marker size (for reference, the standard error on the µ fits is of the order of 10 −4 g).

45 LFIG. 15 .
FIG.15.Dependence of the vison interaction energy on g ′ .a) Scaling of the empirical vison interaction energy Uv − ∂µ ∂g g − ∂µ ∂g ′ g ′ with g ′ .b) Log-log plot of the gradient of the interaction part of the vison energy with respect to g ′ against linear system size L, exhibiting a potential crossover from L −2 to L −1 scaling.The error bars reflect 1σ standard error from linear regression.

Appendix B :
Derivation of the Complex ring flip : p

FIG. 17
FIG. 17. a) Normal S z operators used in the definitions of the ring flip terms appearing in Eq. (B7).b) A pyrochlore plaquette (black dots) and its nearest-neighbor off-ring sites (red and green dots), showing the four symmetry-inequivalent spin pairings.Red nearest neighbours correspond to those shown in panel a).The locations of the S z sites of the four effective interactions ξ1, ξ2A,B, and ξ3 are also indicated in panel b)

FIG. 19 .
FIG.18.An example of leading-order contribution of H±± to complex ring flip processes.

TABLE II .
Careful definitions of the (classical) coordinate quantities used in QSI (see Fig.1).

TABLE III .
Known crystal field doublets in rare-earth pyrochlores under the action of time reversal T and rotation C3 about the local z axis.Rotation C2 about the local x axis gives S z → −S z and S ± → S ∓ for all cases.Here w = e 2πi/3 ; and ζij, γij, and Jz± are defined in Eq. (1).Table adapted from Ref. 75.