Seeing beyond the light: Vison and photon electrodynamics in quantum spin ice

TL;DR

This paper introduces a semiclassical numerical method to study excitations in quantum spin ice, revealing visons as a weak electrolyte and their hybridization with photons, advancing understanding of quantum spin liquids.

Contribution

The paper presents a new semiclassical approach to analyze gapped excitations in quantum spin liquids, specifically characterizing visons and their thermodynamic behavior.

Findings

Visons form a weak electrolyte with low-temperature pair dominance.

Visons strongly hybridize with photon backgrounds, affecting experimental signatures.

The method effectively studies thermodynamic and spectral properties of quantum spin ice.

Abstract

Understanding the nature and behaviour of excitations in quantum spin liquids, and in topological phases of matter in general, is of fundamental importance and has proven crucial for experimental detection and characterisation of candidate materials. Current theoretical and numerical techniques, however, have limited capabilities, especially when it comes to studying gapped excitations. Here, we propose a semiclassical numerical method to study systems whose spin liquid behaviour is underpinned by perturbative ring-exchange Hamiltonians. Our method can readily access both thermodynamic and spectral properties. We focus in particular on quantum spin ice and its photon and vison excitations. After benchmarking the method against existing results on photons, we use it to characterise visons and their thermodynamic behaviour, which remained hitherto largely unexplored. We find that visons, in contrast to spinons in classical spin ice, form a weak electrolyte: vison pairs are the dominant population at low temperatures. This is reflected in the behaviour of thermodynamic quantities, such as pinch point motifs in the relevant correlators. Visons also appear to strongly hybridise with the photon background, a phenomenon that affects the way these quasiparticles may show up in inelastic response measurements. Our results demonstrate that the method, and generalisations thereof, can substantially help our understanding of quasiparticles and their interplay in quantum spin ice and other quantum spin liquids, quantum dimer models, and lattice gauge theories in general.

Equations

$$\mathcal{S}=S\sum _i{\int }_0^T\mathrm{d}t{\int }_0^1\mathrm{d}s{\mathbf{n}}_i\cdot ({\partial }_t{\mathbf{n}}_i)\times ({\partial }_s{\mathbf{n}}_i)-{\int }_0^T\mathrm{d}tH[\mathbf{n}(t)],$$

$$S{\mathbf{n}}_i\times ({\partial }_t{\mathbf{n}}_i)=\frac{\partial H[\mathbf{n}]}{\partial {\mathbf{n}}_i}$$

$${\partial }_t{\mathbf{\sigma }}_i=\frac{\partial H[\mathbf{\sigma }]}{\partial {\mathbf{\sigma }}_i}\times {\mathbf{\sigma }}_i.$$

$$H=J\sum_{\langle ij\rangle}\sigma_{i}^{z}\sigma_{j}^{z}=\frac{J}{2}\sum_{t}\bigg{(}\sum_{i\in t}\sigma_{i}^{z}\bigg{)}^{2}+\mathrm{const.},$$

$$H=\sum _{⟨ij⟩}J{\sigma }_i^z{\sigma }_j^z-\frac{J_{\perp }}{2}({\sigma }_i^{+}{\sigma }_j^{-}+{\sigma }_i^{-}{\sigma }_j^{+}),$$

$$H_{\mathrm{eff}}=-\frac{g}{2}\sum _{\text{hexagon}}\left({\sigma }_1^{+}{\sigma }_2^{-}{\sigma }_3^{+}{\sigma }_4^{-}{\sigma }_5^{+}{\sigma }_6^{-}+H.c.\right),$$

$$A_{ab}\to A_{ab}+{\chi }_a-{\chi }_b,$$

$${\mathbf{h}}_i$$

$$f({\varphi }_i)\propto e^{\beta \mathrm{Re}({\sigma }_i^{+}{h}_i^{\ast })}=\exp \left[|{\sigma }_i^{+}||h_i|\cos ({\varphi }_i-\arg h_i)\right],$$

$$B=\mathrm{curl}A=\arg ({\sigma }_1^{+}{\sigma }_2^{-}{\sigma }_3^{+}{\sigma }_4^{-}{\sigma }_5^{+}{\sigma }_6^{-});$$

$$S^{zz}(\mathbf{q},\omega)=\sum_{\mu}\big{\langle}\sigma^{z}_{\mu}(\mathbf{q},\omega)\sigma^{z}_{\mu}(-\mathbf{q},-\omega)\big{\rangle}$$

$$q=\mathrm{div}B/2\pi$$

$$\epsilon (L)=\mu -\frac{M\alpha }{L{a}_0}$$

$$\mu$$

$$\langle B(-\mathbf{q})B(\mathbf{q})\rangle=C_{1}+C_{2}\sin^{2}(q_{z}/8)+\frac{C}{1+\big{[}8\xi\cos(q_{z}/8)\big{]}^{2}}$$

$${\mu }_E$$

$${\mu }_{Arrh.}$$

$$(\mathrm{grad}U{)}_{\mathbf{r}\to {\mathbf{r}}^{\prime }}$$

$$(\mathrm{div}A{)}_{\mathbf{r}}$$

$$B=-\mathrm{grad}U+\mathrm{curl}A+C^{\mu },$$

$$E-E_0$$

$$\mathrm{div}B=-\mathrm{div}\mathrm{grad}U=2\pi n.$$

$$(\mathrm{grad}U{)}^{\mu }(\mathbf{q})$$

$$(\mathrm{div}B{)}^{\alpha }(\mathbf{q})$$

$$\mathbf{M}(\mathbf{q})=\left(\begin{array}[]{cccc}e^{i\mathbf{q}\cdot\mathbf{r}_{0}/2}&e^{i\mathbf{q}\cdot\mathbf{r}_{1}/2}&e^{i\mathbf{q}\cdot\mathbf{r}_{2}/2}&e^{i\mathbf{q}\cdot\mathbf{r}_{3}/2}\\ -e^{-i\mathbf{q}\cdot\mathbf{r}_{0}/2}&-e^{-i\mathbf{q}\cdot\mathbf{r}_{1}/2}&-e^{-i\mathbf{q}\cdot\mathbf{r}_{2}/2}&-e^{-i\mathbf{q}\cdot\mathbf{r}_{3}/2}\end{array}\right),$$

$$\mathbf{M}{\mathbf{M}}^{†}\mathbf{U}$$

$$\mathbf{U}$$

$$\mathbf{B}={\mathbf{M}}^{†}\mathbf{U}$$

$$E-E_0$$

$$\displaystyle\simeq\frac{(2\pi)^{2}g}{2N}\sum_{\mathbf{q}\neq 0}\mathbf{n}^{\dagger}\big{(}\mathbf{M}\mathbf{M}^{\dagger}\big{)}^{-1}\mathbf{M}\mathbf{M}^{\dagger}\big{(}\mathbf{M}\mathbf{M}^{\dagger}\big{)}^{-1}\mathbf{n}$$

$$\displaystyle\longrightarrow 2\pi^{2}gV_{\mathrm{cell}}\int_{\mathbf{q}\in\mathrm{BZ}}\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\mathbf{n}^{\dagger}\big{(}\mathbf{M}\mathbf{M}^{\dagger}\big{)}^{-1}\mathbf{n}\,,$$

$$\mathbf{M}{\mathbf{M}}^{†}$$

$$\displaystyle\big{(}\mathbf{M}\mathbf{M}^{\dagger}\big{)}^{-1}$$

$$E-E_0$$

Full text

I Introduction

Quantum spin liquids (QSL) are topological phases of frustrated magnetic materials in which quantum fluctuations prevent magnetic order even at zero temperature Savary and Balents (2017); Balents (2010). This phenomenon is often accompanied by exotic behaviour, including emergent gauge symmetries and fractionalised quasiparticle excitations. Understanding and characterising such systems has attracted substantial interest in recent years. A particularly timely and important problem is identifying experimentally relevant signatures of QSL behaviour, in light of several candidate materials that have been proposed of late, but remain to be confirmed and characterised.

Many such phases lend themselves to gauge theoretic descriptions that are reflected in their quasiparticle content, which can in turn give rise to salient features in dynamical spin structure factors and other response properties Knolle and Moessner (2019). Understanding these features is both of fundamental importance and a promising diagnostic tool for QSL materials. Indeed, studying the dynamical structure factor in the Kitaev honeycomb model Kitaev (2006), for example, has recently led to compelling evidence of a possible QSL phase in $\alpha$-RuCl3 in a magnetic field Banerjee et al. (2016, 2017); Do et al. (2017); Banerjee et al. (2018); Kasahara et al. (2018).

In general, however, response and equilibration properties of quantum spin liquids are far from being well understood. Analytical results are a tall order away from exactly solvable models and numerical techniques that access the excitation spectrum of strongly correlated many-body systems in two and three dimensions have limited capability, especially when facing systems with gapped quasiparticles that are well defined only in a coarse-grained sense.

A fruitful avenue that has been used to study quantum magnets, including some QSLs, proceeds by approximating quantum spins as classical ones (see for instance Refs. Moessner and Chalker (1998); Conlon and Chalker (2009, 2010); Taillefumier et al. (2014); Robert et al. (2015); Taillefumier et al. (2017); Robert et al. (2008); Schnabel and Landau (2012) on pyrochlore and kagome antiferromagnets and Refs. Baskaran et al. (2007, 2008); Samarakoon et al. (2017, 2018) on Kitaev honeycomb systems). Removing quantum correlations and uncertainties allows performing Monte Carlo sampling and time evolution more efficiently as well as accessing a wider variety of observables than usually possible with fully quantum methods. The starting points of these studies are nearest neighbour bilinear exchange Hamiltonians. While this is a successful approach for certain spin liquid phases, QSL behaviour is often underpinned by perturbative ring-exchange processes that induce matrix elements between classically degenerate spin configurations (e.g., in quantum spin ice Hermele et al. (2004) and resonant valence bond phases Rokhsar and Kivelson (1988); Moessner and Sondhi (2001)). A direct semiclassical representation of the original Hamiltonian in these cases results in classical phases that do not capture the QSL physics – this renders the approach inapplicable. Indeed, retaining spin liquid behaviour in both semiclassical simulations and large-$S$ path integral calculations depends crucially on a connected continuum of low-energy states between which low-temperature dynamics remains possible.

In this paper, we develop a semiclassical numerical method to investigate quantum spin liquids that hinges upon a large-$S$ description of the effective ring exchange Hamiltonian. Our approach effectively samples the finite-temperature path integral formulation of the problem using classical Monte Carlo and dynamical simulations to obtain both thermodynamic and spectral properties, thus overcoming previous limitations of semiclassical simulations Taillefumier et al. (2017); Onoda and Tanaka (2011).

We demonstrate the validity and capability of our approach by studying a class of highly anisotropic quantum spin liquids on the pyrochlore lattice called quantum spin ice (QSI) Gingras and McClarty (2014). QSI phases are described by a compact $U(1)$ gauge theory analogous to lattice quantum electrodynamics (QED) Hermele et al. (2004). Similarly to ordinary QED, QSI exhibits linearly dispersing gapless photon modes, as well as gapped electric charges (spinons) Note (11). 1111footnotetext: There is some ambiguity around the names of emergent fields in the literature. Following the lattice gauge theory convention, we call the emergent vector potential magnetic, so that the description matches the terminology of standard QED. However, this makes spinons sources of electric field, which is at odds with them being sources of the physical magnetic field in dipolar spin systems and real materials. Furthermore, the compact nature of the theory allows for Dirac quantised Dirac (1931) gapped magnetic monopoles [$U(1)$ visons]. While no material has conclusively been classified as QSI to date, evidence points to, for instance, several praseodymium-based pyrochlores as promising candidates Rau and Gingras (2019).

This picture of QSI is consistent with quantum Monte Carlo (QMC) studies Benton et al. (2012); Shannon et al. (2012); Kato and Onoda (2015); Huang et al. (2018); Banerjee et al. (2008); Lv et al. (2015), but the current understanding of its gapped quasiparticles is limited. Very few experimentally relevant signatures of visons have been proposed Chen (2016, 2017), and they remain elusive to QMC studies, which often suffer from the sign problem Shannon et al. (2012); Kato and Onoda (2015) and have only limited capability to obtain excitation spectra and dynamics Huang et al. (2018). Rigorous analytic treatment of spin-$1/2$ systems is a tall order in all but the simplest cases, and the gauge theory picture is based on analytical soft-spin Hermele et al. (2004); Benton et al. (2012) or large-$S$ Kwasigroch et al. (2017) expansions which in turn struggle to access gapped excitations quantitatively Kwasigroch . Mean field and slave boson approaches Savary and Balents (2012); Lee et al. (2012); Savary and Balents (2013); Hao et al. (2014) as well as numerical linked cluster calculations Benton et al. (2018) have been successful in capturing potential phases exhibited by QSI Hamiltonians and transitions between them; however, they have not been equally informative about their excitation spectra. Finally, exact diagonalisation for three-dimensional strongly correlated quantum systems is limited to very small system sizes Onoda and Tanaka (2011); Udagawa and Moessner (2019). New tools to study excitations in QSI and related systems are therefore in high demand.

Semiclassical techniques have been applied to the archetypal bilinear QSI model Taillefumier et al. (2017); however, these models order in a $\sigma^{z}$-polarised classical spin ice (CSI) state in the QSI limit where the Ising term is dominant Taillefumier et al. (2017); Onoda and Tanaka (2011). By taking perturbative ring-exchange processes Hermele et al. (2004) into account explicitly, our method is able to capture the spin liquid phase of QSI. As a benchmark, we derive the dispersion of its photonic modes and show that it is in excellent agreement with the prediction of large-$S$ field theory Kwasigroch et al. (2017), and with QMC results on the original spin-1/2 system Benton et al. (2012); Huang et al. (2018).

The majority of the paper is devoted to novel results on the $U(1)$ visons that remained elusive in previous studies of quantum spin ice. We consider both the energetics of zero-temperature metastable vison states, as well as thermal ensembles of visons and photonic excitations. We are able to obtain their bare energy cost and long range Coulomb interaction, borne out of the semiclassical equivalent of the quantum kinetic energy in a purely short ranged Hamiltonian. Contrary to spinons, both the energy cost and the interaction strength of visons are controlled by the same energy scale and therefore their relative strength is fixed. We find that visons are in the weak electrolyte limit, where their interaction is strong enough to make nearest neighbour pairs energetically favourable over isolated visons. This has important consequences, for instance, if one aims to develop effective models for visons in QSI Chen (2017, 2016); Kwasigroch , since in thermodynamic equilibrium, dilute isolated visons only occur in a relatively dense environment of vison pairs. We demonstrate that the weak electrolyte behaviour is reflected in thermodynamic properties of the system, such as the blurring of pinch points of magnetic field correlators as a function of temperature. Importantly, we also observe a strong interplay between isolated visons, the aforementioned vison pair plasma, and photons, which significantly affects the thermodynamic vison density. This is possibly a semiclassical reflection of quantum hybridisation of photon and vison excitations in QSI, which would have a substantial impact on the possibility to detect visons experimentally in inelastic response measurements. More specifically, the dressing of isolated visons by vison pairs suggests an intriguing analogy with particle–antiparticle bubbles in the strong coupling problem in QED – an aspect that certainly warrants investigating in future work.

The present work can readily be extended to investigate photon and vison dynamics, e.g., in relation to magnetic noise and transport measurements, as well as to include spinon excitations. Moreover, our approach is not limited to quantum spin ice systems and can be straightforwardly generalised to other QSLs underpinned by perturbative ring exchange processes, including valence bond and quantum dimer models (e.g., following the route proposed in Ref. Szabó et al. (2019)).

The rest of the paper is organised as follows. We discuss our method in general terms in Sec. II. Section III introduces the QSI model and describes how the method was applied to it. Benchmarking results on photon modes are presented in Sec. IV. Metastable vison configurations are studied at low photon densities in Sec. V, while Sec. VI deals with thermodynamic properties of the semiclassical QSI model in the presence of both excitations. Conclusions are drawn in Sec. VII.

II Semiclassical simulation of quantum spin liquids

Consider a Hamiltonian that remains in a QSL phase in the limit of large $S$, that is to say, its eigenstates are massively entangled in any local basis Savary and Balents (2017). In a path integral representation, this entanglement corresponds to the interference of a continuum of equivalent trajectories, related to each other by gauge symmetry Auerbach (1994). The interference itself is a defining feature of spin-1/2 QSLs, since it accounts for differences between related classical and quantum spin liquids Benton et al. (2012). For large $S$, however, quantum fluctuations become unimportant and interference effects are only apparent at the lowest temperatures. Therefore, a large-$S$ QSL is generally indistinguishable from a classical spin liquid (CSL) characterised by a massive degeneracy of the least action trajectories of the quantum path integral.

It is now straightforward to obtain static correlation functions of the large-$S$ QSL by Monte Carlo sampling the CSL Boltzmann distribution given by $e^{-\beta H[\mathbf{n}]}$, with the formally identical Hamiltonian $H[\mathbf{n}]$ understood to act on unit vectors ${\mathbf{n}}\equiv\{{\mathbf{n}}_{i}\}$ rather than quantum spins. Such sampling can be done more efficiently than quantum Monte Carlo and never suffers from a sign problem; furthermore, it naturally captures gapped excitations which are hard to treat in analytic large-$S$ calculations Kwasigroch et al. (2017); Kwasigroch .

The time evolution of the QSL is described by a real time path integral with action

[TABLE]

where the first term is the Berry phase for real time evolution of spins Auerbach (1994), and $\mathbf{n}_{i}$ is a unit vector specifying the spin coherent state. As discussed above, this path integral can be regarded in the large-$S$ limit as an incoherent superposition of independent least action trajectories. For each of these, varying the action (1) gives

[TABLE]

or, equivalently,

[TABLE]

In the second line, we used the fact that $\mathbf{n}_{i}^{2}=1$ at all times and thus $\mathbf{n}_{i}\cdot\partial_{t}\mathbf{n}_{i}=0$ Auerbach (1994), and we re-expressed the equation in terms of the expectation values of the spin operators $\boldsymbol{\sigma}_{i}=S\mathbf{n}_{i}$, as customary. Equation (2) can also be derived from Ehrenfest’s theorem using $[\sigma^{\mu},\sigma^{\nu}]=i\varepsilon^{\mu\nu\lambda}\sigma^{\lambda}$; the semiclassical equations follow by replacing operators with their expectation values.

Simulating the large-$S$ dynamics of the system now involves solving the differential equation (2); however, with no quasiparticles and infinitesimal zero point fluctuations, that solution would be trivial at $T=0$. Using finite temperature Monte Carlo configurations as initial conditions of time evolution is a natural prescription for finding dynamic correlators of a classical spin liquid; in the large-$S$ path integral language, this essentially samples least-action trajectories of a finite-temperature (e.g., Keldysh) path integral.

Finally, it is important to understand how QSL quasiparticles might appear in semiclassical simulations. These are often due to canonical quantisation – for example, emergent photons are quanta of lattice electromagnetic modes and integer spinon numbers are set by the quantisation of angular momentum. In the semiclassical limit, such quantisation is irrelevant, but while individual quasiparticles disappear, we anticipate that equivalent classical normal modes survive and their frequency dispersion remains indicative of the original quasiparticle. By contrast, certain quasiparticles [e.g., $U(1)$ visons in QSI] are due to $2\pi$-ambiguities of phases: Since angles remain quantised even in a semiclassical setting, such particles survive as gapped excitations in our simulations.

III Quantum spin ice

In the rest of the paper, we focus on quantum spin ice (QSI), a QSL model defined on the pyrochlore lattice (see Fig. 1). It is based on the classical spin ice (CSI) model of Anderson Anderson (1956),

[TABLE]

where $t$ runs over the corner-sharing tetrahedra of the pyrochlore lattice. As manifest from the second form, (3) is minimised if the total $\sigma^{z}$ around each tetrahedron is zero. This is analogous to Pauling’s “ice rules” for water ice and gives rise to extensive ground state degeneracy and long-range correlated disorder Bramwell and Gingras (2001). The elementary excitations of the model are individual tetrahedra violating the ice rules that take the form of fractionalised deconfined spinons Castelnovo et al. (2012).

The classical spin ice model (3) is an Ising model and, as such, its quantum dynamics is trivial. Off-diagonal matrix elements between the classical ground states can be introduced by transverse interactions. Perhaps the simplest such model is the anisotropic XXZ Hamiltonian Hermele et al. (2004)

[TABLE]

where $J_{\perp}\ll J$. The ground state of (4) is expected to be a superposition of two-in-two-out spin states, as the Ising term remains dominant; since there is an extensive number of such states, this ground state is massively entangled.

In the no-spinon sector, the transverse term only acts perturbatively, by flipping closed loops of spins: This can be thought of as creating a virtual pair of spinons and carrying them around the loop. Since the shortest closed loops in the pyrochlore lattice are hexagonal plaquettes, the lowest-order nontrivial perturbative Hamiltonian is

[TABLE]

where $g=3J_{\perp}^{3}/J^{2}$ Hermele et al. (2004).

Using soft-spin Hermele et al. (2004) or large-$S$ Kwasigroch et al. (2017) approximations, the ring exchange Hamiltonian (5) can be mapped onto a compact $U(1)$ lattice gauge theory, where $\sigma^{z}$ takes the role of the electric field $E$ and the phase of $\sigma^{+}$ that of the vector potential $A$ Note (11). In particular, the Hamiltonian is invariant under the gauge transformation

[TABLE]

where the $\chi$ are arbitrary angles for each tetrahedron and the subscript $ab$ denotes the spin belonging to neighbouring tetrahedra $a$ and $b$. This mapping indicates that quantum spin ice hosts gapless photon excitations, while spinons become quantised, gapped sources of the electric field. Furthermore, the compact gauge theory also supports $2\pi$-quantised sources of the emergent magnetic field $B=\operatorname{curl}A$, since $B$ itself is only defined modulo $2\pi$. These magnetic charges are called $U(1)$ visons: Like spinons, they are gapped [by an $O(g)$ gap] and interact through a Coulomb interaction mediated by photons Hermele et al. (2004). This picture has partly been confirmed, among others, by quantum Monte Carlo studies of the spin-1/2 model Benton et al. (2012); Shannon et al. (2012); Kato and Onoda (2015); Huang et al. (2018); Banerjee et al. (2008); Lv et al. (2015); however, current analytical and numerical techniques suffer from shortcomings when it comes to investigating its gapped excitations.

It is important to note that the perturbative analysis leading to (5) is only valid for $S=1/2$, where diagonal perturbations to the Hamiltonian are equal for all classical ground states. For $S>1/2$, these terms lift the degeneracy of classical ground states and thus disrupt the spin ice phase. As $S\to\infty$, the Ising and transverse terms favour easy-axis and easy-plane arrangements, respectively, leading to energetic competition: In the limit $J_{\perp}\ll J$, the Ising term wins and the system settles in a classical spin ice state with small transverse spin components Taillefumier et al. (2017); Onoda and Tanaka (2011). These states are stable against transverse fluctuations, removing any dynamics analogous to QSI. The ring exchange Hamiltonian (5), however, remains in a QSL phase for arbitrary $S$ Hermele et al. (2004); Kwasigroch et al. (2017): Since it only contains transverse terms, its classical ground state is easy-plane and the electric field appears as small fluctuations of $\sigma^{z}$, which are amenable to continuous time evolution (cf. Fig. 1).

Using (2), we developed semiclassical equations of motion for the ring exchange Hamiltonian (5). These can be written in the familiar Larmor form $\dot{\boldsymbol{\sigma}}_{i}=\boldsymbol{\sigma}_{i}\times\mathbf{h}_{i}$ with effective field

[TABLE]

where the summation is over the six plaquettes that $i$ belongs to (see Fig. 1) and $i+n$ denotes the $n$th spin counted from $i$ on each plaquette (the direction around the plaquette is immaterial). These equations were integrated numerically using the GNU Scientific Library implementation of the Prince–Dormand (8,9) ODE solver with automatic step size control Galassi et al. (2009).

In order to generate thermally distributed classical spin configurations, we use a Monte Carlo algorithm that samples $\sigma^{z}$ and the phase of $\sigma^{+}$ (for its magnitude is fixed by $\sigma^{z}$) independently – this mirrors the anisotropy of the Hamiltonian (5). Furthermore, we insist that there be no spinons in the system, that is, $\sum\sigma^{z}=0$ for all tetrahedra of the pyrochlore lattice. This can be achieved by updating $\sigma^{z}$ only in closed loops with alternating signs, similarly to typical low-temperature simulations of CSI Barkema and Newman (1998); Melko and Gingras (2004). Unlike CSI, however, $\sigma^{z}$ is now a continuous variable, so loops can always be updated by small amounts. For convenience, we only perform updates around hexagons 313131 One can see that the hexagon updates lead to an ergodic Monte Carlo protocol within each $\sigma^{z}$ magnetisation sector as follows. The no-spinon constraint implies that $\sigma^{z}$ is the sum of a pure lattice curl and a global magnetisation (cf. Appendix B). Therefore, for any two valid configurations of $\sigma^{z}$ with equal overall magnetisation, there is an “electric vector potential” $G$ on the dual pyrochlore lattice such that $\operatorname{curl}G=\Delta\sigma^{z}$: Changing $\sigma^{z}$ around each hexagon by this $G$ on the corresponding dual pyrochlore site takes the system from one configuration to the other. In our work, we focus on the entropically dominant zero magnetisation sector only since none of the properties we look at depend on changes in the macroscopic magnetisation of the system. If needed, loop updates to sample different sectors could easily be introduced in the algorithm to remove this limitation. ; the proposed change is drawn from a Gaussian distribution whose variance is proportional to temperature (this ensures a large acceptance rate at all temperatures); and acceptance is decided using the Metropolis method.

On the other hand, there are no conservation laws to be obeyed by the phase $\phi$ of $\sigma^{+}$, so it can be updated independently for each spin. In particular, since the magnitude of $\sigma^{+}$ is kept constant, its distribution is given by

[TABLE]

where the effective exchange field $h_{i}$ is given by (7). This is a von Mises distribution which can be efficiently sampled directly Best and Fisher (1979); Devroye (1986), that is, $\phi_{i}$ can be sampled without rejection. Furthermore, to eliminate spurious correlations due to not sampling the gauge freedom (6), each Monte Carlo step includes rotating $\phi$ for the four spins of each tetrahedron by an angle $\chi$ drawn from a uniform distribution.

It is important to note that, in the semiclassical picture, quantum correlations and uncertainties are irrelevant and therefore all quantities can be represented by pure numbers rather than quantum operators. That is, given a configuration of spins, the value of any observable can be determined unambiguously and straightforwardly. For intricate quantities such as the magnetic field $B$, this is a substantial improvement over standard methods, e.g., quantum Monte Carlo. In particular, we follow the soft spin prescription to identify $\sigma^{\pm}$ with $e^{\pm iA}$ Hermele et al. (2004) and take the vector potential $A$ to be the complex argument of $\sigma^{+}=\sigma^{x}+i\sigma^{y}$. Now, the magnetic field $B$ follows as

[TABLE]

the argument function is restricted to the interval $[-\pi,\pi)$ for a unique $B$ with the smallest possible modulus, consistently with Ref. Hermele et al. (2004). Furthermore, the choice of $\sigma_{1}$ for each plaquette affects the sign of $B$; this choice must be made with reference to the pyrochlore lattice geometry to obtain a self-consistent lattice vector field (see also Ref. Benton et al. (2012) and Appendix A).

IV Photons

We demonstrate the validity and benchmark the accuracy of our method by investigating photon modes in the low-temperature dynamics of the system. Indeed, clean photon modes arise if their interaction with each other and with other excitations is minimised. Spinons are excluded altogether by the loop update Monte Carlo algorithm; visons and photon–photon interactions can be eliminated by reducing the temperature.

We generated 4096 stochastically independent Monte Carlo configurations of a cubic sample of size $L=24$ unit cells at temperature $T=10^{-4}g$ and calculated the time evolution of each for 2048 time steps of size $\delta t=(16g)^{-1}$. The results were Fourier transformed using the FFTW library Frigo and Johnson (2005) in time and space, separately for the four fcc sublattices of the pyrochlore lattice. Following Ref. Huang et al. (2018), we evaluated the correlator

[TABLE]

along high symmetry directions, where the summation runs over the sublattices; the results are plotted in Fig. 2. A single set of remarkably sharp normal modes appear in the data 414141The period of oscillations is not necessarily commensurate with the simulation time window. After Fourier transforming, this leads to the slight broadening of the dispersion seen in Fig. 2.. The frequencies of the numerically obtained modes match perfectly with analytic results for the large-$S$ photon dispersion Kwasigroch et al. (2017) (green line in Fig. 2), confirming that the CSL simulated by our method is indeed equivalent to large-$S$ QSI.

We note that the integrated structure factor $S^{zz}(\mathbf{q})=\int\!\mathrm{d}\omega\,S^{zz}(\mathbf{q},\omega)$ is independent of $\mathbf{q}$. This contrasts spin-1/2 QSI where $S^{zz}(\mathbf{q})\propto\omega(\mathbf{q})$ Benton et al. (2012). This discrepancy is caused by the different realisation of photons in the two systems. In spin-1/2 QSI, zero point fluctuations of photon modes of energy $\hbar\omega(\mathbf{q})/2$ give rise to dynamic correlators even at zero temperature; in the large-$S$ classical spin liquid, photons are classical normal modes with energy determined by equipartition, leading to $\mathbf{q}$-independent correlators 494949In terms of a large-$S$ expansion, we have set the magnitude of our spins, $\hbar S$, to 1. Therefore, $\hbar\sim 1/S$ and thus the energy of individual photons, $\hbar\omega(\mathbf{q})$, also scales as $1/S$. In our simulations, $S\to\infty$ and so $T\gg\hbar\omega$ at any nonzero temperature. Photon populations are thus always large and can be thought of as being in a coherent, classical state..

V Energetics of bare visons

$U(1)$ visons in QSI are $2\pi$-quantised sources of the emergent magnetic field $B$. Their existence and quantisation is due to the $2\pi$-ambiguity of the transverse phases $A$ that are promoted to a vector potential in the gauge theory description. For the same reason, however, specifying vison numbers unambiguously is far from trivial. Typically, visons are understood through their far-field effects, where $B\sim 1/r^{2}$ is small and thus well defined, and the total flux across a large, closed surface gives a unique vison charge Savary and Balents (2017). In principle, one can define a $U(1)$ vison charge operator

[TABLE]

for each dual diamond site by giving $B$ on each plaquette a unique value (e.g., by restricting it to between $-\pi$ and $\pi$). However, since visons are dynamical, the system will not normally be in an eigenstate of $q$, which makes pinpointing visons complicated. A great advantage of the semiclassical method is that observables like $B$ and $q$ are pure numbers rather than quantum operators: This means that the vison charge as defined above is always an unambiguous integer for all dual diamond sites.

The semiclassical simulation also allows us to introduce visons at will. We achieve this via a fundamental step that consists of rotating the transverse components of the six spins around a given plaquette by $\pi/4$ in alternating directions. Doing so changes $B$ on the chosen plaquette by $-3\pi/2$ and by $\pi/2$ on its neighbours. Assuming that $B$ was small initially, one can regard the former as a change of $B$ by $\pi/2$ together with a $2\pi$ phase slip that inserts a vison pair across the plaquette in question. That is, this elementary step inserts a pair of visons with a symmetric near-field pattern around them. Repeating it along a chain of sites amounts to moving the visons apart from one another, similarly to the insertion of spinons in classical spin ice.

The state generated by these local operations, however, is not the least energetic one, for they do not capture the change in $B$-field away from the visons or the “Dirac string” connecting them. While it may be possible to construct operators acting on all spins that have a larger overlap with the “true vison creation operator” Hermele et al. (2004), we adopt here a simpler and more straightforwardly reliable approach. We equilibrate the photons generated by the local vison creation operation using the Monte Carlo algorithm and gradually reduce the effective temperature until the remaining photon population can be ignored, leaving behind a two-vison metastable state. The only issue of this method is vison movement: Visons in the semiclassical model are not inherently mobile Kwasigroch (there are no explicit vison hopping terms in the Hamiltonian), but a high-temperature photon cloud can move them around. We find that starting the annealing procedure at a sufficiently low temperature prevents such motion unless the visons are introduced within a distance of about $2a_{0}$ (where $a_{0}$ is the fcc lattice parameter).

The energy difference between the resulting configuration and the ground state can be regarded as the energy cost of two visons plus their Ewald summed interaction energy. The latter, however, contains a surface term de Leeuw et al. (1980) which makes interpreting the results complicated, see Appendix B.2. This can be remedied by using an arrangement of visons with no net dipole moment, such as the rock salt configuration in Fig. 3(a). This arrangement was set up using the local vison insertion protocol described above; photons were equilibrated for $16L$ Monte Carlo steps 515151A “Monte Carlo step” in this paper consists of the following: sampling the $xy$ phase angle $\phi$ of each spin; a Metropolis attempt to change $\sigma^{z}$ around each hexagonal plaquette; and sampling the gauge freedom of $\phi$ on each tetrahedron. These elementary steps are described in more detail in Sec. III. at $T=g/256$; then cooled 32 times by a factor of 2 and equilibrated for $8L$ steps each time. In the end, the temperature of the photon cloud reached $2^{-40}g\approx 9\cdot 10^{-13}g$; equilibration was monitored through the acceptance rate of Metropolis steps which remained consistently high throughout the process. The so obtained metastable energies are plotted in Fig. 3(a).

Modelling visons as Coulomb interacting charges, the energy of the rock salt configuration is

[TABLE]

per vison, where $\mu$ is the energy cost of an isolated vison, $L$ is the size of the cubic simulation box in units of the fcc lattice parameter $a_{0}$, $\alpha$ is the Coulomb interaction strength, and $M=1.74756459\dots$ is the appropriate Madelung constant. Fitting this form to the numerical results plotted in Fig. 3(a) yields

[TABLE]

The numerical results can be compared to an analytic estimate of the energy cost and interaction strength of visons from a quadratic approximation to the energy of the magnetic field, $-g\cos B$, see Appendix B.1. The quadratic estimate of $\mu$, $8.858g$, is substantially different from the numerical result (13); however, the interaction strength, $\alpha=\pi ga_{0}$, matches excellently. This is so because the vison energy cost includes that of its immediate neighbourhood, where the quadratic theory breaks down; on the contrary, most of the interaction energy is due to the overlap of far fields for which the quadratic theory is accurate.

We finally consider how accurately the interaction of nearby visons is described by the asymptotic Coulomb law. Since these visons were unstable against the photon cooling protocol, an alternative technique had to be used. Within quadratic approximation, the $B$-field configuration that minimises the magnetic energy for a given arrangement of visons can be found explicitly, as described in Appendix B.1: The sum of $-g\cos B$ for all plaquettes in this configuration is an upper bound on the true energy of the vison configuration. We benchmarked this estimate against the photon cooling technique at vison separations where the latter is applicable; because of the good agreement there, we decided to use this method to estimate the energy of vison dipoles with separations smaller than $2a_{0}$. The results are plotted in Fig. 3(b) together with the Coulomb energy estimate $2\mu-\alpha/r$, where $\mu$ and $\alpha$ are taken from (13). The latter is a remarkably good approximation even at third nearest neighbour distance; there is a discrepancy of about $0.5g$ for next-nearest neighbours and about $2g$ for nearest neighbours.

Most notably, the energy of a nearest neighbour dipole is smaller than that of a single isolated vison, $\mu$. Visons in QSI thus form a weak electrolyte. While deconfined, their dissociation is so energetically unfavourable that most visons at low temperatures remain associated with an oppositely charged one. This behaviour is quite different from that of spinons in CSI, whose energy cost is set independently by the dominant Ising exchange interaction, and they remain energetically dissociated even at nearest-neighbour distance in the presence of entropic and dipolar interactions.

VI Thermodynamics of photons and visons

A natural way to introduce gapped excitations in numerical simulations is via thermal fluctuations at finite temperature where an equilibrium population of such excitations arises. Thermodynamic quantities like heat capacity or thermal conductivity are promising signatures of gapped quasiparticles and QSL behaviour in general Knolle and Moessner (2019); Rau and Gingras (2019); Tokiwa et al. (2018). Therefore, we studied the interactions of visons and photonic modes in a finite temperature ensemble. This was greatly aided by the ability of our method to directly access observables such as the emergent magnetic field and the vison charge operator (11).

VI.1 Magnetic field pinch points

We evaluated static correlation functions of the emergent magnetic field, $\langle B(-\mathbf{q})B(\mathbf{q})\rangle$, at temperatures between $0.4g$ and $1g$ using Monte Carlo spin configurations of a cubic sample of size $L=20$ unit cells 61616155 temperature points were used, uniformly distributed in $1/T$. For temperatures above $0.5g$, 131 072 stochastically independent Monte Carlo samples were generated; for those between $0.4g$ and $0.5g$, 262 144 samples were used.. These correlators are plotted in the $(hhk)$ plane for three temperatures in Fig. 4; the correlators of the emergent electric field $\sigma^{z}$ are also shown for comparison. At low temperatures, both correlators exhibit sharp pinch points at the $\Gamma$ points in the pattern familiar from classical spin ice Isakov et al. (2004); Henley (2005). The pinch points remain sharp for the electric field at all temperatures as spinons are excluded by the Monte Carlo algorithm. On the contrary, the introduction of visons blurs the $B$-field pinch points in much the same way monopoles blur CSI pinch points Fennell et al. (2009); Sen et al. (2013). This picture is in some departure from spin-$1/2$ QSI in which pinch points are suppressed near the $\Gamma$ points and thus no sharp features appear Shannon et al. (2012); Benton et al. (2012). This is due to the different way in which photons appear in the two systems. In the semiclassical case, they are classical normal modes, hence their energy content is constant as per equipartition Note (49); in the spin-$1/2$ case, low-temperature physics is dominated by quantum zero-point fluctuations of photon modes which give rise to correlators proportional to the photon dispersion $\omega(\mathbf{q})$, thus suppressing the pinch points.

Pinch point blurring is a common experimental diagnostic of spinons in CSI Fennell et al. (2009); likewise, we were able to extract quantitative information about the visons from the blurring of $B$-field pinch points. We focus on the $(00k)$ axis (green line in Fig. 4), where the photon contribution to the correlator vanishes Isakov et al. (2004); Henley (2005), that is, the signal is entirely due to visons. Cuts of the correlator along this axis are plotted in Fig. 5(a) for four temperatures. These show an apparently Lorentzian peak at the $(002)$ pinch point, indicating a Debye plasma of visons. Unlike pinch point blurring in CSI, however, our blurring pattern is not explained by a Lorentzian peak in itself; in particular, the peak appears on top of a substantial constant background, as demonstrated in Fig. 5(b). This constant correlator can be ascribed to a large population of nearest neighbour vison dipoles that are independent of the Debye plasma mentioned above; for a detailed derivation, see Appendix C. This underlines the observation that visons in QSI form a weak electrolyte, that is, they interact strongly enough that a large fraction of their thermal population remains associated, as discussed in Sec. V.

Furthermore, the Coulomb energy formula for vison pairs substantially overestimates the energy cost of next-nearest-neighbour vison dipoles, see Fig. 3(b). Therefore, we anticipate an excess population of them compared to that predicted by the Debye–Hückel approximation. While this effect is not qualitative, it does introduce a correction to the $B$-field correlator that is proportional to $\sin^{2}(q_{z}/8)$ along the $(00k)$ axis, see Appendix C.2. To take into account the effects of these closely associated dipoles, we fitted the functional form

[TABLE]

to the data at all temperature points; $\xi\propto\rho_{\mathrm{free}}^{-1/2}$ is the Debye length of the plasma formed by dissociated dipoles, where $\rho_{\mathrm{free}}$ is the density of dissociated visons; while $C_{1}$ is proportional to the density of nearest neighbour vison pairs, $\rho_{\mathrm{nn}}$ 626262The complicated form of the Lorentzian is to account for the periodicity of the data imposed by the lattice; $C_{2}$ is related to the excess density of next-nearest-neighbour dipoles, a quantity that is hard to get a direct handle on. Equation (14) fits the data excellently throughout the investigated temperature range; Fig. 5(b) demonstrates that all three terms are necessary to achieve this, although a very good agreement is already obtained without the $C_{2}$ contribution.

VI.2 Temperature dependence of vison density

Visons can also be counted explicitly in our Monte Carlo simulations. Besides the number $N$ of all visons, the number $N_{\mathrm{nn}}$ of dual diamond lattice bonds with two visons of opposite charge on their ends was obtained; from these, the number of dissociated visons was estimated as $N_{\mathrm{free}}=N-2N_{\mathrm{nn}}$. The density of bound dipoles and free visons was plotted as a function of inverse temperature in Fig. 6, together with the constant background $C_{1}$ of the $B$-field correlator and the Debye length $\xi$, respectively; we expect $\xi\propto\rho_{\mathrm{free}}^{-1/2}$ Levin (2002) and $C_{1}\propto\rho_{\mathrm{nn}}$ (see Appendix C.2). These relations hold quite well throughout the temperature range, confirming that magnetic field pinch point blurring is a good measure of vison populations.

The densities of both dissociated and bound visons follow an approximate Arrhenius law at low temperatures with an apparent gap close to the bare vison energy (13) and the nearest neighbour dipole energy shown in Fig. 3(b), respectively. However, these vison densities saturate at $T\approx g$, an order of magnitude below the bare energy costs, and correspondingly, densities at $T\lesssim g$ are much larger than predicted by a simple Boltzmann factor with quasi-equilibrium gaps (red lines in Fig. 6).

VI.3 Thermodynamic effects of quasiparticle interactions

In order to identify the origin of this excess density of visons, we considered two different estimates of their energy cost in the thermal ensemble: the dependence of energy on vison number, quantified by $\mathrm{d}\langle E\rangle/\mathrm{d}N$, and the slope of the Arrhenius plot $\log N$ versus $1/T$. Both quantities were evaluated using data recorded at single temperature points, similarly to heat capacity estimation using the fluctuation–dissipation theorem Landau and Lifshitz (1980):

[TABLE]

for a derivation, see Appendix D. These estimates are plotted in Fig. 7, together with $\mu_{\mathrm{Arrh.}}/\mu_{E}=\operatorname{var}N/N$. At low temperatures, the distribution of $N$ is generated by a Poisson distribution of thermal excitations which may well be a collective one made up of several visons: It follows (see Appendix D) that $\operatorname{var}N/N$ gives the typical vison content of such a collective excitation.

At the lowest temperatures, $\operatorname{var}N/N$ tends to 2; this again demonstrates that the dominant vison species at low temperatures are bound dipoles. At intermediate temperatures, $\operatorname{var}N/N$ increases further, to about 2.5, indicating collective excitations consisting of more than two visons. Salient examples of such clusters could be ‘polarons’ consisting of a free vison and nearest-neighbour dipoles. As the $B$-field emanating from the former is quite large in its immediate neighbourhood, the energy cost of appropriately aligned dipoles decreases substantially, causing them to proliferate. An effective phenomenological model based on this mechanism that allows us to gain a qualitative understanding is presented in Appendix E. At even higher temperatures, $\operatorname{var}N/N$ decreases as $N$ approaches saturation.

We also note that the energy cost of visons in the thermal ensemble declines steeply for $T\gtrsim 0.7g$. This is due to interactions between visons and photons arising from corrections to the gauge theory beyond quadratic order. The next order in the Villain expansion Kwasigroch et al. (2017) consists of negative quartic terms which, in the presence of a (thermal) population of photons, lead to a negative renormalisation of vison energies and interactions; see Appendix E for a more detailed discussion of both the mechanism and its effect on the effective vison gap. This mechanism is at work at intermediate temperatures, reducing the gap of visons compared to the zero-temperature value, resulting in larger thermal populations. This may also explain the excess vison density observed in Fig. 6.

As temperature is increased, the energy cost of visons drops substantially and their density saturates. Eventually, at temperatures $T\gg g$, the visons cease to be useful quasiparticles as it is to be expected from the perspective of the individual spins: The contribution of each hexagonal plaquette to the Hamiltonian (5) is $O(g)$, therefore, all spin configurations satisfying the no-spinon constraint (which is still enforced by the Monte Carlo algorithm) become roughly equally likely, regardless of vison content. The situation is similar to the crossover of classical spin ice into a high-temperature paramagnetic phase as spinons cease to be useful quasiparticles to describe the system. In our case, visons and photons are washed out at high temperatures, giving rise to a classical spin ice phase with the spinon as its only excitation. Indeed, the heat capacity of the system, plotted in Fig. 7, displays a Schottky peak at $T\approx 0.7g$ above which it drops below the equipartition heat capacity of photon modes, $k_{\mathrm{B}}/2$ per spin, and tends to zero as $T\to\infty$. This indicates that photons break down as quasiparticles together with visons.

VII Conclusion

We developed a semiclassical numerical technique to simulate quantum spin liquids whose physics is underpinned by perturbative ring exchange processes. Thus far, these systems remained elusive to classical simulations, as their native bilinear exchange Hamiltonians tend to predict ordered or Ising spin liquid phases, devoid of QSL behaviour. In contrast, we studied effective ring-exchange Hamiltonians directly, which are expected to remain in the quantum spin liquid phase in the large-$S$ limit. We formulated our method in terms of a large-$S$ path integral formalism and demonstrated that a combination of classical Monte Carlo sampling and Landau–Lifshitz dynamical evolution captures the least action trajectories of the finite-temperature path integral of the system. Treating spins semiclassically allowed us to access a variety of observables not readily available to other simulation methods. We thus gained a detailed insight into the behaviour of quasiparticles in these quantum spin liquids, including gapped excitations which are not yet amenable to other computational and analytic techniques.

We demonstrated the potential of our method on pyrochlore quantum spin ice, a paradigmatic $U(1)$ QSL Hermele et al. (2004). We identified a gapless, linearly dispersing branch of classical normal modes in spin dynamics. At low temperatures, these photonic modes are remarkably sharp and their dispersion matches analytic results from large-$S$ path integral calculations extremely well Kwasigroch et al. (2017). We also observed gapped, quantised vortices of the emergent magnetic field with the phenomenology expected for the elusive $U(1)$ vison quasiparticles of QSI. In a showcase of the capabilities of our method, we were able to introduce these visons in a controlled way in the system and study their interactions in vacuum (that is, in an effectively zero-temperature photon background). In contrast with the behaviour of spinons, the energy cost and interaction strength of visons are set by a single parameter in the original Hamiltonian, that is, they are not separated parametrically. We found that the interaction of a nearest-neighbour pair of visons is strong enough that the energy cost of the resulting dipole is less than that of an isolated vison; namely, the visons form a weak electrolyte. At low temperatures, this leads to a relatively large population of closely associated visons which must be taken into account in modelling and understanding their behaviour, in particular in effective descriptions of vison dynamics Chen (2016, 2017); Kwasigroch .

We also studied the interplay of photonic modes and visons in thermal equilibrium at finite temperature. We developed a detailed understanding of visons through directly tracking them and studying their effect on pinch point patterns of the relevant static correlation functions. Similarly to the well understood case of spinons in classical spin ice, these pinch points blur due to the presence of visons; by contrast, however, they are also uniformly lifted by a population of associated visons. These effects can be used to reliably measure the density of associated and dissociated visons separately. We find that they are thermally activated with effective gaps comparable to their respective energy costs in the zero-temperature limit.

A simple Boltzmann factor calculation, however, substantially underestimates vison populations at intermediate temperatures, which also saturate at temperatures far smaller than the $T\to 0$ vison cost. We understand this to be a consequence of strong interactions between visons and a highly excited photon background. The latter introduces fluctuations in the emergent electric and magnetic fields that blur visons, reducing their energy and leading to larger populations than naively expected. In spin-1/2 QSI, vacuum fluctuations of the photon modes occur even at zero temperature; we anticipate that these would similarly reduce the chemical potential of visons, possibly bringing it within the photon band. Since vison–photon interactions are controlled by the same ring-exchange energy scale as the energy of both quasiparticles, visons and photons of similar energy may then readily hybridise, leading to composite particles with more elaborate and exciting features than those of either visons or photons. It would be interesting to further test and confirm this behaviour in a spin-1/2 quantum setting, e.g., using quantum Monte Carlo. In particular, our work suggests that static correlators of the emergent magnetic field, possibly cumbersome but likely accessible to QMC, present a promising angle to study visons. Such studies would also be instrumental in confirming the weak electrolyte behaviour found in the semiclassical limit.

One direct experimental signature that was identified in our work is the Schottky peak in the specific heat due to visons (see Fig. 7). However, one must keep in mind that quantum photon excitations also contribute a Debye term to the specific heat, which is absent in semiclassical photonic normal modes. Hybridisation between photons and visons may well merge the two contributions and alter the shape of the Schottky anomaly.

The semiclassical numerical technique can naturally be applied to the dynamics of a high-temperature ensemble that contains thermally generated visons. However, the dynamics of visons is far from trivial, since they are immobile at zero temperature Kwasigroch and their motion at finite temperature is due to being “tugged” by the photon background. Further work is needed to gain better insight into this behaviour.

A particularly interesting experimentally relevant direction is understanding magnetic noise in QSI Dusad et al. (2019): Here, our method could be readily used to investigate the low temperature regime where the leading magnetic noise contribution is due to photons and spinon excitations can be neglected. Similarly, our approach could be applied to better understand the photon and vison contribution – and their interplay, in particular the weak electrolyte behaviour – to thermal transport in candidate QSI materials Tokiwa et al. (2018).

Whereas the gauge charge and fields due to (static) visons are not directly accessible in experiments, it is often the case that gauge degrees of freedom also carry irrational physical ones Moessner and Sondhi (2010). For example, magnetostriction or Dzyaloshinskii–Moriya effects can associate real electric polarisation with the emergent gauge fields of QSI Khomskii (2012); Lantagne-Hurtubise et al. (2017); Nakosai and Onoda (2019). As a result, one may expect, for instance, the vison pairs to have a signature in the dielectric response and electric susceptibility of QSI materials.

Our method may be extended to include terms in the simulated Hamiltonian that enable introducing static or dynamical spinons. Since spinons appear to have salient experimental signatures in, for instance, magnetisation and neutron scattering measurements Huang et al. (2018), a better understanding of spinon–vison interactions through our simulations may provide experimentally accessible handles to study collective photon and vison behaviour. It is important to note, however, that spinons are not quantised in the semiclassical setting (see Sec. II); therefore, understanding what their behaviour tells about the original quantum problem requires some care.

Beyond QSI, our approach is manifestly applicable to other quantum spin liquids underpinned by ring exchange processes. Quantum dimer models are a salient example, especially in light of novel large-$S$ analytical approaches Szabó et al. (2019) that apply also to non-bipartite lattices and their $\mathbb{Z}_{2}$ RVB phases Moessner and Sondhi (2001). Our technique could add to the understanding of these systems, in particular to the interplay of the $\mathbb{Z}_{2}$ vison excitations with the dimer liquid background as correlations develop, e.g., upon approaching quantum critical points out of the $\mathbb{Z}_{2}$ RVB phase.

Finally, the ability of the method to naturally include quantised charges [namely, U(1) visons] in an effectively classical system may have ramifications for lattice gauge theories in general. Vison excitations in our semiclassical model are quantised solitons, similarly to Dirac monopoles in quantum electrodynamics, and like these monopoles, they are likely to be in a strong coupling regime. The reduction of vison energy by a thermally fluctuating background may indeed be a semiclassical analogue of running couplings that are brought about by virtual particle–antiparticle bubbles in QED. This is but a speculative yet intriguing potential connection between QSI and the strong coupling problem in QED, which warrants further investigation in future work.

Acknowledgements

We thank Cristian Batista, John Chalker, Felix Flicker, Roderich Moessner, Nic Shannon, and David Tong for stimulating discussions. We are particularly grateful to David Ho for bringing the idea of a possible analogy to the strong coupling problem in QED monopole production to our attention, and to Michał Kwasigroch for sharing preliminary results of his work Kwasigroch and for engaging in fruitful conversations with us on numerous occasions. Fig. 1 was prepared using VESTA 3 Momma and Izumi (2011). Figs. 2, 4, and 9 use perceptually uniform color maps developed in Ref. Kovesi (2015). This work was supported in part by Engineering and Physical Sciences Research Council (EPSRC) Grants No. EP/P034616/1, EP/M007065/1, and EP/K028960/1.

Appendix A

Scalar and vector fields on the diamond lattice

Lattice vector fields can be defined on the bonds of any bipartite lattice; in order to imbue them with a sense of directionality, their sign is taken to depend on the orientation $\mathbf{r}\to\mathbf{r}^{\prime}$ of each bond: $A_{\mathbf{r}\to\mathbf{r}^{\prime}}=-A_{\mathbf{r}^{\prime}\to\mathbf{r}}$. Since the lattice is bipartite, all bonds can be oriented consistently from one sublattice towards the other, giving rise to a natural and consistent sign convention.

Scalar fields can also be introduced. They live on the sites of the lattice and, like scalar fields in $\mathbb{R}^{d}$, their sign is uniquely defined. The lattice gradient of a scalar field and the lattice divergence of a vector field can then be written as

[TABLE]

It is easy to see that the gradient is a well-defined vector field while the divergence is a scalar field.

The spins in quantum spin ice live on a pyrochlore lattice formed by the bond midpoints of a bipartite diamond lattice (cf. Fig. 1). By employing the sign convention mentioned above, $\sigma^{z}$ and the $xy$ phase angle of the spins can be turned into a lattice electric field and vector potential, respectively Hermele et al. (2004). A peculiarity of the pyrochlore lattice is that the centres of its hexagonal plaquettes form a dual pyrochlore lattice, whose dual is in turn the original lattice. This allows for defining the lattice curl as the sum of a vector field around a hexagonal plaquette, traversed according to the right hand rule with respect to the direction of the dual diamond lattice bond corresponding to that plaquette. Upon flipping this direction, the sense of circulation is changed, which in turn flips the sign of all the terms involved. Therefore, the lattice curl is a well-defined vector field of the dual lattice. Furthermore, the diamond sites around each plaquette alternate between the two sublattices, and the lattice curl includes the field variables with alternating signs as compared to the unique sign convention. Ultimately, this leads to the familiar vector calculus identities $\operatorname{div}\operatorname{curl}A=0$ and $\operatorname{curl}\operatorname{grad}U=0$.

Appendix B Quadratic estimates of the vison cost and interaction

Following Ref. Kwasigroch et al. (2017), we estimate the zero-temperature gap and interaction strength of two visons in a quadratic theory where the energy of the magnetic field on each plaquette, $-g\cos B$, is approximated as $-g+gB^{2}/2$. Let the number of visons on each site of the dual diamond lattice be $n(\mathbf{r})$. We want to find the lowest energy configuration of the magnetic field for which $\operatorname{div}B=2\pi n(\mathbf{r})$. To do so, we consider a lattice version of Helmholtz’s theorem: Any vector field $B(\mathbf{r})$ on the dual pyrochlore lattice with periodic boundary conditions can be written as

[TABLE]

where $U$ is a scalar potential defined on the dual diamond lattice, $A$ is a vector potential defined on the original pyrochlore lattice (both with the same periodic boundary conditions as $B$), and $C^{\mu}$ is a constant that can vary amongst the four fcc sublattices $\mu=0,1,2,3$ of the dual pyrochlore lattice. Now, the cross-terms between the three components vanish in $\sum B^{2}$, that is, the approximate total energy of a field configuration can be written as

[TABLE]

Since $U$ is uniquely determined by $n$, $E$ is minimised if $\operatorname{curl}A$ is identically zero, that is, there are no non-gradient components of $B$.

B.1 Coulomb contribution

We focus initially on the first two terms in Eq. (18) and we shall return to $C^{\mu}$ later in Sec. B.2. The scalar potential $U$ can be obtained by requiring that

[TABLE]

The divergence and the gradient can be expressed in reciprocal space as

[TABLE]

where the indices $\alpha=\pm$ and $\mu=0,1,2,3$ indicate the fcc sublattices of the diamond and pyrochlore lattices, respectively, and

[TABLE]

where $\mathbf{r}_{0}=a_{0}[111]/4$, $\mathbf{r}_{1}=a_{0}[1\overline{11}]/4$, $\mathbf{r}_{2}=a_{0}[\overline{1}1\overline{1}]/4$, and $\mathbf{r}_{3}=a_{0}[\overline{11}1]/4$ are the vectors pointing from a “+” diamond lattice site to its nearest neighbours. It then follows that (20) can be written as (we suppress the arguments $\mathbf{q}$ from now on)

[TABLE]

where we introduced the vectors $\mathbf{n}=[n^{+},n^{-}]$, $\mathbf{U}=[U^{+},U^{-}]$, and $\mathbf{B}=[B^{0},B^{1},B^{2},B^{3}]$ for convenience. Equation (24) gives all $\mathbf{q}\neq 0$ Fourier components of the ground state magnetic field configuration. [$\mathbf{M}\mathbf{M}^{\dagger}$ is singular at the $\Gamma$ point: this is the origin of the constant term in (18).] Now, the total energy of this configuration in the quadratic approximation is

[TABLE]

where $N$ is the number of fcc unit cells of volume $V_{\mathrm{cell}}=a_{0}^{3}/4$. We also note that

[TABLE]

where $\gamma=\frac{1}{4}\sum e^{i\mathbf{q}\cdot\mathbf{r}_{\mu}}=\cos(q_{x}/4)\cos(q_{y}/4)\cos(q_{z}/4)-i\sin(q_{x}/4)\sin(q_{y}/4)\sin(q_{z}/4)$ Henley (2005).

To obtain the interaction energy of two visons from the general form (25), consider a positive and a negative vison, both on the “+” sublattice, a distance $\mathbf{R}$ away from each other. That is, let $n(0)=+1$, $n(\mathbf{R})=-1$, and $n(\mathbf{r})=0$ otherwise. In reciprocal space, this yields $n^{+}(\mathbf{q})=1-e^{i\mathbf{q}\cdot\mathbf{R}}$ and $n^{-}(\mathbf{q})=0$. Substituting this into (25) gives

[TABLE]

The two terms of (32) can be regarded as the energies of the two isolated visons and their interaction, respectively. The bare vison energy works out to be $\mu\approx 8.848g$. For large $\mathbf{R}$, the second term integrates to zero over most of the Brillouin zone due to the rapidly oscillating cosine factor; the only exception is the vicinity of the $\Gamma$ point where $\gamma=1$, and so the integrand diverges. There, $|\gamma|^{2}\approx 1-\mathbf{q}^{2}a_{0}^{2}/16$, and so the integral becomes

[TABLE]

which demonstrates the effective Coulomb interaction between visons at large distances, arising from a purely short-ranged Hamiltonian.

B.2 Ewald surface term

In the discussion above, we did not consider the $\mathbf{q}=0$ term of the Fourier transformed energy (25) because it is not uniquely determined by the $\operatorname{div}B=2\pi n$ condition. This ambiguity is manifest in the constants $C^{\mu}$ in the Helmholtz decomposition (18). To derive these constants, we have to note that the origin of visons is a $2\pi$ phase ambiguity in the emergent vector potential $\phi$, the curl of which is the magnetic field $B$. We can fix this ambiguity by giving each $\phi$ a unique value and by introducing an integer-valued field $\mathcal{B}$ that ensures $|B|<\pi$ for each plaquette Hermele et al. (2004):

[TABLE]

As $\phi$ is now uniquely defined, $\operatorname{div}\operatorname{curl}\phi=0$, and so

[TABLE]

This means that visons can only be introduced by changing $\mathcal{B}$ on some plaquettes, that is, by phase slips.

It is easy to see that the lattice curl of any field with wave vector $\mathbf{q}=0$, that is, one that is constant on each pyrochlore sublattice, is uniformly zero. As a result, the $\mathbf{q}=0$ component of the magnetic field (34) is entirely due to the $2\pi\mathcal{B}$ term, and we now show that it depends only on the position of the visons. To do so, consider introducing a pair of visons with (positive to negative) separation $\Delta\mathbf{r}$. This involves adding $\pm 1$ to $\mathcal{B}$ on plaquettes along some path between the two end points: Heading from the positive to the negative vison, a plaquette encountered going from a “$+$” to a “$-$” tetrahedron gets a $+1$, while a “$-$” to a “$+$” one gets a $-1$. Since each of these steps is associated with a vison movement $\pm\mathbf{r}_{\mu}$, respectively, we have

[TABLE]

where $\mathcal{B}^{\mu}$ is the total $\mathcal{B}$ on pyrochlore sublattice $\mu$ (that is, its $\mathbf{q}=0$ Fourier component). Equation (36) gives three equations for the four $\mathcal{B}^{\mu}$. A fourth one can be obtained by realising that the total $\mathcal{B}$ changes by $+1$ when a vison moves onto the “$+$” diamond sublattice and by $-1$ when a vison moves away from it, that is, the total $\Delta\mathcal{B}$ upon inserting a vison pair equals the total change in $q_{+}$, the net vison charge in the “$+$” sublattice:

[TABLE]

Since all visons are located at the ends of Dirac strings, one could introduce them from a setup containing no visons by repeating this operation. Summing (36, 37) for all of these vison creation events gives

[TABLE]

where $\mathbf{P}=-\sum\Delta\mathbf{r}$ is the total dipole moment of the vison configuration. From (34, 38), the contribution of the $\mathbf{q}=0$ component to the total energy can be written as

[TABLE]

where $V$ is the volume of the system with periodic boundary conditions and we used that $\mathbf{N}/2$ is orthonormal.

This result can be connected to the standard surface term in Ewald summation de Leeuw et al. (1980), which is of the form

[TABLE]

where $\kappa$ is the Coulomb constant defined by $V_{ij}=\kappa q_{i}q_{j}/(4\pi r)$ and $\varepsilon$ is the relative permittivity of the medium surrounding the system. In our case, $\kappa=4\pi^{2}ga_{0}$ and $\varepsilon=0$, for the emergent magnetic field does not even exist outside of quantum spin ice. Substituting this into (40) reproduces the $\mathbf{P}$-dependent part of (39).

B.3 Comparison to numerical simulations

In order to confirm the validity of the quadratic theory for interactions, we used the photon cooling protocol described in Sec. V to measure the energy of two visons separated by a range of distances in a cubic simulation box of size $32a_{0}$ in the zero-temperature limit. The results are plotted in Fig. 8 together with the predictions from Ewald summation of the Coulomb interaction (33) with surface term (39), and the bare energy cost of visons taken from (13). This prediction has no free parameters and agrees excellently with the numerical simulations.

Appendix C Pinch point blurring due to free and bound visons

We derive the contribution of visons to the equal time correlator of the emergent magnetic field, $\langle B(\mathbf{q})B(-\mathbf{q})\rangle$, for several arrangements of visons that are realised in the semiclassical quantum spin ice model. We assume throughout that there is no interaction between photons and visons, that is, the overall correlator can be written as the sum of independent vison and photon contributions.

Let the vison number on each site of the dual diamond lattice be $n(\mathbf{r})$, which can be Fourier transformed into $n^{\alpha}(\mathbf{q})$, as done in Appendix B. The contribution of these visons to the magnetic field is given by (24), and therefore

[TABLE]

where $\mathbf{1}$ is the row 4-vector all entries of which are 1, and we introduced the row 2-vector $\mathbf{v}=2\pi\mathbf{1}\mathbf{M}^{\dagger}(\mathbf{M}\mathbf{M}^{\dagger})^{-1}$. The latter can be written as

[TABLE]

where $\kappa=\frac{1}{4}\sum_{\mu}e^{i\mathbf{q}\cdot\mathbf{r}_{\mu}/2}$. Therefore, the vison contribution to the correlator of $B(\mathbf{q})$ follows from that of $n(\mathbf{q})$ as

[TABLE]

Below we derive the correlator of vison number, and hence its contribution to that of the magnetic field, for vison arrangements relevant to our model. The most salient feature of the full $\langle B(-\mathbf{q})B(\mathbf{q})\rangle$ correlator are the pinch points. For convenience, we focus on the behaviour along the $q_{z}$ axis, where the photon contribution vanishes Henley (2005), and so the $\langle BB\rangle$ correlator measured along it is entirely due to visons and allows for direct comparison with the result in Eq. (43).

C.1 Debye plasma of dissociated visons

Let the visons interact through the reduced Coulomb interaction $\beta V(r)=Kq_{i}q_{j}/(4\pi r)$. Assume that the density $\rho$ of visons is small, that is, both their typical separation and the Debye screening length of the resulting plasma is much larger than the lattice spacing. Accordingly, we focus on the vicinity of $\Gamma$ points in reciprocal space. From standard Debye–Hückel theory Levin (2002), the pair correlation function of visons is

[TABLE]

where $\xi_{\mathrm{D}}=(K\rho)^{-1/2}$ is the Debye screening length and $a,b=\pm$ denote the positive and negative vison species. To go from (44) to the correlators $\langle n^{\alpha}n^{\beta}\rangle$, we note that the long-wavelength theory does not discriminate between the two sublattices. Therefore, the density of visons in each is $\rho/2$ and $\langle n^{\alpha}n^{\beta}\rangle$ does not depend on $\alpha$ and $\beta$ in the long wavelength limit. In this approximation,

[TABLE]

The above derivation captures the correlation between pairs of two visons; however, the correlation of visons with themselves also contributes to the $\langle nn\rangle$ correlator. This contribution is clearly a $\delta$-function in real space, and so only couples each sublattice to itself; after Fourier transforming, we get

[TABLE]

We now calculate the magnetic field correlator $\langle B(-\mathbf{q})B(\mathbf{q})\rangle$ using (43) along the $q_{z}$ axis near the $\mathbf{G}=(002)$ pinch point. Since $\mathbf{G}$ is a $\Gamma$ point of the fcc lattice, both sublattices of the diamond lattice behave the same way as at $\mathbf{q}=0$, but a relative phase $\mathbf{G}\cdot\mathbf{r}_{\mu}=\pi$ is introduced between them. In practice, this means that the $\langle n^{\pm}n^{\mp}\rangle$ correlators are the opposite of what they were near $\mathbf{q}=0$, while the sign of $\langle n^{\pm}n^{\pm}\rangle$ remains unaffected. Equations (43, 46, 47) now yield

[TABLE]

where $\mathbf{q}=(0,0,4\pi/a_{0}+k)$. That is, the Debye plasma of visons introduces a Lorentzian blurring of the pinch points. The width of this blurring is $\xi^{-1}_{\mathrm{D}}\propto\sqrt{\rho}$.

C.2 Tightly bound dipoles

Since the dominant vison species at low temperatures is not the isolated vison, but a dipole of nearest neighbour visons, we need to derive $\langle B(-\mathbf{q})B(\mathbf{q})\rangle$ for a gas of nearest neighbour dipoles. Furthermore, the energy of a second neighbour dipole is significantly smaller (by about $0.48g$) than predicted by the simple Coulomb approximation, see Fig. 3(b). For these reasons, the Debye plasma approximation will significantly underestimate the population of such dipoles, which must be corrected for explicitly. The only effect of these associated dipoles is correction to the polarisability of the emergent magnetic field. Therefore, their interactions with each other can be neglected, and the only contribution to $\langle nn\rangle$ correlators comes from visons within the same dipole.

For nearest neighbour ($\mathbf{r}_{1}-\mathbf{r}_{2}=a_{0}\langle 111\rangle/4$) dipoles, the real space vison correlators are thus

[TABLE]

where $\rho$ is the density of dipoles; the factor of $1/4$ is due to the four possible orientations of the dipole. Clearly, the two visons are on the same sublattice in the first line in (49) and on different ones in the second. Therefore, the reciprocal space correlators are

[TABLE]

from which the magnetic field correlator follows as

[TABLE]

which is constant along the $q_{z}$ axis. That is, the contribution of nearest neighbour dipoles to the magnetic field correlator is a uniform background that gradually submerges the pinch point, eventually washing it out.

For second neighbour dipoles ($\mathbf{r}_{1}-\mathbf{r}_{2}=a_{0}\langle 110\rangle/2$), the real space vison correlators are

[TABLE]

Since the two components of the dipole are on the same sublattice, the reciprocal space correlators are

[TABLE]

which is proportional to $\sin^{2}(q_{z}a_{0}/8)$ along the $q_{z}$ axis. This explains the small but significant cosine modulation of the lifting of the pinch point.

Besides the $q_{z}$ axis, we considered the behaviour of the $\langle BB\rangle$ correlators around the pinch points, especially the one at $\mathbf{G}=(002)$. Vison number correlators are analytic in all cases; however, the “response function” $\mathbf{v}$ gives rise to pinch points perpendicular to those due to photons (see Fig. 4). The Debye plasma contribution is peaked at the pinch point. Its maximum matches the intensity of the photon contribution there, resulting in Lorentzian pinch point blurring Henley (2005); Sen et al. (2013). On the other hand, closely associated dipoles result in contributions proportional to the density throughout reciprocal space, which reduce the overall contrast of the pinch points. For completeness, the $\langle BB\rangle$ correlators produced by first and second neighbour dipoles, (51) and (54), are plotted on the $(hhk)$ plane in Fig. 9.

Appendix D Estimates of the vison gap from thermal statistics

We discuss two methods to estimate the chemical potential of visons in a thermodynamic set-up. We also derive an expression for both of them in terms of statistics of the vison number and the energy at a single temperature, in a similar fashion to the fluctuation–dissipation theorem. We treat our model as a classical thermodynamic system with microstates of well-defined energy $E_{\alpha}$ and vison number $N_{\alpha}$, weighted according to the partition function

[TABLE]

where $\zeta$ is a fictitious chemical potential introduced to keep track of vison number; in the physical partition function, $\zeta=0$.

The first approach estimates the energy cost of a single vison directly, which can formally be written as $\mu_{E}=\mathrm{d}E(N)/\mathrm{d}N$, where $E(N)$ is the mean energy of the system constrained to $N$ visons. In the thermodynamic limit, this derivative is equivalent to the ratio of variations in $\langle E\rangle$ and $\langle N\rangle$ due to a changing chemical potential:

[TABLE]

Each of the two derivatives in (56) can be expressed in terms of derivatives of $Z$, which can in turn be rewritten in terms of statistics of energies:

[TABLE]

and therefore

[TABLE]

Another estimate of the excitation gap is the local slope of the Arrhenius plot $\log N$ vs. $1/T$, at least at low temperatures where $N$ has not saturated. However, this is not necessarily the gap of a single vison, but of whatever (possibly multi-vison) excitations are created thermally in the system. This slope is given by

[TABLE]

Finally, it is instructive to consider $\mu_{\mathrm{Arrh.}}/\mu_{E}=\operatorname{var}N/N$. As discussed before, if the dominant thermal (collective) excitation of the system consists of $m$ visons, $\mu_{\mathrm{Arrh.}}$ is expected to be $m\mu_{E}$ at low temperatures, and so $\operatorname{var}N/N\approx m$. We can obtain this last result directly by considering that the number $\tilde{N}$ of collective excitations obeys a Poisson distribution at low temperatures and so $\operatorname{var}\tilde{N}=\tilde{N}$. The result then follows from $N=m\tilde{N}$.

Appendix E

Semiquantitative model of the partition function

The behaviour of vison number found in the simulations can be explained by the following three observations:

The energy cost of a nearest-neighbour dipole of visons is smaller than that of a single isolated vison. 2. 2.

An isolated vison polarises its surroundings, reducing the energy cost of nearby, aligned vison dipoles even further. 3. 3.

The energy cost of visons reduces as temperature increases because highly excited photon modes “wash them out.”

To demonstrate this, we construct a simplistic model of the partition function (55), which is semiquantitative at low temperatures and captures the salient features at large $T$. If we ignore interactions, other than the ones within two visons in a pair, the partition function factorises by diamond lattice sites: $Z=Z_{1}^{\mathcal{V}}$, where $\mathcal{V}$ is the number of such sites. Now, each site can host an isolated vison of either charge, or one half of a nearest-neighbour dipole. Assuming their energy cost is $E_{1}$ and $E_{2}$, respectively, we can write down a first approximation to $Z_{1}$ as

[TABLE]

where $E_{2}<E_{1}$ in line with observation 1.

Observation 2 states that the energy cost $E_{2}$ of a dipole near an isolated vison with the right alignment is reduced to $E_{2}-E_{\mathrm{int}}$ while that of a dipole on the same bond but with opposite orientation is increased to $E_{2}+E_{\mathrm{int}}$. This changes the partition function of the bond by

[TABLE]

In a crude approximation, we assume that each isolated vison introduces a fixed $E_{\mathrm{int}}$ to $m$ nearby bonds and that the resulting $\Delta Z_{1}$ can be factored into the partition function of the isolated vison. This gives

[TABLE]

Observation 3 concerns the strong interactions between a highly excited photon bath and the visons, a full treatment of which is a tall order. To estimate its effect on the vison thermodynamics, we propose a “Hartree–Fock approximation” where the effect of visons on the photon cloud is neglected, and we assume that the quadratic theory governing the photon modes at low temperatures Kwasigroch et al. (2017) holds at arbitrary $T$. In this approximation, the energy associated with the gradient component $B_{\mathrm{grad}}$ of the magnetic field is

[TABLE]

since $B_{\mathrm{curl}}$ on each site has a Gaussian distribution of variance $T/(2g)$. That is, all energy scales associated with $B_{\mathrm{grad}}$, and hence with the visons, are exponentially suppressed at high temperatures.

The partition function (62–64) can now be used to derive thermodynamic quantities. We plotted $\operatorname{var}N/N$ and the slope of the Arrhenius curve $\mathrm{d}(\log N)/\mathrm{d}\beta$ in Fig. 10; the results are in good qualitative agreement with the numerical simulation, cf. Fig. 7. Saturation occurs at a higher temperature than in the full large-$S$ treatment, albeit well below the zero-temperature energy cost of visons. More accurate estimates would likely follow from taking the emergent electric field and photon–photon interactions into account, but this is beyond the scope of the present paper.