Valence Bonds in Random Quantum Magnets: Theory and Application to YbMgGaO4

We analyze the effect of quenched disorder on spin-1/2 quantum magnets in which magnetic frustration promotes the formation of local singlets. Our results include a theory for 2d valence-bond solids subject to weak bond randomness, as well as extensions to stronger disorder regimes where we make connections with quantum spin liquids. We find, on various lattices, that the destruction of a valence-bond solid phase by weak quenched disorder leads inevitably to the nucleation of topological defects carrying spin-1/2 moments. This renormalizes the lattice into a strongly random spin network with interesting low-energy excitations. Similarly when short-ranged valence bonds would be pinned by stronger disorder, we find that this putative glass is unstable to defects that carry spin-1/2 magnetic moments, and whose residual interactions decide the ultimate low energy fate. Motivated by these results we conjecture Lieb-Schultz-Mattis-like restrictions on ground states for disordered magnets with spin-1/2 per statistical unit cell. These conjectures are supported by an argument for 1d spin chains. We apply insights from this study to the phenomenology of YbMgGaO$_4$, a recently discovered triangular lattice spin-1/2 insulator which was proposed to be a quantum spin liquid. We instead explore a description based on the present theory. Experimental signatures, including unusual specific heat, thermal conductivity, and dynamical structure factor, and their behavior in a magnetic field, are predicted from the theory, and compare favorably with existing measurements on YbMgGaO$_4$ and related materials.

Magnetic insulators often exhibit quenched disorder from material defects, including randomness in the strengths of magnetic exchanges. How this bond randomness interacts with geometrical frustration and quantum fluctuations of the moments is an interesting and largely open question, especially urgent in view of perplexing experiments that seem related to both features.
intrinsic for YbMgGaO 4 : its magnetic Yb 3+ layers are separated by layers of non-magnetic Mg 2+ and Ga 3+ ions which randomly occupy a single R3m crystallographic position [1,2,9], forming a triangular lattice with Mg/Ga occupancy described by a disordered Ising variable. The random electric field from the Mg/Ga site (which sits directly above or below a Yb-O-Yb oxygen site) can modify the Yb-O-Yb magnetic superexchange as well as the Yb effective spin-1/2 g-factor; experimentally, the distribution of g-factors [9] shows direct evidence for Hamiltonian randomness, though exchange randomness [12] has been difficult to quantify. . Each vortex hosts a protected spin-1/2 in its core. A particular type of domain wall, here shown separating the sets of domains on the left and the right, also carries dangling spin-1/2. When bond randomness destroys long-range VBS lattice-symmetry-breaking order it also nucleates a random network of these defects [bottom]: RG flows at low energies can produce longer-range random singlets, larger-spin clusters and spin glass freezing of the nucleated spins.
In YbMgGaO 4 , neutron scattering [5][6][7][8] and µSR [3] studies at temperatures down to around 50 mK (less than two percent of the Curie-Weiss temperature θ CW ) found no signs of magnetic ordering or of frozen moments. Even more unusually, careful comparisons with the nonmagnetic analogue LuMgGaO 4 exposed a peculiar power-law C(T ) ∼ T 0.7 heat capacity of the magnetic moments, 1 extending over a decade in temperature from 1 K down to 60 mK [1]. This fractional power law inspired theoretical work [6,[13][14][15][16] that interpreted these observations as evidence for a spin-liquid phase with a spinon Fermi surface. However a comparison of inelastic neutron scattering [5][6][7][8]10] at low and high temperatures shows little direct evidence 2 of a Fermi surface, and even more strikingly measurements of thermal conductivity [4] found that κ/T vanishes with temperature, substantially complicating any interpretation in terms of itinerant spinons.
This complex phenomenology brings to mind doped semiconductors [17,18] such as Si:P, where a broad distribution of couplings between magnetic moments, at randomly located dopant sites, was argued to generate a random singlet phase [19][20][21][22] -in which each spin i forms a singlet with another spin j, with weak singlets forming across arbitrarily large distances. The random singlet phase exists in 1d, where it is tractable using the strong disorder renormalization group, which iteratively integrates out the strongest couplings. In higher dimensions the above cartoon is not accurate for disordered spin networks: numerical strong disorder renormalization group calculations show that there is formation of higherspin clusters as well as long-range singlets, and likely spin glass freezing at the lowest energies (where strong disorder RG breaks down). [20,[23][24][25] However, if the initial coupling distribution is parametrically broad, strong disorder RG will capture the physics over a parametrically large range of lengthscales. Thus this is a mechanism for power law heat capacity without frozen moments over a large energy range, as is indeed experimentally observed in doped semiconductors.
Motivated by these developments, in this paper we study the effects of disorder on quantum paramagnetic phases that arise in frustrated quantum magnets. Our main concern is the interplay between (1) valence bond physics and (2) randomness in the exchange couplings.

A. Setting
Consider a magnet which, in the clean limit, is in a spin-gapped paramagnetic phase. It is useful to have a heuristic picture in mind for the various energy scales. 3 The spin gap ∆ S , which is the energy scale for creating spin-carrying excitations, may be loosely associated with the energy scale of formation of singlet valence bonds between pairs of local moments. At a second energy scale ∆ V B , which we may loosely identify with the energy scale for valence-bond-type excitations that do not carry spin, these valence bonds may freeze into a Valence Bond Solid (VBS) phase that breaks lattice symmetries, or they may form a quantum liquid. The latter phase is a gapped quantum spin liquid known as a short-ranged Resonating Valence Bond (RVB) state. The ratio ∆ V B /∆ S is typically of order one, but it will sometimes be instructive to consider the limit where the scales become well separated, ∆ V B /∆ S 1. What is the fate of such a system when the exchange couplings are random? We will restrict attention to the situation in which the width ∆ of the distribution of the random exchanges is smaller than the spin gap ∆ S of the clean magnet, but will consider both the case where randomness is parametrically weak (∆ ∆ V B , ∆ S ), and the case where randomness is weak compared to the spin gap, but has a strong pinning effect on the valence bonds (loosely speaking, ∆ V B ∆ ∆ S ). If the randomness ∆ is the weakest energy scale in the problem (i.e ∆ ∆ V B , ∆ S ) then it has strikingly different effects on VBS and RVB phases. RVB phases are stable. However VBS phases in d ≤ 2 are unstable, as we review below, because disorder couples linearly to the VBS order parameter. We will study the eventual fate of the disordered VBS state in various examples in We find that in the thermodynamic limit such a valence bond glass is unstable to the nucleation of defect monomers that carry spin-1/2. The RG flow from the resulting random spin network then determines the low energy physics, as in Fig. 1. both 1d and 2d.
In the weak disorder limit, we introduce a theoretical mechanism by which a regime of strong randomness -involving broad coupling distributions and random connectivities -can arise in 2d quantum magnets even when disorder is weak at the lattice scale. The lowenergy physics is determined by a multi-stage RG flow that begins with the nucleation of topological defects. These carry spin, despite the fact that the disorder scale is small compared to the spin gap of the clean system, ∆ ∆ S . This flow has several well-defined regimes when the initial disorder is weak. We characterize these regimes on the triangular and square lattices, showing that gapless spinful excitations necessarily emerge. At the very largest scales, the defect spins can break spin symmetry, for example via spin-glass order.
We will also study the case ∆ V B ∆ ∆ S , where randomness strongly pins the valence bonds. This limit is well modeled by a dimer model with random energies. We will show that broadly similar results obtain in both the weakly disordered VBS model and the random dimer model. At first sight one might have guessed that strong pinning would allow a distinct paramagnetic phase made only of randomly pinned singlets (a 'valence bond glass'), but we show that this phase is unstable to the nucleation of spinful defects in 2d.
These analyses raise interesting questions about what kinds of states can exist, even in principle, for disordered magnets on various lattices. Our results in the two limits of disorder strength motivate conjectures in the spirit of the Lieb-Schultz-Mattis theorem of clean magnets. Our conjectures rule out certain kinds of spin-gapped states for disordered Hamiltonians that preserve spin symmetry and preserve lattice symmetry on average, with spin-1/2 per unit cell of the statistical translation symmetry. We describe some constraints that we expect in physical terms, without aspiring to mathematical rigour. We propose that a useful perspective on Lieb-Schultz-Mattis-like constraints in disordered systems is from the standpoint of the entanglement structure of ground states.
The picture that emerges for the pinning of valence bonds suggests an application of the above ideas to YbMgGaO 4 (which we discuss in detail) and other frustrated materials.
Indeed an isostructural compound, YbZnGaO 4 , was recently studied [12] and shown to exhibit phenomenology closely related to that of YbMgGaO 4 . The YbZnGaO 4 dynamical structure factor S(q, ω) shows similar features to YbMgGaO 4 ; thermal conductivity κ(T ) is close to indistinguishable; and, remarkably, YbZnGaO 4 was also found to show a power law magnetic specific heat with an anomalous exponent, here C(T ) ∼ T 0.59 . Our random-singlet-inspired picture, in which the anomalous C(T ) exponent is an effective exponent over a range of temperature scales, and can take different values depending on the disorder distribution, appears to be a good candidate for the phenomenological description.
The importance of disorder in the YbMgGaO 4 family of compounds was discussed for g-factor randomness [9] as well as for exchange randomness: in particular the study of YbZnGaO 4 also reported new measurements of AC susceptibility for both compounds [12], which find broad excitation spectra suggested to be associated with spin freezing at ultralow temperatures T f < 0.1 K. This suggests that the interplay of disorder with frustration should be part of the description of these magnets. Here we consider this interplay together with the additional ingredient of spin-1/2 quantum fluctuations and singlet formation. Interestingly, the signatures of spin freezing in the two Yb magnets were reported to be distinct from conventional spin glasses: T f in both compounds is ∼ 20 times smaller than the temperature of the peak in specific heat, in contrast to conventional spin glasses, and the spin freezing entropy was also reported to be anomalously small [12]. As elaborated on below, ultralow temperature spin freezing of remnant defect spins is indeed the likely RG fate for the defect spin network that arises on the triangular lattice, suggesting the above ideas may be on the right track for these materials.
In the literature, short-ranged random singlet type phases have been proposed [26][27][28][29] for lattice magnets in numerical studies that found an apparent suppression of ordered phases under strong bond randomness. The suppression occurs when disorder strength becomes of order unity, i.e. couplings distributed across {J±∆} seem to require 3 a distribution width 2∆ J. Bipartite antiferromagnets are argued to remain stable through large ∆. [30,31] We note that recent work [32] on Hamiltonians relevant to YbMgGaO 4 found that adding a particular spin-orbit-coupled term (J ±± ) with small magnitude but random sign {±∆} melts the orientational part of the magnetic order. Relatedly, phenomenologically adding disorder with 2∆∼J in a spin wave theory was found [12] to capture the continuum features of the structure factor. The structure factor features were also discussed in the context of short ranged singlets [7]. Based on our analysis and the conjectured disordered-LSM restriction, if a short ranged valence bond glass is formed, for example by adding strong randomness to a magnetic phase, then in the thermodynamic limit it would become unstable to spin-1/2 defects, which would exhibit the strong-disorder physics discussed above.
Random singlet physics has been proposed to describe the dilute impurity spins (Cu-Zn substitution sites) in the spin liquid candidate Herbertsmithite [33,34]. It was also studied theoretically [35] on random graphs with fixed connectivity (related to the Bethe lattice), in a large-N limit, which found that spin excitations were gapless ("pseudo-gapped"). It is not obvious how a parametrically broad distribution of 'bare' coupling strengths can naturally arise in crystalline lattice magnets without dilution, so at first glance it is not clear how this kind of strong disorder physics can play a role. Counterintuitively, we will show that a broad-randomness phenomenology can arise even in a tractable weak disorder regime for certain 2D lattice magnets.

B. Overview
We will start in the following section (Sec. II) with VBS phases subjected to weak disorder. The weak disorder limit ∆ ∆ S , ∆ V B provides a clear separation of length scales, and allows predictions that are independent of the precise choice of microscopic Hamiltonian. We will argue subsequently that many phenomenological conclusions are relevant also in the intermediate disorder limit ∆ V B ∆ ∆ S . Recall that VBS phases have no magnetic moment and preserve time reversal symmetry, but spontaneously break some crystal symmetries 4 (for a simple example see Fig. 1). Various VBS phases have been observed numerically [37][38][39] as well as experimentally [40][41][42][43][44] in quantum magnets whose classical magnetic configurations are frustrated.
Since the VBS order parameter transforms under lattice symmetries, bond randomness couples to it as a random field. Consequently VBS phases are highly sensitive to disorder. The fate of the system may be analyzed in successive steps, associated with distinct scales of renormalization group flow.
The first step is the destruction of long range VBS order by disorder. The second step, at lower energies, involves topological defects in the VBS order parameter: in the simplest cases these are vortices which carry spin-1/2 (in some cases low energy states in domain walls are also important). These moments are nucleated even though the scale ∆ of disorder is much smaller than the spin gap ∆ S : the local disorder is not strong enough to break a singlet, but the rigidity of the VBS order parameter allows the effects of disorder to accumulate over large scales. The spin-1/2 defects form a random network with broadly distributed couplings, leading to formation of long range singlets and spin clusters which is captured, at least up to a large lengthscale, by a strong disorder RG.
In Sec. II A we warm up by reviewing, and clarifying some aspects of, this physics in 1d where it is known that the random singlet phase is one possible fate. In 2d, the final long length scale physics is subtle and depends on the details of the lattice, as we discuss for the square lattice in Sec. II B and the triangular lattice in Secs. II C and II D. At the longest length scales the strong disorder RG breaks down (in 2d), with the likely fate of the system being very weak spin glass order (on the triangular lattice) or Néel order (on the square lattice). We will substantiate this picture in the following. The theory enables strong randomness to arise from weak disorder through a multi-scale renormalization group flow, which is initiated by magnetic frustration at the lattice scale and remains protected by geometrical frustration down to ultra-low energy scales.
In Sec. II D we provide an alternate point of view on these results by formulating a quantum Landau-Ginzburg-like framework to discuss the effects of disorder on a valence-bond solid. This theory is not a standard Landau-Ginzburg theory, which is an expansion in powers of the order parameter and its gradients. Such a 'naive' Landau-Ginzburg expansion should really be thought of as an expansion about a trivial symmetrypreserving state. In quantum magnets with, say, an odd number of spin-1/2 per unit cell, Lieb-Schultz-Mattis (LSM) theorems [45][46][47] prohibit such a trivial symmetry preserving phase. Thus the standard Landau-Ginzburg expansion is not available for such a quantum magnet.
How can we construct a useful effective theory for a quantum magnet in the absence of a trivial symmetry preserving phase? There is a well-known alternative [48] that we will briefly describe and utilize. Though there is no trivial symmetry-preserving state for a magnet with spin-1/2 per unit cell, there do exist symmetry-preserving gapped phases. These states necessarily have topological order. The simplest one that is consistent with LSM restrictions is a Z 2 quantum spin liquid. 5 The best we can then do, in the spirit of Landau-Ginzburg theory, is to develop an effective 'quantum' Landau-Ginzburg theory for the quantum magnet in terms of the excitations of this Z 2 quantum spin liquid. This takes the form of a Z 2 gauge theory coupled to Z 2 charges. In the presence of spin SO(3) and time reversal symmetries, the symmetry properties of the excitations of the Z 2 spin liquid are severely constrained. [49][50][51][52][53] In Sec. II D we use the resulting 'quantum' Landau-Ginzburg theory as a framework to discuss disorder effects.
In Sec. III we address the intermediate-disorder regime (Fig. 2). We use the dimer model alluded to above to discuss pinning of singlets (Sec. III A) and to argue that a defect-free 'valence-bond glass' is unstable. We confirm using the 'quantum' Landau-Ginzburg theory that the results hold more generally for a strongly pinned VBS order parameter (Sec. III B). We study the random dimer model on the triangular lattice numerically, and provide a simple explanation for the necessity of a finite density of monomers in this limit. Our results for the effect of 5 Its universal physics is described by a deconfined Z 2 gauge theory. disorder on triangular-lattice classical dimers may also be of independent interest.
Using insights from these analyses, in Sec. IV we present some natural conjectures regarding LSM-like restrictions on the possible ground states of disordered magnets with spin symmetry and statistical translation symmetry.
In Sec. V we discuss the phenomenology of experimental observables from the standpoint described in this paper. For concreteness we focus on YbMgGaO 4 though the results are more general. This material or any other of course may not be in the controlled theoretical limit ∆ ∆ S , and indeed here there is no experimental access to the clean limit that is the conceptually useful starting point in the theory. However we note that a quantum paramagnetic regime has been found [32,54,55] in numerical work on the triangular lattice Heisenberg model for weak second-neighbor exchange J 2 0.07J 1 . Assuming a qualitative theoretical description motivated by our results in the theoretically controlled limits, in terms of randomly pinned short range singlets together with larger scale excitations arising from defects, we discuss the expected phenomenology for a variety of experimental probes: specific heat C(T ), magnetic susceptibility χ(T ), thermal conductivity κ(T ), and dynamical spin structure factor S(q, ω), as well as their behavior in an applied magnetic field. We also briefly discuss NMR and µSR. We compare all of these theoretical expectations to available experiments. We also discuss the isostructural compound YbZnGaO 4 and various other material candidates for investigating random valence bond physics.
We conclude with a summary of our results, and a discussion of the questions raised, in Sec. VI. Several Appendices contain additional details.

II. DISORDERING A VALENCE BOND SOLID
We consider spin Hamiltonians featuring both frustration and bond randomness. We begin by considering the simple case where the Hamiltonian has spin SO(3) invariance, 6 the bond randomness scale ∆J is much weaker than the mean valueJ.
We assume that the ground state is a paramagnetic VBS state when ∆J = 0. In contrast to conventional magnetically ordered phases, where bond randomness is irrelevant, the linear coupling to the VBS order parameter means that VBS order is destroyed even by arbitrarily weak randomness in dimensions d ≤ 2. The Imry-Ma argument [56] shows that domain walls in the VBS appear on a length scale ξ ∼ (J 2 /∆J 2 ) 1/(2−d) when d < 2. In d = 2 a more sophisticated real space renormalization group treatment of the domain walls [57,58] yields a length scale ξ 2d ∼ exp[J 2 /∆J 2 ] at which the long range order is broken up.
Let us now consider the physics that results at lower energies or on length scales larger than ξ. Since the order parameter is pinned by the coarse-grained random field on these scales, one might at first sight expect that there will be no remaining modes at low energies. This would be the case if the order parameter was, say, the Ising spin in the random field Ising model (modulo rare region effects). But, we argue below, this 'trivial' RG endpoint is forbidden on topological grounds for a VBS order parameter on the square or triangular lattice. The random field necessarily introduces topological defects in the VBS order parameter, in particular certain vortex defects. Crucially, these topological defects carry spin-1/2 degrees of freedom that are not bound into short-range singlets. The physics on scales ξ 2d is that of a random network of these defect spins, leading to a strongly disordered regime even when the bare disorder strength ∆J is small. The low temperature response is then dominated by this network of broadly distributed defect spins.
Before describing this physics for triangular lattice magnets, we consider 1d chains and 2d square lattice magnets, which are interesting in their own right.
A. Spin-1/2 chain A spin-1/2 chain in the spontaneously dimerized phase shows the simplest version of this physics. [59][60][61][62] We imagine perturbing such a Hamiltonian (e.g. the J 1 -J 2 chain with J 2 > 0.2411J 1 [63,64]) by adding weak, shortrange correlated bond randomness. A given realization of the randomness will break translational symmetry, but we assume that translational symmetry is preserved on average: i.e. that the probability distribution of the disorder is translation invariant.
Weak disorder induces static domain wall defects in the VBS order parameter with typical separation [56] ξ 1d ∼ J 2 /∆J 2 , so VBS order is lost at this scale. However each domain wall carries a single unpaired spin-1/2 moment. Virtual processes induce random couplings between the domain wall spins, with the strength of the coupling falling off exponentially in the domain wall separation r: schematically, |J eff | ∼ e −r/η , where η is a spin correlation length associated with the dimerized region. This leads to a 'renormalized' spin chain with random couplings.
The exponential sensitivity to separation, combined with the random locations of the defects, means that these couplings are very broadly distributed. A strong disorder RG approach is therefore appropriate for the next stage of RG flow. [22,65,66] What happens in this flow depends on the sign structure of the effective couplings, as we clarify below. In the simplest case all the effective couplings between adjacent defects are antiferromagnetic. As is well known, such a random AF chain flows to the random singlet fixed point. [22,65] Roughly speaking, this is a state where each spin is paired into a singlet with another spin which may be arbitrarily far away. This yields excitations at arbitrarily low energies that are associated with breaking weak, long-distance singlets.
Note that the above is in stark contrast to the fate of a 2-leg ladder with antiferromagnetic couplings. Again let us start in a columnar VBS phase, where parallel valence bonds form within each chain of the ladder. The VBS order parameter is again Ising-like, since there are two degenerate VBS tilings. However domain walls between these two tilings now host two unpaired spins. These spins will form a singlet, so that no modes survive to longer lengthscales. This difference between the single chain and the 2-leg ladder is indicative of a Lieb-Schultz-Mattis-like constraint on disordered magnets that we will return to later. Now let us discuss the signs of the effective couplings J eff in more detail. The essential question is whether the signs of the effective couplings J eff are deterministic or random in the weak disorder limit ∆J/J 1. In fact, both scenarios are possible without fine-tuning. Which one occurs depends on the nature of the spin correlations in the VBS phase of the clean system, as we discuss in Appendix 1. The amplitude of these spin correlations of course decays exponentially with distance, but the sign structure may be either commensurate or incommensurate with the lattice. The two possibilities are exemplified in the J 1 − J 2 chain. In the region of the VBS phase close to the phase boundary to the gapless phase, for J 2 /J 1 0.2411, the sign of S(0). S(r) alternates with period 2. On the other hand spin correlations become incommensurate for J 2 /J 1 0.52. [63] In the former case J eff between two adjacent defects, which necessarily occupy opposite sublattices, will always be antiferromagnetic. Interestingly, these couplings are also guaranteed to be antiferromagnetic for a class of Hamiltonians that permit a sign-free Monte-Carlo treatment: a general argument for this is in Appendix 1. On the other hand, when the spin correlations of the clean system are incommensurate, the large random separations of defects imply that the signs of the couplings J eff will be independently random in the limit of weak disorder, i.e. the limit of large ξ 1d (Appendix. 1).
The two possibilities lead to different behaviour at large lengthscales when the VBS is perturbed by weak disorder. As noted above, the antiferromagnetic case leads to the random singlet phase. Conversely, when ferromagnetic couplings are present, strong disorder RG shows that a different fixed point is reached, involving the generation of large effective spins during the RG procedure. [67,68] (This fixed point is not at infinite randomness, so the SDRG treatment of it is not strictly controlled.) Physically, the resulting phase should perhaps be thought of as having 'quasi-long range' spin glass order. We will encounter related issues in two dimensions.
Some of these theoretical expectations for J eff are manifested experimentally in quasi-one dimensional S=1/2 two-leg-ladder materials [69,70] such as SrCu 2 O 3 and BiCu 2 PO 6 . Hole doping via Cu to Zn substitution, at the level of less than a few percent, shows a signal of spin freezing into a pattern with three dimensional antiferromagnetic Néel correlations. The ordering occurs among the moments that are induced by a nonmagnetic impurity, which roughly speaking removes one of the Cu electrons from a singlet bond. Due to the well-defined sublattice sign structure, the coupling between these dangling moments is unfrustrated, producing antiferromagnetic correlations with spins up on one sublattice and down on the other.
B. Spin-1/2 square lattice Now let us consider square lattice spin-1/2 magnets in a columnar VBS phase. This fourfold-degenerate VBS pattern can be associated with four cardinal directions of a planar vector ϕ. In addition to domain walls, it admits a discrete Z 4 vortex defect which carries an unbound spin-1/2 in its core [71]. This vortex defect is a junction between four 'elementary' domain walls, across which ϕ rotates by π 2 . The VBS order will be disrupted at the Imry-Ma lengthscale ξ 2d where such domain walls appear. We make the natural assumption that 'elementary' domain walls are less energetically costly than 'composite' domain walls where ϕ flips sign.
One may argue that the breakup of the VBS into domains necessarily also introduces VBS vortices with typical separation ξ 2d . This is because, when ξ 2d is large, the core energy cost of a vortex is negligible in comparison with the typical energy cost of domain rearrangements on this scale. See Appendix 2 for more detail (similar arguments have been made in the context of the random field XY model [72]).
The spin−1/2 moments in the vortex cores will then determine the eventual fate of the system. As in 1d, a key feature of this system is that the distribution of couplings J eff for adjacent vortex spins is extremely broad, despite the fact that the bare disorder ∆J is weak. Since the magnitude |J eff | depends exponentially on the separation of the vortices, while the distribution of separations has mean and width both of order ξ 2d , the distribution of J eff is broad even on a logarithmic scale.
On the square lattice, a key geometrical fact is that core spins of Z 4 vortices and anti-vortices are associated with opposite sublattices. Therefore a given defect has a well-defined sublattice assignment, even when there are quantum fluctuations in its precise position. Further, we show in Appendix. 1 that for a natural class of models the sign of the effective interaction is determined solely by whether the two core spins are on the same or opposite sublattices. (This class includes, but is not limited to, sign-problem-free models for VBS phases such as the much-studied J-Q model [37].) The effective Hamiltonian for the defect spins then takes the form where r, r form a random (but correlated) selection of sites of the square lattice with typical separation ξ 2d , and J eff > 0 is exponentially small in the separation between r and r , with the sign r = ±1 for the two sublattices/vorticities. Note that, despite the strongly random sign and magnitude, this interaction is unfrustrated. For energetic reasons, the closest defect pairs will predominantly be of opposite vorticity (and thus antiferromagnetically coupled), so a significant fraction of the defect spins will be bound into singlets of size ∼ ξ 2d in the very first step of strong disorder RG. But in this 2d problem, the strong disorder RG is expected eventually to break down, due to a flow away from infinite randomness. [20,[23][24][25] The breakdown of the strong disorder RG, at late stages of the singlet and cluster formation process, does not automatically rule out nontrivial disordered fixed points at finite randomness. [73] 7 [74] However since the interaction here is unfrustrated, a very simple alternative scenario on the square lattice is that at long scales the disordered VBS could transform itself into a dilute very weakly Néel ordered state. (On frustrated lattices very weak spin-glass order is the natural possibility at the very longest scales, beyond the randomsinglet-like regime, as we discuss below.) As in 1d spin chains, depending on the initial clean Hamiltonian on the square lattice, we can also in principle obtain effective couplings with fully randomized signs and no sublattice sign structure (Appendix. 1); we discuss this case in the following section.
C. Spin-1/2 triangular lattice: columnar VBS We are now ready to consider the triangular lattice. We start with columnar VBS order, since this is the simplest to visualize. The infra-red fate in this case is fairly intricate, involving several lengthscales. In Sec. II D we comment on another VBS order with a larger unit cell. For the columnar VBS there are 12 symmetry-related ground states. These can be labelled by a pair of vectors (a : b). The first vector a specifies the direction of the valence bond out of some fixed base point r 0 on the lattice, and can take the values a i for i = 1, 2, 3 (the three vectors shown in Fig. 3a) and alsoā i ≡ −a i for i = 1, 2, 3. The second vector b specifies the axis along which the columns line up, and along which the pattern is invariant by a unit translation. Since the sign of b has no meaning, it is enough to let b be one of the three vectors a 1,2,3 . We will denote the various possibilities by (i : j) or (ī : j) for (a i : a j ) or (ā i : a j ) respectively. Since j = i, there are 12 columnar VBS patterns in total. For example the domains in Fig. 3(a,b,c) are labelled (1:3), (1:3), (1:2) respectively.
As we will see, it is natural to group the 12 VBS patterns into three sets of 4. We denote these sets [1,2], [2,3], [3,1]. These groupings have a simple geometric meaning: each set corresponds to one of three ways of viewing the triangular lattice as a square lattice with an extra diagonal bond. For example, removing all the bonds parallel to a 3 yields one possible square lattice. The corresponding VBS patterns (a : b) are those that can be drawn on this square lattice: they make use only of the vectors a 1 and a 2 . This gives the set [1,2], which includes the 4 VBS patterns (1:2), (2:1), (1:2) and (2:1).
We distinguish two kinds of domain walls. By the above construction, domain walls between patterns within a given set [i, j] are mapped to the VBS domain walls we already encountered on the square lattice. They can similarly be identified as composite and ele-mentary, with four elementary domain walls emanating from a VBS vortex. Such a vortex is shown in Fig. 3(e). Each type of domain allows a straight boundary with two possible orientations. These allowed orientations are the same for all the domains belonging to a given set [i, j], and are parallel to the vectors a i and a j . This means that an elementary 'intra-set' domain wall can alternate between these two orientations without incurring unpaired dangling spins. As a result, after coarse-graining these domain walls can take an arbitrary path without incurring 'frustration'.
Domain walls between patterns from different sets, which we dub super domain walls, are more complicated. An example can be seen in Fig. 3(d). Each superdomain wall generically has a finite density of dangling unpaired spins along its length. More precisely, it can avoid this fate only when it lies exactly parallel to a single preferred orientation. A superdomain wall between superdomains [i, j] and [i, k] can be 'unfrustrated' only if it lies parallel to the lattice vector a i that is shared by the two superdomains. For the superdomain wall in Fig. 3(d), the unfrustrated orientation would be along a 2 . A nonstraight path, as is required generically between two superdomains, must deviate from this single allowed orientation, and hence will nucleate a finite density of dangling spins along the path. This contrasts with the intra-set domain walls, which are able to avoid dangling spins even if they are not straight. Now let us turn on bond randomness. The VBS pattern will fragment into domains at a scale ξ 2d . This is the lengthscale at which the 'cheapest' domain walls proliferate. The most natural possibility is that these are the 'unfrustrated' elementary domain walls within a given set [i, j]. The superdomain walls will then proliferate only at a parametrically longer length scale ξ S ξ 2d . More precisely, the line tension of domain walls decreases under coarse graining, due to energy optimization from short-scale disorder pinning. [58] (A detailed RG treatment in the present case would have to take account of the anisotropy of the line tension, with different preferred directions for different domain walls.) ξ 2d is the lengthscale where the line tension of the cheapest domain walls can no longer compete with the energy gain from disorder when domains are introduced. On longer scales, the line tension of the more expensive superdomain walls continues to renormalize downwards, and superdomains eventually appear on the lengthscale ξ S . Since the Imry-Ma lengthscale depends exponentially on the ratio of the bare line tension to disorder, the larger bare line tension for superdomain walls means that ξ S /ξ 2d will be exponentially large in the limit of weak disorder.
This results in a range of scales ξ 2d L ξ S where (after averaging over scales of order ξ 2d ) translation and reflection symmetries are restored, but 60 degree lattice rotation remains broken. On this scale there is a patchwork of domains all belonging to one set [i, j], for example the set [1,2]. In this regime there is therefore long range order for a lattice nematic order parameter which picks out one of three unoriented vectors a i (in this example a 3 , i.e. the vector not included in the set [1,2]). This phenomenon, where lattice nematic order is longer ranged than lattice translation order, has been dubbed 'vestigial order'. [75] The patchwork of domains from a single set [i, j] can be related to domains on the square lattice, as noted above. As on the square lattice, this leads to vortices which have spin-1/2 moments at their cores, and a broad distribution of exchange strengths for the effective interactions J eff between defects, What are the signs of these couplings? For the triangular lattice, fully randomized signs seem the most natural possibility, since this is what happens when the VBS order in the clean system is not too strong: this is argued for heuristically in Appendix 1. In this regime, the nonbipartiteness of the triangular lattice ensures the signs are scrambled for large ξ 2d . We will focus on this case in the following. (Somewhat surprisingly, there is in principle another possibility: for an appropriate choice of clean Hamiltonian, the couplings can have a definite sign structure within a given superdomain. 8 ) For the case of random signs, within a superdomain, We may contemplate a strong-disorder RG treatment for the vortex defects, in which the strongest couplings (either F or AF) are sequentially integrated out. The breadth of the distribution of J eff will ensure that this procedure is accurate up to a parametrically large lengthscale.
During the procedure the presence of F couplings will lead to the formation of large, randomly coupled effective spins [20,23,67,68], and the eventual fate is very likely 8 At a minimum this requires the VBS order to be strong, strongly breaking the symmetry of the triangular lattice locally. As we noted above, a given superdomain type (e.g. [1,2]) corresponds to a way of embedding a square lattice into the triangular lattice. If the 'extra' diagonal bonds were deleted we would be left with the square lattice problem. In this case there is a mechanism for stabilizing an unfrustrated sign structure for J eff which we discussed in Sec. II B. (This mechanism relied on the fact that vortices and antivortices are assigned to definite sublattices of the bipartite square lattice.) There is no obstacle in principle to this sign structure surviving when the diagonal bonds are reintroduced; this depends on nonuniversal features of the Hamiltonian (Appendix. 1). In this situation the ultimate fate of a superdomain may be 'Néel' ordered, with the Néel pattern inherited from the effective square lattice. The definition of the Néel order parameter differs between distinct superdomains, so on scales larger than ξ S the Néel-ordered superdomains are coupled together in a frustrated manner, likely leading eventually to a spin glass fixed point.
to be spin glass order. The corresponding lengthscales would require more detailed analysis.
In the present context (columnar VBS) we must also remember that on the scale ξ S the superdomain walls appear. These make up a network of 1d spin-1/2 chains, subject locally to weak disorder. These chains will also be gapless at the lowest energies, for example exhibiting the 1d 'large spin' phase discussed in Sec. II A. The relative abundance of different types of excitations (associated with vortex spins versus domain walls) at different energy scales is likely to depend on nonuniversal parameters such as the ratio ξ S /ξ 2 2d and the energetics on the domain walls.
Just as on the square lattice, the basic point is that the Imry-Ma instability of the VBS led to a proliferation of spinful defects. This is true also for other VBS patterns on the triangular lattice which are not so easily visualized.

D. Other VBS states via proximate Z2 spin liquid
The necessity of 'spinful' topological defects, which we saw directly at the lattice level for the columnar VBS, is a more general phenomenon. It is useful to develop a coarse-grained approach that can be applied to VBS orders with a more complex structure, e.g. a large unit cell, or strong quantum fluctuations, for which defects are not easily visualized at the microscopic level.
For ordered phases in classical magnets, the natural coarse-grained language is Landau-Ginzburg theory which is an expansion of the free energy density in powers of the order parameter and its gradients. As explained in the Introduction, for quantum magnets with spin-1/2 per unit cell, an analogous Landau-Ginzburg theory is not available. The best we can do is a more sophisticated 'quantum' Landau-Ginzburg theory 9 which uses a description in terms of the excitations of a gapped topologically ordered Z 2 quantum spin liquid.
Consider a gapped Z 2 spin liquid on the (clean) triangular lattice, with a spin-1/2 'spinon' excitation and a spinless 'vison' excitation v, both of which are bosonic. Each sees the other as a π flux (they are mutual semions). Condensing either of these excitations confines the other one, and destroys the topological order. The spinon and vison, together with their bound state, which is a fermion, exhaust the three nontrivial anyon types in the spin liquid.
In any such spin liquid the vison field v must transform nontrivially under lattice symmetries, for a reason touched on above. If v was trivial under lattice symmetries, then a symmetric gapped state, without topological order, could be obtained by condensing v. This would be inconsistent with the LSM theorem. When the spin Hamiltonian has spin SO(3) symmetry, it is known, in fact, that the vison -defined to live on the dual honeycomb lattice -sees a background π flux through each plaquette. [48][49][50][51][52][53]76] Physically one can simply think of the microscopic spin-1/2 magnet as having one 'background' spinon at each site of the triangular lattice. The visons see this spinon as π flux. This restriction then completely fixes the action of lattice symmetries on the vison field. Lattice translation, for example, acts projectively.
An immediate result of the nontrivial action of lattice symmetry on v is that its condensation breaks spatial symmetry to give a VBS phase. Therefore a quantum Landau-Ginzburg theory for v allows us to discuss the VBS order and its defects. (We will put this theory on a lattice, since this makes it easier to handle the Z 2 gauge structure.) Starting from the spin liquid, a variety of VBS phases can be described by various condensation patterns for v (See Appendix 3 for details). For concreteness we will restrict to a specific one.
A mean field treatment of the vison field on the triangular lattice suggests that one natural possibility for the condensed state is a plaquette VBS with a 12 site unit cell. [76,77] In this continuum treatment, the vison field is viewed as a 4-component vector v, whose four components arise from four low-lying modes in the vison's microscopic dispersion. 10 The field v changes sign under Z 2 gauge transformations, and the physical VBS order parameter is a bilinear in v. There is a discrete set of 48 possible directions for ordering of the vector v and 24 for the physical VBS order parameter.
Let us write down a schematic energy functional for the vison. We are interested in a regime where quantum fluctuations of the (pinned) VBS order are irrelevant, so it suffices to think about minimizing a classical energy. For convenience we regularize the vison condensate vector v on a (fictitious) coarse-grained lattice, giving a description for the low energy vison field v through a lattice gauge theory: The t term captures the stiffness of the VBS order (and is invariant under the gauge The d term is the simplest gauge invariant coupling to weak disorder. The precise form of the interactions is not important, but the ellipses . . . must include the anisotropy terms which select out a set of discrete ordering directions for v, corresponding to the discrete set of degenerate VBS states. The . . . may also include an energy penalty for gauge flux excitations, which are plaquettes of the lattice in which the background flux seen by the bare vison field (π flux per spin-1/2 site of the original lattice) is modified by the presence of an additional π flux, seen by the low-energy vison condensate field v. In the gauge theory language these gauge flux excitations are plaquettes of the lattice where the product of the gauge fields on the links is equal to −1: These gauge fluxes, seen by the vison, are precisely the spinon excitations of the spin liquid. Therefore Eq. 5 must be supplemented with the crucial information that these spinons carry spin-1/2. The interactions of these spinful degrees of freedom are of course neglected in the classical Hamiltonian above. We now consider topological defects in the VBS order parameter, showing that appropriate defects bind a spinon (see also Appendix 3). It is simplest first to consider an artificial continuum limit where the anisotropies are switched off, and v is allowed to 'order' anywhere on the sphere S 3 . The VBS order parameter then lives in projective space, RP 3 = S 3 /Z 2 . (The above energy functional is then simply the standard Ising gauge theory representation for an order parameter in projective space [79].) When the disorder d is small in comparison to the stiffness t, this order parameter varies slowly and smoothly, except at the locations of point vortex defects. These are allowed as a result of the nontrivial homotopy group π 1 (RP 3 ) = Z 2 .
As we traverse a loop surrounding such a defect, the VBS order parameter smoothly traverses the topologically nontrivial cycle in RP 3 . But the hallmark of such a topologically nontrivial trajectory is that the vison field v acquires a minus sign on traversing the loop around the defect -i.e. that v has a branch cut ending at the defect. Energetics dictates that this branch point terminates at a gauge flux. This is because, for the terms σ ij v i . v j to remain positive on the links ij that cross the branch cut, σ ij must be −1. The termination of this line of negative σs is a plaquette where σ = −1. Therefore the presence of RP 3 vortex defects is equivalent to the presence of gauge fluxes, i.e. spin-1/2 spinons.
Returning to the case with nonzero anisotropy, these gauge fluxes are inevitable on the Imry-Ma lengthscale. Fluxes necessarily accompany appropriate discrete 'vortices' in the VBS order parameter, see below. These vortices will necessarily be nucleated since as noted above (Appendix 2; see also Ref. [72]) the core energy of a defect cannot compete with the large energy scale associated with domain rearrangements on the Imry-Ma lengthscale. 11 Roughly speaking, a vortex will be nucle-ated when the random field, coarse-grained on this scale, itself displays such a vortex configuration.
In the weak disorder limit, the discreteness of the VBS order parameter must be taken into account in classifying defects. (In a given microscopic model the anisotropies that select the discrete ordering directions might happen to be weak, manifesting themselves only on a large lengthscale. But in the weak disorder limit this lengthscale is still much smaller than the size of domains, so discreteness is important.) In the discrete setting a defect is a junction between k > 2 domains, and types of defect are labelled by the list of ordering directions v (1) , . . . , v (k) that surround the defect, in anticlockwise order. The set of possible defects depends on the set of domain walls that are generated in the Imry-Ma process. So long as each type of domain wall has a unique gapped ground state when regarded as an effective 1D system, one may argue heuristically that the flux trapped at the defect is fully determined by the defect type. In the above model, a nontrivial flux is trapped whenever the sign of the gauge invariant quantity One may check explicitly in this model that such defects indeed exist. These fluxes are the discrete analogues of the RP 3 vortices, and trap localized spinons.

III. STRONG PINNING OF VALENCE BONDS
Unpaired spins are energetically expensive. Nevertheless, in the controlled limit of weak disorder we saw that defects in the singlet pattern were inevitable on the Imry-Ma lengthscale. As a result, the ground state could not remain paramagnetic on the longest scales. A ground state made only of short-range, static singlets was impossible.
In this section we argue that this remains true, in two dimensions, in the regime where the disorder is no longer weak, and has a strong pinning effect on valence bonds. More precisely, we argue that a state of static short-range singlets without defects -which is what we refer to as a 'valence bond glass' -is not a possible ground state for a square or triangular lattice magnet with short-range correlated disorder. Defects are again nucleated, and the natural possibilities for the physics at the longest lengthscales are again those discussed in the previous section.
In Sec. IV we use this result and that of the previous section to motivate Lieb-Schultz-Mattis-like conjectures on allowed ground states for magnets with 1/2-odd integer spin per unit cell.

A. Instability of valence bond glass to defects
When the spin gap is finite, we expect that the ground state can be described in terms of valence bonds with a finite typical size. In general the ground state is a superposition of valence bond configurations. Here we discuss the extreme limit of valence-bond pinning, in which quantum fluctuations are completely suppressed, and the ground state is a single frozen configuration of nearestneighbour singlets. In this limit there is no meaningful local VBS order parameter (but see also Sec. III B).
In the present regime the system can be mapped into a classical dimer model with random bond energies, at zero temperature. Hard-core dimers living on bonds of the lattice represent the singlets. Each lattice bond is assigned a random energy, which we draw from a distribution (e.g. uniform and bounded) with width of order ∆. In the putative 'paramagnetic' state, which would correspond to a valence bond glass, the dimer covering is complete (every site belongs to a dimer) and the energy is the sum of energies of dimer-occupied bonds. We also allow unpaired monomer sites, at a large energy cost K.
In a given finite sample, the ground state when K = ∞ is a unique complete dimer covering selected by the disorder. This energy minimization problem is nontrivial due to the dimer constraint. We ask whether this state is stable when monomers are allowed with a large but finite cost K.
The universal properties of the K = ∞ state depend on whether the lattice is bipartite. For bipartite lattices such as the square lattice the question of stability has been addressed, [80,81] and is closely connected to the stability of the random field XY model to vortices [82] and the stability of the Bragg glass [83] describing an elastic medium subject to pinning. Numerical studies in this context [80,81] found that, while introducing a fixed monomer in a finite system costs a positive average energy δE > 0 which grows without bound with system size, the standard deviation of δE also grows. The net result is that the typical energy cost of an optimally placed defect in a system of size L is negative at large size and of order δE opt ∼ −(log L) 3/2 . At large L this overwhelms the core cost K, so that defects are nucleated and the 'Bragg glass' is destroyed.
For nonbipartite lattices the mapping to an elastic medium does not apply, and stability of the pinned 'dimer glass' does not appear to have been studied. We first give a numerical treatment and then a very simple theoretical explanation.
To determine the stability of the dimerized state we must determine the probability distribution of the energy cost δE for introducing a defect (δE is measured relative to the defect-free ground state, and may be negative). We have simulated the model on the triangular lattice and studied the behavior of monomers. Results are shown in Fig. 4 (see Appendix 4 for details and comparison with the square lattice).
Fixed defects on the triangular lattice are found to have an energy distribution which converges to a fixed form as L → ∞, with mean and width of order one, unlike the square lattice case where the mean and width diverge with L. Varying K simply shifts the mean of this distribution. If the tail of the distribution extends to arbitrarily negative values, this will be sufficient to guarantee nucleation of defects, as we discuss below.
Defect nucleation in the random dimer model. We obtained optimal dimer coverings numerically for an L × L triangular lattice, with random bond energies drawn from the interval [−1, 1]. Monomer defects were introduced in two ways: (a) at fixed sites (blue triangles, top); and (b) at sites chosen, from L possibilities, to optimize the defect energy cost δE (red circles, bottom). The histograms of the δE distributions for various system sizes are shown on the right. Left: their means (solid curves) and standard deviations (dashed). (a) Fixed defects typically cost energy δE > 0: the distribution of δE converges to a fixed L = ∞ shape, with a nonzero tail at δE < 0 (shaded). (b) Optimized defects occupy these rare regions, giving an energy gain −δE > 0 which appears to increase without bound as L → ∞.
For a direct check, we compute the energy of a 'partially optimized' defect pair, where the optimization is over O(L) configurations. We find that whereas the mean energy of unoptimized defects is finite as L → ∞, the optimized energy becomes increasingly negative at large L, in a manner consistent with the form δE opt ∼ −(log L) 1/2 .
Thus despite the differences between the two classes of geometries, the ultimate conclusion is the same. By allowing previously-impossible dimer configurations that better conform to the disordered potential landscape, optimally placed monomers reduce the system's energy at sufficiently large sizes (no matter how large the core energy) yielding a finite density of spin-1/2 defects in the thermodynamic limit.
These results can be simply interpreted. We conjecture that the universal properties of the defect free K = ∞ state are related by duality to those of the 2D Ising spin glass, which exists at zero temperature though not above. [84] One way to motivate this conjecture is to note that if we relax the constraint that one dimer touches each site, so as to allow any odd number of dimers, there is an exact duality mapping to a standard nearest-neigbour Ising spin glass on the honeycomb lattice. 12 More heuristically, note that, as in the Ising spin glass, the elementary excitations of the ground state are 12 Let the dimer occupation number on a link of the triangular lattice be represented with a variable µ = ±1 which is −1 if the link is occupied. Take the energy to be E = µ , with the random bond energies independent and symmetrically distributed around zero. Let µ 0 be an arbitrary fixed reference configuration. We may define an Ising spin configuration on the sites of supported on loops. In the Ising spin glass these (unoriented) loops are the boundaries of domains which are flipped relative to the ground state. Here, a loop excitation is supported on a chain of bonds which are alternately occupied and unoccupied, and consists in reversing the occupation of each bond on the loop compared to its value in the ground state.
In a system of size L, the lowest energy excitation is expected to be a loop of linear extent ∼ L which visits ∼ L d f bonds (with d f > 1), and which has energy L θ for negative θ. (We have computed the exponent d f = 1.28(1) numerically to confirm it is consistent with the Ising spin glass value [85][86][87].) The decrease of the energy L −|θ| with L is equivalent to the instability of this glass to finite temperature. [88,89] The negative value of θ also explains the L-independence of the distribution of energy of defects, as explained in Appendix 5. Introducing a defect introduces an open excitation string, whose energy can be parsed into contributions at different distances from the defect; the negative value of θ means that contributions from large distances are subleading.
The convergence of the distribution at large L means that introducing a defect at a typical site costs an energy ∼ K with fluctuations only of order one. At first sight one might think that this implied the state was stable for large K. However this is not correct, because we must take into account rare locations where the energy gain |δE| from the defect is much larger than average. If there exists a finite density ρ of rare locations with δE ∼ −K, then no matter how small this density, the glass will be destroyed on scales 1/ρ.
A simple rare region argument (Appendix 5) indeed shows that such locations always exist, with a density which is at least of order ρ min ∼ exp[−c(K/∆) 2 ] at large K, with c a numerical constant. This scaling arises because a rare region of area ∼ (K/∆) 2 is required to overcome the core energy of the defect. (This simple picture gives δE opt ∼ −(log L) 1/2 at large L.) As a result of the nonzero value of ρ, the glass is ultimately unstable. In our toy model the lengthscale for this instability can be made large by tuning K, but it is finite for any K < ∞.
In conclusion, the natural attempt to use disorder to produce a paramagnetic state, a valence bond glass without topological order or long range VBS order, does not work.
the dual honeycomb lattice (which are located at the centres of triangular plaquettes) satisfying S i S j = µ 0 µ , where i and j are the dual sites either side of link . If the number of dimers touching each triangular lattice site is constrained only to be odd, then S i = ±1 is unconstrained. Defining similarly J ij = µ 0 , the en- For the above distribution of bond energies, the J ij are independent and symmetrically distributed around zero, giving a standard nearest-neigbour spin glass.
In the previous section disorder was applied at the level of individual singlets. In this section we point out that the same conclusions hold in a slightly different limit, where local VBS order is formed on short scales and then pinned on slightly longer scales.
As discussed above, various VBS order parameters on the triangular lattice can be viewed as living on a manifold of the form RP n−1 = S n−1 /Z 2 . Let us consider the pinning of such an order parameter using the Ising gauge theory formalism described in Sec. II D [79]: As above, the n-component unit vector v ∈ S n−1 is a redundant parametrization of RP n−1 , and the classical energy functional must be supplemented with the information that a gauge flux (a plaquette where ij σ = −1) is accompanied by a spinon. We assign an energy cost K to such plaquettes. In contrast to Sec. II D we now take the disorder strength d to be large. 13 Consider the limit of strong pinning on the sites, so When spinons are forbidden by hand (at K = ∞) we can choose the gauge σ ij = 1. We then have an Ising model with effective couplings J ij = h i · h j . These couplings are randomly frustrated due to plaquettes where the gaugeinvariant quantity ij ( h i · h j ) has negative sign. Similar models, without the gauge field, have been studied in the context of Heisenberg models with random anisotropy [90,91]. Ising spin glass behaviour is expected at zero temperature [91]. This agrees with the above picture for the dimer model.
The stability to introduction of spinons is then the question of whether a large core energy cost K for gauge fluxes in σ can be outweighed by their efficacy (when optimally placed) in relieving frustration in J ij . As above, a simple rare region argument indicates that spinons will be nucleated, with a density at least exp In principle we could wonder about the possibility of nontrivial phases at weaker pinning that are not captured by the above treatment. But it seems unlikely that other stable phases exist in 2d.
We have found, in three separate regimes, that pinning of singlets by disorder inevitably nucleates spinful defects in 2d. In the previous sections we considered a VBS-ordered state subject to arbitrarily weak disorder.
In this section we have considered a limit with strong pinning of singlets on individual bonds, and also a limit in which there is a well-defined VBS order parameter on short lengthscales that is strongly pinned on lengthscales of the same order. In Sec. IV we will synthesize these observations.

IV. LIEB-SCHULTZ-MATTIS-LIKE CONJECTURES FOR DISORDERED MAGNETS
Above we started with a clean system in either a VBS phase, breaking translation symmetry, or a spin liquid phase with topological order. We then added enough disorder to remove symmetry-breaking VBS order via Imry-Ma, or to destroy topological order via vison condensation. Naively one might have expected the resulting phases to be spin-gapped and featureless (with no symmetry breaking or topological order). But this fate was averted by the appearance of the spin-1/2 defects.
In this section we suggest that this is a consequence of more general constraints on ground states of random magnets which can prevent them from having a 'trivial' disordered ground state. That is, some Lieb-Schultz-Mattis-like constraints survive in the random setting, so long as the Hamiltonian preserves translation symmetry 'on average'. Recall that for non-random (i.e. translationally invariant) magnets with an odd number of spin-1/2 per unit cell, the Lieb-Schultz-Mattis (LSM) theorem [45], extended by Oshikawa [46,92] and Hastings [47], implies that a unique gapped ground state is not possible in a system with periodic boundary conditions in the thermodynamic limit. In effect, the ground state must have a broken symmetry or topological order with its associated ground state degeneracy on the torus; or, otherwise, it must have gapless excitations.
Here we make some initial conjectures for magnets with quenched randomness in one, two and three dimensions, with spin-1/2 per unit cell. After defining the setting in Sec. IV A, in Sec. IV B we discuss the spectrum and correlation functions. We conjecture that in 2d a disordered paramagnet, without topological order or translational symmetry breaking, must have gapless spin excitations. We also show that in 1d, if statistical translational symmetry is not broken, the spin correlation function must decay sufficiently slowly in space (1/r α with α ≤ 2). We then suggest that 3d admits states that are rather different to those in lower dimensions. Finally, in Sec. IV D we propose that the entanglement structure of the ground state provides a useful alternative perspective.
The importance of 'statistical' lattice symmetries in protecting topological insulator surface states from localization was appreciated in Refs. [93,94]. Also, Ref. [95] considered a different approach to Lieb-Schultz-Mattis like constraints in disordered systems. 14

14
A. Setting Consider Hamiltonians, or rather statistical ensembles of Hamiltonians, that preserve exact spin symmetry and statistical translation symmetry. For the purposes of this discussion the spin symmetry will be taken to be SO(3). 15 'Statistical' translation symmetry means that the probability distribution of the disorder, which is assumed to be short-range-correlated, is translationally invariant. This statistical translation symmetry allows us to define a unit cell. Crucially, we assume that there is half-odd-integer spin per unit cell.
It is important that the spin symmetry is exact rather than merely statistical: for example the Hamiltonian H = − i h i · S i , with the random fields h i uniformly distributed on the sphere | h i | = h, has a manifestly trivial ground state.
We consider systems on an L × L × . . . × L torus of even size L. First we conjecture that the averaged energy gap ∆E must vanish with L at least as a power law in this situation. Note that this conjecture is only possible with both statistical translation symmetry and shortrange correlated disorder: without these conditions it is easy to construct counterexamples. For example, we may obtain a gapped Hamiltonian by explicitly dimerizing the couplings in the pattern of a defect-free dimer configuration; but this does not correspond to short-range correlated disorder. 16 We now consider more detailed constraints. To avoid trivial counterexamples to some of the statements, let us restrict to ensembles for which the ground state on the torus is unique with probability one (and therefore a singlet). From now on, let us also restrict to ground states that do not spontaneously break either the spin symmetry or the statistical translational symmetry. We discuss how to define this in the disordered setting below using appropriate local order parameters. Note that in a generic 1d or 2d disordered system translational symmetry breaking is prevented by Imry-Ma.
Lieb-Schultz-Mattis-Hastings-Oshikawa theorems may in certain cases give a gapless spectrum but do not give constraints for its spin/charge quantum numbers; issues related to such quantum number restrictions were explored in Ref. [96]. 15 We expect these statements also to apply for O(2) symmetry i.e.
Sz rotations plus a discrete spin flip that reverses Sz. Similar restrictions should also apply for discrete symmetries, e.g. Z 2 × Z 2 , or time reversal. 16 It is also useful to check that the case of a VBS state in 3d, with weak disorder, is consistent with our conjecture. In 3d the VBS order is stable to weak disorder. In the clean system the reason for the small gap is the existence of multiple symmetry-equivalent VBS states. In the disordered system these states are split by an amount scaling as the square root of the volume. Instead the gap is closed by rare region effects: a small density of patches where the VBS order parameter can be flipped with low energy.

B. Spin gap and correlation functions
First consider the spectrum in 1d and 2d. Strengthening the above conjecture regarding the vanishing of the gap ∆E, we conjecture that in 1d the average gap to spinful excitations, ∆E S , also vanishes at least as fast as a power law at large L. This is natural given the argument in Sec. IV C for power law spin correlations.
In 2d we must allow for the possibility of topological order: a spin liquid, with a spin gap, is of course a possible ground state. However we conjecture that a spin gap is possible only in the presence of topological order. In the absence of topological order we conjecture that ∆E S again vanishes at least as fast as a power law at large L. Therefore, either ∆E S is bounded by L −a for some a > 0, or the gap to all excitations, ∆E, is exponentially small in L as a result of topological ground state degeneracy on the torus. In a disordered system, especially when the disorder distribution is not bounded, the gap may of course close for relatively trivial reasons as a result of Griffiths (rare region) effects. But we may contrast the above with the case of a weakly disordered magnet with an even number of spin-1/2 per unit cell, for example a ladder with dimers on the rungs, or more nontrivially, the featureless gapped states that can arise on e.g. the honeycomb lattice. [97,98] For weak, bounded disorder, these states remain gapped.
Nevertheless it is worth also looking for additional diagnostics in the disordered setting. One possibility is to examine correlation functions. See also Sec. IV D where we discuss entanglement.
In fact, for the 1d chain it is possible to show that if statistical translation symmetry is unbroken (this is guaranteed by Imry-Ma for a generic 1d disordered system) then the disorder-averaged S z S z correlation function decays no faster than 1/(distance) 2 . A version of this follows from a general theorem about classical probability distributions for spins or charges on the line [99], and for ground states of the translationally-invariant (clean) Heisenberg chain this statement is proved in Ref. [100,101]. Related ideas have also been discussed in Ref. [95]. In Section IV C we give a self-contained argument for any singlet states -which applies even if the states are not ground states of local Hamiltonianswhich conforms to the definition of statistical translation breaking given below.
Since the Imry-Ma mechanism prohibits translation breaking in a generic random 1d model, this implies that the random singlet phase, for which | S z i S z j | ∼ |i−j| −2 , exhibits the fastest possible power law for decay of spin correlations.
Given the above, a natural question is whether in 2d, if the system is not a spin liquid, there must again be a spinful local operator O S (x) whose two-point function decays slower than exponentially with distance.
In 3d, we again conjecture that, in the absence of topological order, the spin gap ∆E S must vanish as L → ∞.
(Recall that we are restricting to states that do not break statistical translation symmetry, so for example the VBS state is excluded.) However we speculate that a much weaker scaling is possible than in lower dimensions, with ∆E S tending to zero only logarithmically with L. The mechanism for this is as follows. We have seen that a valence bond glass is unstable in 2d. In 3d magnets, however, it is possible to consider a wider range of paramagnetic states built from singlets. We will discuss these in Ref. [102]. Unlike in 2d, stable glassy states can be constructed. These allow for rare regions in the form of defect loops, with long loops being exponentially rare. These loops carry spinful excitations whose gap is powerlaw small in the loop length. This necessarily closes the spin gap. But if the ground state of each loop is a singlet, the gap closes weakly. The longest loops in a sample of size L are logarithmically large in L, giving a logarithmically small spin gap.
Finally let us consider what it means to break spin or translational symmetry in the disordered system. We will not attempt to be mathematically precise (for example we neglect the possibility that the large L limits below are ill-defined). Spin symmetry is unbroken if for every 'spinful' local 17 operator O S (x), transforming in a nontrivial representation of SO(3), the following two-point function tends to zero for large |x − y|, no matter how the limit of large |x − y| is taken: Here the operator O S (x) depends on x only by simple translations. The angle brackets are the zerotemperature quantum average and the overline is the disorder average. The above rules out for example spin-glass order as well as uniform orders. Translational symmetry is broken if there exists a local operator O(x) (again, depending on x only via simple translations) for which the correlation function has non-decaying oscillations, with a nonzero wavevector, at large distance. 18

C. Restriction on spin correlations in 1d chains
In this section we construct an argument restricting the spin correlations in one-dimensional spin chains that 17 O S (x) is supported on a finite number of sites around x. 18 In principle these oscillations may be present only along certain directions. For example consider a clean 2d system that splits into decoupled 1d chains that are in the VBS phase. The ground state is unique since the ground state of each chain is the superposition of its two VBS configurations. Spatial symmetry is broken but this can only be detected by C(r) if r is parallel to the chains.
have an odd number of spin-1/2 sites per unit cell of statistical translation symmetries. Take a 1D spin-1/2 chain of arbitrary even length L with periodic BCs. We consider an arbitrary probability distribution on states |ψ that are global singlets under U (1) spin rotations around the Z axis, where Z i is a Pauli matrix (σ z ) for the spin on site i. The distribution could be the distribution of ground states of a local random Hamiltonian (in which case the distribution is translationally invariant although the states are not) but it does not need to be. For example we could also consider a distribution which is supported on a single state for each L. In this case the averages below can be dropped. We show that if the averaged Z two-point function, decays sufficiently rapidly with distance, then we can define a strictly local operator V η (x) (supported on a finite number of sites) whose averaged two-point function has oscillations at nonzero wavevector (namely π) which do not decay at large distance, so that statistical translation symmetry is broken according to our definition.
To allow for periodic BCs, let us write (k) L = k + mL, with the integer m chosen so that |k + mL| is minimized, and similarly |k| L = |(k) L |. Then for the above to hold it is sufficient if the correlator is bounded by for L-independent c > 0 and α > 2. The random singlet phase has |C Z (i, j)| ∼ 1/|i−j| 2 , [22] so exhibits the largest α possible for a system that does not break translational symmetry.
The present argument is straightforward using rotations around the Z axis. A similar statement also follows from a general theorem in Ref. [99]. 19 Related ideas are discussed in Ref. [100,101] for clean systems and in Ref. [95], which shows that for a disordered spin Hamiltonian with a 'mobility gap', a certain flux insertion operator (which does not however appear to be a local operator) has oscillations at momentum π. A local twist operator related to the one considered here was used in Ref. [103] for a different construction in a clean system.
Label sites i = 1, . . . L. Letting x ∈ Z + 1/2 label a 16 bond, and fixing an integer η, define the unitary operator is supported on a finite number (2η) of sites (so it is a local operator) and it rotates spins by an angle which jumps from approximately π to −π as the bond x is crossed, and which decays gradually to zero away from bond x. We will show that whenever (13) holds with α > 2, it is possible to choose η large enough so that the averaged two-point function of V shows nondecaying oscillations at wavevector π. In fact for any desired > 0, it is possible to choose η large enough so that, for any sufficiently large L, is within of (−1) y−x . Consider The phase of the rotation in V V † now jumps at both x and y. We can eliminate both of these jumps by rotating the spins in between x and y by an additional 2π. But such a 2π rotation is equivalent to the (−1) y−x which is included explicitly above. That is, where θ j increases from 0 to 2π over a region of length 2η around x, and then decreases back to 0 over a similar region around y. Let us define so that A(0) = 1 and A(1) = (−1) y−x C V (x, y). Differentiating, The RHS is proportional to the inner product of ψ| exp iλ 2 j θ j Z j and k θ k Z k |ψ . The former has norm one and the latter has norm N , where we used j Z j |ψ = 0. Therefore Averaging, The RHS involves It is now easy to check that if C Z obeys (13) with α > 2, the above expression is of order η 2−a for large η. (We examine separately the cases where one, both, or neither of j and k lie within a given region of size 2η where θ is varying.) Therefore |C V (x, y) − (−1) y−x |, and a fortiori |C V (x, y) − (−1) y−x |, can be made as small as desired by choosing η large enough.

D. Entanglement structure of ground states
In this section we give a heuristic picture for the kinds of states that are allowed. (From now on when we refer to 'states' we implicitly assume the absence of long-range order.) We would like to classify ground states according to how 'nonlocal' they are. In a disordered system there are more possibilities than in a clean one.
First, it may or may not be possible to construct the state using a 'finite' depth unitary circuit (in a disordered system this concept must be extended to allow for rare regions 20 ). For example, a ground state which is a product of nearest-neighbour singlets can be constructed from a product of 'up' spins by a quantum circuit of depth one. 21 Let us call states that can be constructed using finite-depth circuits 'constructible'. An example of a state which is not constructible is the ground state of a quantum spin liquid: this requires a quantum circuit of depth proportional to the system size. [104,105] Another example is a random singlet ground state in 1d, which 20 The depth of the circuit would not be strictly finite. Instead there would be a requirement on how the accuracy with which the ground state could be approximated converged with the circuit depth. In a disordered system it is important that such a definition should allow for rare (Griffiths) regions: locally the required circuit depth should be allowed to be arbitrarily large, so long as the probability for this is sufficiently small. 21 Here we are not imposing conditions on the transformation properties of the initial product state, or the unitaries, under spinwe are only constraining the final state. Other protocols (e.g. starting from a product of spin singlets) could be considered.
requires a circuit of O(L) depth to produce the largest singlets. Second, if the state is 'constructible', there is still a distinction to be made according to whether or not the local structure of the circuit depends on the disorder far away. If it does not, we call the state 'constructible with local information'. For a disorder ensemble whose ground state is constructible with local information, there exists a protocol for defining the unitary circuit, given the Hamiltonian, in which the local structure of the unitary circuit depends exponentially weakly on values of the random couplings far away.
We only attempt to make this distinction at a heuristic level. Consider first the example of a 2-leg ladder with nonzero, weakly random AF couplings on the rungs, and all other couplings zero. This has two spin-1/2s per unit cell and we do not have any LSM restriction. The ground state is the product of singlets on the rungs and is evidently constructible with local information: the unitary circuit that constructs this state has no dependence on disorder at all. By contrast, consider a magnet with spin-1/2 per unit cell in high spatial dimension. Assume that in 3d and above, as discussed in Sec. IV B, there exist stable 'glassy' states built only out of short-range singlets (neglecting rare regions). Since these states are close to a product of local singlets, they are also constructible. However, unlike the example of the ladder, they are not constructible with local information. A slight change to the disorder at one location can affect the singlet pattern at a distant location (with a probability that is small, but not exponentially small).
Let us now consider one, two and three spatial dimensions in turn.
In 1d, we propose that there are no constructible states for magnets with spin-1/2 per unit cell -even if we do not require them to be ground states of local Hamiltonians. As we have discussed (Section IV C) it is impossible to write down a singlet wavefunction in which spin correlations are rapidly decaying and in which translational symmetry is unbroken. Therefore there are no constructible ground states for disordered Hamiltonians with statistical translation symmetry. Recall that we are restricting here to states without LRO.
In 1d we also make a conjecture in terms of the amount of entanglement entropy between two halves of a system of size L. We conjecture that the disorder-averaged entanglement entropy necessarily diverges at least logarithmically with L in any ground state without LRO. The random singlet phase [106], and the gapless ground state of the clean antiferromagnetic chain, are examples consistent with this conjecture.
In 2d we conjecture that there are again no constructible ground states for magnets with spin-1/2 per unit cell. In Sec. III A we discussed a putative 'valence bond glass' in 2d. This discussion shows that constructible wavefunctions can be written down (unlike in 1d), but we argued that they could not be ground states of local Hamiltonians, as they are unstable to the nucle-ation of spinful defects.
In 3d and above we suggested the possibility of constructible 'valence bond glass' ground states for spin-1/2 magnets. However, these states are not constructible with local information. This distinguishes them from more 'trivial' ground states that are constructible with local information. We propose that the latter can only exist for magnets with integer spin in the unit cell.

V. APPLICATION TO EXPERIMENTS
Let us now apply these theoretical ideas to an experimental setting, focusing on the particular case of YbMgGaO 4 . A schematic partial summary is offered in Table II.
Before comparing with experimental data, we note two differences from the theoretical scenario described above. First, in YbMgGaO 4 bond randomness is likely not weak, since the differing charges on Mg 2+ and Ga 3+ modify the oxygen charge-transfer energy as well as the Yb-O-Yb bond angle [9]; in iridium oxides with equivalent edgesharing oxygen octahedra such variations are known [107,108] to impact the SOC-based magnetic exchanges. Of course since the material is intrinsically disordered, the energy scales ∆ S and ∆ V B , which are defined only in the theoretical clean limit with vanishing disorder strength ∆ → 0, cannot be accessed directly.
Second, we note that strong SOC breaks spin SU(2) rotation symmetries; 22 however VBS phases remain well defined as quantum paramagnets that preserve timereversal symmetry but break some lattice symmetries. It is also conceptually helpful to consider two-spin singlet states, which remain eigenstates even with SOC as long as J µν ij is symmetric, as may be enforced by inversion symmetry across a bond midpoint, since the singlet state is odd under inversion (see Appendix 6 for details).More generally, time reversal symmetry and Kramer's degeneracy is enough to protect the S=1/2 Kramer's doublet at the core of an isolated defect, even without any spin rotation symmetry.
To discuss the consequences of the pinned singlets and their instabilities, let us consider in turn two distinct regimes of temperature and energy. 22 Spin-orbit coupling can produce multiple spatially-varying anisotropic exchanges. In YbMgGaO 4 , crystal symmetries restrict nearest-neighbor interactions J µν ij to four independent parameters; these are sometimes denoted by J z , J xy , J z± , J ±± in the literature [2], but may be fruitfully rewritten in terms of the Kitaev-type exchanges that occur in iridates with analogous edge-sharing oxygen octahedra, namely J z , J xy , Kitaev K, and pseudo-dipole I interactions; see Appendix 8 and also Ref . 32.

S[H; q, ω]
Eq. 29, Fig. 8 TABLE II. Shorthand summary of theoretical results for a few experimental observables within the framework of a d > 1 lattice-emergent random-singlet regime following the ideas discussed above for a random frustrated Heisenberg model Eq. 1. These predictions, namely heat capacity, spin susceptibility, thermal conductivity, scaling of heat capacity in magnetic field, and in-field dynamical structure factor, are expected to be seen in YbMgGaO4 and related systems at low temperatures.

A. Structure factor at low energies
At low or intermediate temperatures and energies of order J, inelastic neutron scattering data [5][6][7][8] suggests that the random singlets scenario is on the right track. Two distinctive features appear at the lowest frequencies ω 0.1 meV in low temperature: (a) increasing intensity from the BZ center to the BZ edge, and (b) a broad maximum at the BZ edge midpoint, position M . The features persist across the energies 0.1 < E (meV) < 1; for comparison the magnetic exchange scale has been estimated[1-3] to be roughly 0.2 meV. Randomly oriented short-ranged singlets capture both features, as we show in Fig. 5. Note that feature (b) clearly implies some second-neighbor correlations; the wavevector M is oriented towards second neighbors on the triangular lat-

FIG. 5.
Spin structure factor for pinned singlets. Equal time spin correlations S + S − for spins frozen into randomly-oriented short-ranged valence bonds show a deep minimum at the Brillouin zone center; a 4:1 ratio of first/second-neighbor bonds also give a broad maximum at the edge midpoint M . Both features are seen in neutron scattering [5][6][7][8].
tice. Since the distribution of short ranged S = 0 resonances can vary even within a particular VBS phase through different resonances, we take the ratio of first-to second-neighbor singlets as a free fitting parameter. In this short-ranged-singlets phenomenological model, the equal time spin correlations S + S − (for dynamics see Eq. 29 below) are easily computed and are given by for a fraction f 1 (f 2 ) of spins participating in first (second) neighbor singlets; a i (b i ) are first (second) neighbor vectors. A good fit is found for a ratio f 1 :f 2 = 4:1, as used (with f 1 +f 2 =1) for Fig. 5. The experimental observations at T J are well reproduced by random short-ranged-singlets.

B. Heat capacity at ultralow energies
At ultralow energies and temperatures below J, the magnetic specific heat C(T ) ∼ T 0.7 provides direct evidence for a power-law density of states of strong randomness. To compare specific heat as well as magnetic susceptibility across a wider temperature range, we introduce a two-regime model which, despite its simplicity, can still capture the experimental data semi-quantitatively as seen in Fig. 6. The model is defined by considering a density of states consisting of two distinct regimes, whose integrated density of states is plotted in Fig. 6(a). For this model we let a fraction f of the sites form singlet pairs at energies tightly distributed (in a triangular distribution) within the narrow range (J 0 −∆, J 0 +∆), while at lower energies the remaining 1−f sites form either singlet pairs, or alternatively coexisting singlets and largerspin triplets, in either case with energies drawn from the power-law distribution ρ(E) ∼ E α−1 with α = 0.7. (The exponent α is non-universal in the theory; for example, a phenomenological fit for YbZnGaO 4 would instead entail α = 0.59.) The effective average magnetic energy scale 23 J 0 may be estimated [1][2][3] as J 0 ∼ 1.5K. The results of the model are largely insensitive to the value of ∆; here we took ∆ = 0.15J 0 to match the g-factor broadening extracted from crystalline electric field excitations [9]. By construction, the model captures the observed magnitude and scaling of specific heat, C = 1.41 T 0.7 J/mol K = 0.17 T 0.7 k B with temperature measured in Kelvin. Within a zeroth-order approximation of frozen isolated singlets (see Appendix 7 for details), the computed specific heat captures the experiments at the parameter value f =4/5; a better approximation including (model-dependent) singlet fluctuations would likely increase C(T ) and require a larger f for fitting the experimental C(T ) data. When the renormalized spin coupling is antiferromagnetic such that only singlet pairs form at low temperature, the susceptibility χ(T ) is tied to C(T )/T and shows a soft divergence χ ∼ T α−1 at ultralow temperatures below a pronounced peak (Fig. 6 (c)). Admixing a fraction of large-spin clusters, likely to develop from any ferromagnetic bonds in the low-energy distribution, will contribute an additional 1/T Curie term, flattening out the peak to produce an apparent saturation before the rise (Fig. 6 (b)). Experimentally an apparent saturation is seen [1,3] around T = 0.3 K. The lowest temperature measurements of susceptibility are thus consistent with the random singlet mechanism; the random singlet phase would predict a slow rise of susceptibility at very low temperatures. 24 At even lower temperatures, if large-spin clusters develop and freeze into a spin glass -as is the likely T =0 RG fixed point fate of the random spin network -then the resulting glassy moment formation would serve as a cutoff for an otherwise-divergent magnetic susceptibility, and simultaneously would be visible for standard probes of spin glasses.

C. Thermal conductivity
Next let us consider the ultralow temperature thermal conductivity. Since the gapless spin excitations are all localized, heat is carried primarily by phonons. At ultralow temperatures the distribution of random singlets (as well as any frozen spin-glass moments, if those arise at lower temperatures) both produce a collection of quantum twolevel systems, associated with local rearrangements of the singlets and their levels or spin configurations, which 24 Ideally the q = 0 susceptibility at ultralow temperatures should be measured directly; measurements of local suscepibility, through the µSR Knight shift, may be further complicated by interactions of the muons with the charged Mg/Ga disorder.
scatter acoustic phonons via resonant absorption at a rate Γ that is linear in the phonon frequency [109,110]. The leading contribution to low temperature thermal conductivity κ = Cv 2 /3Γ (with v the acoustic velocity) arises from this phonon scattering mechanism, leading to the anomalous scaling κ ∼ T 2 . This ultralow-T thermal conductivity, in particular the vanishing κ/T → 0 as T → 0, presents a contrast with the finite κ/T at zero temperature that would be expected from a spinon-Fermi-surface with itinerant spinons. Experimentally [4] κ/T was found to vanish with T → 0. The κ ∼ T 2 power law naively expected here is a familiar theoretical result for structural as well as spin glasses, where log corrections and other effects are known [110] 25 to result in experimentally observed effective power laws κ ∼ T α with a slightly reduced exponent 1.8 α 2. Thermal conductivity measurements [4] in YbMgGaO 4 found a power law behavior κ ∼ T α with a α = 1.85, and measurements [12] in YbZnGaO 4 found κ ∼ T α with a α = 1.97, both in excellent agreement with this full theoretical description.

D. Signatures in a magnetic field
Applying an external magnetic field H will produce signatures in both low-temperature thermal transport as well as heat capacity. We will discuss the effects of an external magnetic field H on three physical observables: low-temperature thermal transport; heat capacity; and the dynamical structure factor.
First we consider how the form κ ∼ T 2 for thermal conductivity is modified by a magnetic field. When the Zeeman energy grows larger than the singlet gap for a given set of singlets or frozen moments, their ground state becomes spin polarized and no longer participates in scattering acoustic phonons. This reduction in scattering mechanism will be directly observable as an en- Specific heat in a magnetic field. Specific heat C(T ) is computed within the two-regime model with low energy singlets, as in Fig. 6(c,d), here with an application of a magnetic field H with dimensionless energy h ≡ gµBH/J0 taking the values h = 0, 0.1, 0.2, 0.5, 1, 1.5, 2, 3 (red to violet colors). The anomalous curve at h = 1 is due to a resonance with the narrow distribution chosen for Fig. 6(a). Applying a small field changes the low temperature scaling of specific heat from C(T ) ∼ T 0.7 to a linear form ∼ T /H 0.3 , then turning to an exponential form at larger fields. (For modifications of the scaling by spin-orbit coupling see also Ref. [111]). For fields h < 1 smaller than the singlet-triplet splitting, low intensity features are visible at the Γ point (q = 0): fixing frequency at ω = ω± ≡ 1 ± h (here ω± = 0.5, 1.5) and taking q → 0 shows the intensity decreasing. The spread of intensity in q and ω away from these points can manifest as "constriction points" of minimum intensity. (b): For fields h > 1 that start to polarize a given singlet, a q = 0 peak is observed at a frequency ω0 ≡ h, which scales linearly with field. Coupling the isolated singlets will connect the ω = h, h − 1 features into a magnon dispersion. (c): Given some broad distribution of energy splittings, e.g. as in Fig. 6, the structure factor at q = 0 shows a single peak, at ω = h, with an intensity that increases with h and saturates for h 1. This q = 0 feature, as well as the q → 0 decrease in intensity approaching the "constriction points", are observed in experiments [10,11]. hancement of thermal conductivity in a magnetic field. For large enough fields of order a few times the average exchange J, the spins will all be essentially polarized and κ will saturate at a value substantially larger than its zero-field form. This enhancement and saturation is indeed seen [4] in YbMgGaO 4 at fields of a few tesla.
Second we consider effects of magnetic field on the heat capacity and magnetic susceptibility. The behavior of specific heat in a magnetic field, computed within the two-regime model discussed above, for the case of low energy singlet formation, is plotted in Fig. 7. Within the low-temperature power law tail, applying a small magnetic field H changes the scaling of C(T ) from T 0.7 to linear in T , following the scaling function which holds for a power law distribution of pure spinsinglets. (A modification of this scaling function under particular forms of spin-orbit coupling will be discussed elsewhere [111]).
Scaling in a magnetic field is also seen in other observables. Susceptibility χ[H, T ] is cut off to scale as a constant proportional to H −0.3 at low T . The scaling limit requires H and T to both be sufficiently smaller than the lattice energy scales, so as to have a chance of being in the emergent power law regime. Fields with Zeeman energies larger than the lattice exchange energy scales J 0 polarize the spins, leading to a gapped spin polarized state with asymptotically exponential form (C(T ) ≈ (h − 1) 2 /2T 2 exp[(1 − h)/T ] with T in units of J 0 , with dimensionless Zeeman energy h ≡ gµ B H/J 0 ). Adding triplets or larger spin clusters into the low energy spectrum would modify the susceptibility sharply and specific heat somewhat more mildly. In YbMgGaO 4 , a scaling regime may be difficult to reach, experimental studies of C(T ) in fields above 1 tesla have been reported [1] and are consistent with this sharp decrease of C(T ) in a field.
Third we briefly consider the dynamical spin structure factor in a field. A spin singlet with some effective splitting J 0 will transition from the singlet ground state into a polarized state when the Zeeman energy overcomes the splitting J 0 , at H = J 0 (here H is the Zeeman energy, i.e. the magnetic moment is set to unity). The spin dynamics S + S − transverse to the applied field, for a valence bond involving two sites separated by the vector R, are straightforwardly computed by considering the dynamics in the singlet or spin-polarized ground state of the two-spin Hilbert space. This yields the structure factor contribution of a singlet bond with energy and separation {J 0 , R}. The expression for the structure factor, which holds for finite field H > 0 as well as zero field H = 0, is with θ[x] the Heaviside step function, assuming full spin rotation symmetry. In the presence of spin-orbit coupling, the H = J transition is smeared and full polarization is achieved only asymptotically at large fields. The signatures in the spin dynamics are shown in Fig. 8. The distribution of singlets used for Fig. 8, namely 4:1 first:second neighbor singlet bonds, is the same as that extracted from the zero-field structure factor as shown in Fig. 5. At wavevectors near the Brillouin zone center, q = 0, the dynamical response shows a constriction at low fields into a feature of minimum intensity ( Fig. 8 (a)). Sitting at ω fixed to the frequency of this "constriction point", and taking q from the BZ edge to q = 0, shows the intensity decreasing. Interestingly within a spinon Fermi surface [10,16] interpretation the expectation for this spectral-crossing-like "constriction point" is for an opposite feature, namely a sharp increase in intensity as q → 0. Experimentally, the intensity approaching the "constriction point" or crossinglike feature is observed to decrease [10]. At large enough fields ( Fig. 8 (b)) the dynamics show a delta-function peak at a frequency proportional to the magnetic field, S(q=0, ω) = δ(ω − h) for h > J. The distribution of effective J implies that this field-gapped delta function response, S(q=0, ω) = δ(ω − h), will be observed at any magnetic field, but with an intensity that grows with field and then saturates above h ≈ J 0 (Fig. 8 (c)), again in apparent agreement with recent experiments [10,11]. At very low frequencies, the expression S[J 0 , R](q, ω) of Eq. 29 must be integrated against the density of states ρ[E] to yield the power-law form of the low frequency structure factor, where R(E) is the singlet size as a function of its energy splitting, which can depend on the lattice and other details; in the 1d spin chain random singlet fixed point [22], log E ∼ √ R.

E. Other experimental signatures
Let us also consider signatures of this theoretical scenario in other measurements. Inelastic neutron scattering at ultralow frequencies would distinguish the present picture from a spinon-Fermi-surface, expected to show some 2k F structure, as well as from any magnetically ordered domains, which would have magnetic Bragg peaks visible in scattering and ordered moments visible in e.g. µSR. Short ranged VBS order may be observed, through high-order crystalline symmetry breaking effects which can be seen in sensitive structural probes; for columnar VBS, the lattice modulation would be visible at wavevector M , though we note that the (non-infinitesimal) microscopic randomness present in YbMgGaO 4 means that short-ranged VBS order is not required by the theory. NMR spin-lattice relaxation is expected to exhibit a stretched exponential associated with a distribution of 1/T 1 decay coefficients. Finally, though the ultimate fate of the RG flow for the nucleated defect spins is unknown, on the triangular lattice as discussed above one reasonable expectation is that large ferromagnetic spin clusters form, and then freeze into a spin glass; signatures of spin freezing, with several unusual properties, were recently observed [12] in both YbMgGaO 4 and YbZnGaO 4 .

F. Relevance to other materials
The simplest starting point for the theoretical treatment presented above is a clean Hamiltonian in a VBS phase. Such a setting is well approximated experimentally by the organic material EtMe3P[Pd(dmit)2]2, a triangular lattice S = 1/2 magnet with an antiferromagnetic Heisenberg coupling of order J=250 K, which has been shown [43] to exhibit valence bond solid order below an ordering temperature T =25 K. Given such a clean system in a VBS phase, one can hope to introduce quenched disorder in a controlled manner. Indeed irradiation by x-rays has been shown [112] to produce nonmagnetic disorder that can nevertheless substantially modify the magnetic properties of organic compounds. For EtMe3P[Pd(dmit)2]2 or similar systems, controlling the irradiation time would allow experimentalists to interpolate between the well-understood theoretical limit of weak disorder and the relatively strong disorder seen in YbMgGaO 4 and YbZnGaO 4 .
The physics discussed here is expected to be quite general. For very weak disorder, the very small fraction of spins that participate in the emergent power laws may make detection difficult even as the materials tend closer to the theoretically controlled limits; and in other compounds, the presence of site disorder in addition to bond randomness, such as dilution or random fields, gives an obvious source of magnetic impurities that does not require the analysis provided here. But in many materials time-reversal-symmetry is well preserved microscopically even by disorder, and there are numerous possible microscopic sources of spatial variations in the energy, which corresponds to bond randomness, that then permit an analysis based on the theory provided here.
Indeed since this work was first released a variety of other magnetic insulators have been found which appear to be well described by the present theory. The essence is a heat capacity that exhibits power law scaling with a fractional exponent, with a magnitude consistent with a contribution arising from a portion of the local moments in the material. At the time of writing these compounds already include the layered spin-1/2 Mott insulators H 3 LiIr 2 O 6 ,[113] LiZn 2 Mo 3 O 8 , [111,[114][115][116] synthetic herbertsmithite [34,111,117] and 1T-TaS 2 . [118] All of these compounds show power law heat capacity, and moreover a data collapse of heat capacity in a magnetic field, into a scaling function of the single variable T /H, analogous to Eq. 28. Indeed 1T-TaS 2 appears [118,119] to exhibit exactly the C ∼ T /H 1−α scaling function of Eq. 28. The other three materials exhibit a similar scaling function but with an additional factor of T /H at low temperature; the complete analysis of this modified scaling, based on a particular role of spin-orbit coupling, will be discussed elsewhere [111].

VI. DISCUSSION
Quenched randomness allows for new kinds of quantum ground state with interesting entanglement structures, and destroys states that would otherwise be stable. One of the main results of this paper is a study of the path by which weak disorder transforms a paramagnetic valencebond solid into a state with spinful excitations, starting with the nucleation of spin-1/2 vortex defects. Spin-1/2 defects were also inevitable in a stronger disorder regime where the starting point was a 'glassy' covering of short-range valence bonds rather than an ordered one. Together, both of these results motivated our disordered-LSM conjectures.
The starting point for this analysis requires us to make the distinction between a state of frozen short ranged singlets, which we denote a valence bond glass, and valence bond states with gapless spin excitations such as random singlet phases with singlets at arbitrarily long length scales and low energies. We exploited a weak disorder limit in which the renormalization group flow could be partly tracked. In this limit, the relevant energy scales are exponentially small, likely precluding any observation of this regime in experiments. Here we chose to focus on experimental applications to YbMgGaO 4 , which shows complex phenomenology that is captured by the theory for stronger disorder. Indeed disorder arises in YbMgGaO 4 , as well as the isostructural and phenomenologically analogous compound YbZnGaO 4 , through the nonmagnetic layers which in both cases appear to be fully amorphous. We expect the theory to also describe other strongly frustrated compounds with much milder bond randomness; in those cases, the small emergent energy scales may require further care to observe.
We have focused here on systems in one and two spatial dimensions. In 3d a VBS is stable to weak disorder, but similar issues could arise in the presence of sufficiently strong disorder. Distinct types of paramagnetic states are also possible in 3d: we will discuss these separately [102]. In the context of 3d systems with random valence bond patterns it is interesting to consider the 3d magnet Ba 2 YMoO 6 , where degenerate orbital and spin-1/2 moments residing on a face-centered cubic lattice are seen to freeze into a disordered pattern of orbital dimers and spin singlets. [120] Defect spin-1/2 moments are observed at finite density within the dimerization pattern, and these defects appear to freeze into a dilute spin glass at ultralow temperatures T g ∼ 0.01θ CW with θ CW =160 K the Curie-Weiss temperature. [121] A recent theory [122] suggests that this dimer-singlet phase can arise from an unusually robust extensive degeneracy due to the orbital dimers; these properties of the orbital moments, together with a lack of observed bond randomness but apparent presence of some magnetic site dilution disorder [120,123], suggest that the physics in Ba 2 YMoO 6 may involve additional ingredients beyond those considered here.
Returning to 2d magnets, it is worth noting that even for a layered material with vanishing magnetic coupling between two layers, some three-dimensionality will be required in describing the VBS order in the real material, due to phonons. The 3d phonons couple to the lattice modulations of the VBS order, thereby giving an effective coupling between VBS order parameters across different layers. This gives a three dimensional VBS order in the clean system, which then melts only under finite, but very weak, disorder.
Some of the phenomena discussed here could be studied numerically. Sign-free Quantum Monte Carlo can be used to numerically study the fate of the weakly disordered VBS on the unfrustrated square lattice [37,62]. However the very small energy scale of defect interactions -in the theoretically controlled ultraweak disorder limit of Sec. II, it is exponentially small in ξ 2d -could be an obstacle to numerical studies of the low energy regime. 26 (In principle, an indirect 'multi-scale' approach might be possible, with separate simulations to establish the geometry of the defect array and to treat their interactions.) Further studies of the phase diagrams for experimentally relevant Hamiltonians, even in the clean limit, would also be useful as a starting point for understanding the disordered materials. We note that when bond randomness suppresses a given magnetic ordering without immediately producing a spin glass, the resulting low energy excitations may be well described by the RG flow from pinned singlets discussed here.
A particularly interesting aspect of our considerations is the restriction on ground states of disordered quantum magnets with spin-1/2 per unit cell and with statistical lattice symmetries. The restriction naturally explains the emergence of low energy excitations triggered by disordering the VBS symmetry breaking order. We hope that the physical arguments presented in this paper lead to mathematically rigorous scrutiny of such Lieb-Schultz-Mattis like restrictions in disordered systems. The specific conjectures we formulated should be a good target for such future studies.
Acknowledgements. We thank Leon Balents,  We consider the effective interactions J eff between the spin-1/2s in the cores of defects, in either 1d or 2d. We restrict to the limit of weak disorder ∆, where the typical defect separation ξ diverges, and where J eff is exponentially small in ξ. In particular we are interested in whether the signs of these effective couplings are necessarily randomized due to disorder in this ∆ → 0, ξ → ∞ limit or whether they can retain a definite sign structure.
For simplicity, let us assume that the vortex spins are well-localized at a single site. In the large ξ limit, interactions involving more than two vortex spins are strongly suppressed relative to the two-spin interactions, so it suffices to focus on a pair of vortices in an infinite system, with Pauli operators σ 1 and σ 2 . A convenient way to define their effective Hamiltonian is The resulting H eff will be approximately β-independent for sufficiently large β (apart from a trivial constant term) and by SO(3) symmetry will be of the form The effective coupling J eff is determined by the exponentially decaying spin correlations of the 'other spins' in the background configuration of pinned VBS domains. This can be seen in linear response, where we treat the coupling between each defect spin and its neighbours as a small parameter (removing this approximation would only dress the operators appearing in the correlator below, and would not change the basic behaviour). Let us write the couplings in H that involve σ 1 as 1 σ 1 . S 1 and similarly for σ 2 . Here S 1 is a weighted sum of the magnetizations of the sites surrounding the defects (the coefficients depending on the microscopic Hamiltonian) and the coupling strength 1 is treated as a small param-eter. We have From (31), expanding in gives 27 Here it must be remembered that the expectation value is taken not in the clean system, but in a configuration in which VBS domains of scale ξ are pinned in a specific pattern that depends on the disorder realization.
In Sec. 1 a we consider a class of 'sign-free' models on bipartite lattices, where J eff can be shown directly to have a non-random sublattice sign structure, in any disorder realization. While this class includes various important models, sign-freeness is a fine-tuned property which is not present in generic models. Therefore we must ask whether the sign structure of J eff can survive (in the limit of interest where disorder is present but small, and ξ is parametrically large) when the sign-free property is broken. In Sec. 1 b (which is independent of Sec. 1 a) we argue that the sign structure found for the sign-free models can survive in generic models in the relevant limit ∆ → 0, ξ → ∞.
Our discussion of generic models uses the fact that correlation functions of gapped degrees of freedom can be understood, very generally, in terms of directed path ensembles. [124] We also point out that generic models exist where the sign of J eff is completely randomized at weak disorder.

a. Sign-free models on bipartite lattices
The sign-free class we consider includes, among many other models, Heisenberg models on bipartite lattices with (perhaps random) nearest-neighbour AF couplings. 28 The class also includes the (randomized) 'JQ model', [37] which includes a 4-spin resonance term which can be used to drive the clean system into a VBS phase. For these models J eff has a definite sublattice sign structure: sign J eff = −χ, where we define χ to be +1 if the two defects are on the same sublattice and −1 if they are on opposite sublattices. This sign structure is well 24 known in the context of strong-disorder RG for the 1d AF Heisenberg chain and generalizations. [19,62] For various models it also follows directly from Eq. 33 and standard 'loop' representations (see e.g. Ref. [125]) of correlation functions. However, let us give a simple general argument which holds for arbitrary disorder strength and dimensionality and does not require the loop representation.
The defining property of this class of sign-free models is that all off-diagonal H elements are negative when the Hamiltonian is written using the single-site basis |α , α = 1, 2, defined by (so that a singlet between opposite sublattice sites is . This means that Trotterizing Tr e −βH in this basis maps the partition function to a classical statistical mechanics problem, with positive weights, for degrees of freedom α = 1, 2 in spacetime. By Eq. 31, the partition function of interest to us, which yields matrix elements of H eff , is Tr other spins e −βH α β ,αβ . This has periodic BCs in time at all sites except for the defect sites, where there are temporal 'boundaries' at t = 0, β where the 'colour index' is fixed by the matrix element considered. Consider two defects, one on the A and one on the B sublattice, and label their states by α and β respectively. There is only one nontrivial term in H eff , which can be written in terms of the singlet projector P, so we can always write (for some constants A, B) e −βH eff α β ,αβ = Tr other spins e −βH α β ,αβ = Aδ α α δ β β + Bδ α β δ αβ .
The first term is the identity and the second term is proportional to P. If we instead consider defects on the same sublattice, with states labelled by α 1 , α 2 , then the singlet projector may be written P = (1 − S)/2, where S is the swap operator, which in components is the final term below: If we take β β min and βJ eff 1, then B/A is proportional to βJ eff in the first case and −βJ eff in the second case. So to show that the coupling is fixed by the sublattice (AF for opposite sublattices and F for same sublatice) it is sufficient to show that the constants A and B are always positive for the sign-free models considered.
The LHS of each formula maps to a 'classical' partition function with fixed boundary conditions at t = 0, β for the defect sites, and is manifestly positive. This establishes the claim about the sign structure of J eff for this class of sign free models.

b. Generic models
It is clear from Eq. 33 that models exist where sign J eff is completely randomized in the weak disorder limit. As mentioned in the text, this will certainly occur in 1d whenever the exponentially decaying spin correlations in the clean VBS have incommensurate sign structure, as a result of the large random separations between defects. For a 2d model where the spin correlations in the clean VBS are incommensurate, randomized sign J eff is also the natural possibility (although it does not strictly follow because of a caveat discussed at the end of this section).
It is less obvious whether the non-random sign structure for certain bipartite models described in Sec. 1 a can survive in models that are not fine-tuned (in the limit of interest, where ∆ → 0 and ξ → 0 in the manner prescribed by Imry-Ma). First recall that the exponentially decaying correlations of gapped degrees of freedom (here the spin) can be understood in terms of directed path expansions. [124] A simple classical example is the Ising model in the disordered phase, where the hightemperature expansion relates the two-point function to an effective partition function for a path which connects the two points. This formalism is essential for understanding the effect of disorder on correlation functions of disordered degrees of freedom. Disorder is an RGrelevant perturbation for sums over paths of this type.
Similar path expansions may be written for the quantum magnet. Heuristically, these paths are tunneling trajectories for gapped spin-1 excitations. The details will not concern us, but it is useful to have in mind a toy model. Let us take the valence bond covering to be determined by strong explicit dimerization in the Hamiltonian. This is not of course the regime of interest to us, but it illustrates the basic features and makes the path expansion simple. Fix a dimer covering that is complete except for the two defects, and take the couplings on dimerized bonds to be of order J strong and couplings on undimerized bonds to be of order J weak with J weak /J strong 1. In the limit of small J weak , the perturbative expression for J eff is dominated by directed paths which traverse the minimal number k of weak bonds, and is of order J weak (J weak /J strong ) k−1 . More precisely, in this limit where we sum over paths of minimal k. In the numerator the product is over all weak bonds traversed by a given path. In the denominator the product is over all the strong bonds that the path visits: the path can visit a strong bond either by traversing it, or by taking two successive steps on weak bonds that are adjacent to the strong bond. 'Length' is the total number of bonds traversed by the path. The simplest case is 1d, where the above sum reduces to a single term. On a bipartite lattice in any number of dimensions, every term in the above sum has the same sign (if the couplings are all positive). This is because on a bipartite lattice all paths between two points have the same length modulo two. In this case J eff is given by a sum of positive terms, once we extract the sign factor sign J eff = −χ which depends only on the sublattices of the two defects. After extracting this sign factor, the expression for J eff above defines an effective classical partition function for a path or 'string' whose energy depends on which links it visits.
A similar picture, in terms of a partition function for a directed path, will be valid at a coarse-grained level even away from the artificial limit above, so long as spin correlations are exponentially decaying. [124] For a signfree model, the Boltzmann weights for the effective classical partition function will always be positive. The exponential decay of J eff with separation is equivalent to the nonzero free energy per unit length, or 'line tension', of this path. In a clean 2d system, this line tension ensures that the sum is dominated by trajectories that are straight on scales of order ξ 2d . 29 In a given disordered VBS background, the geometry of the path may not be straight on this scale, but a section of the path away from VBS domain walls will be approximately straight. 30 By standard results on paths in random media, such a section will deviate from being straight on a parametrically smaller scale ∼ ξ 2/3 2d (neglecting logs 31 ) but this will not concern us. [126][127][128] In a model that is not sign-free some of the paths that are summed over in the effective classical partition function determining J eff will acquire negative 'Boltzmann weights'. The question is whether J eff retains a fixed sign when the breaking of sign-freeness is weak. There are two distinct ways to 'weakly' break the sign-free condition on the off-diagonal matrix elements of H: (1) First, we can introduce weak sign-rule violating couplings everywhere. For example, in our toy model we may introduce AF second-neighbour couplings with a small magnitude J 2 J weak . Such weak violations of the sign rule will be present generically even in the clean system so we must consider them. 29 If the VBS breaks rotational symmetry the line tension can depend on orientation. 30 Rotational symmetry can be broken in different ways in different VBS domains, so that the orientation-dependence of the line tension differs in different domains. The path is of length ∼ ξ 2d , so it can visit more than one (but order one) VBS domains. In a given VBS background its energy may be minimized by a nonstraight trajectory crossing through more than one domain. The path could also prefer to be 'glued' to a VBS domain wall, as discussed below. 31 Logarithms in ξ 2d arise because the limit of interest is where disorder ∆ → 0 and ξ 2d → ∞ with ξ 2d ∼ e c/∆ 2 , i.e. ∆ ∼ (ln ξ 2d ) −1/2 .
(2) Second, we can introduce rare locations where the sign rule is strongly violated. For concreteness, we can imagine flipping the sign of J weak in our toy model for a small fraction p of bonds. This is a spatially random effect associated with bonds where the disorder |∆J| is locally strong. More generally, in the following we may take p to be the probability of having a sign-rule violating bond whose strength is above some order one threshold.
We take the weak disorder limit by rescaling the probability density for ∆J on each bond, so that it has standard deviation ∆.
That is we take P (∆J) = ∆ −1 f (∆J/∆), where the distribution f has standard deviation 1. If the distribution f has bounded support, then effect (2) cannot occur in the weak disorder limit of interest to us: the probability p above is strictly zero for small ∆. However, if the distribution f does not have bounded support, then p(∆) tends to zero in some way as ∆ tends to zero. This gives a small but nonzero density of negative bonds in the directed path partition function. In this case we must check whether this small density of negative bonds can affect the very long paths which are relevant in the weak disorder limit.
In 1d, for any distribution with a finite variance, p(∆) scales to zero faster than ∆ 2 . This means that rare strong bonds are not seen on the Imry-Ma scale ξ 1d ∼ const./∆ 2 and effect (2) does not play a role.
In 2d, the much larger Imry-Ma lengthscale means that, depending on the disorder distribution, there may be many negative bonds in a patch of size ξ 2d ∼ e −const./∆ 2 . If the tail of f (u) decays as e −const.u 2 with a sufficiently large constant, this is narrow enough to banish negative bonds on this scale. This is certainly sufficient to ensure that effect (2) does not play a role. This is not in fact necessary -a weaker condition on f would certainly also be sufficient as we discuss below. However, it may be natural to impose this stronger condition on the distribution f anyway. If a patch of size ξ 2d contains a significant number of rare bonds with large |∆J|, there will likely be other disorder effects that make the signs of J eff moot. For example, a strong ferromagnetic bond can nucleate a spin-1, and when these spin-1s have larger density than the vortex defects the large-scale physics is no longer dominated by the latter.
Nevertheless, let us briefly consider the possible effect of (2). The effect of rare negative bonds (with density p 1) on directed path partition functions has been studied extensively, with motivations from various contexts. [129][130][131][132][133][134][135][136] In 1d the sign of the partition function is of course randomized whenever the length ξ of the path is much greater than the crossover scale ξ * = 1/p, i.e. whenever the path encounters a large number of negative bonds. The 2d case is considerably more subtle. The sign is again randomized at large scales whenever p > 0. The crossover scale ξ * depends on what other disorder is present, in addition to the negative bonds. The scaling is ξ * ∼ 1/p if, in addition to the negative bonds, there is disorder of O(1) strength in the positive bonds. [135,136] The case where the negative bonds are the only source of disorder has recently been addressed in Ref. [136], and an exponentially larger lengthscale, p −const.p −1/3 , was found. Here we are interested in yet another case, where there is disorder of parametrically small strength ∆ on the positive bonds, in addition to the parametrically small density p(∆) of negative bonds. By the rare-region logic of Ref. [136], it is clear that the crossover scale ξ * [which we must compare with the path length ξ 2d (∆)] will be exponentially large at small ∆, with the functional form depending on the breadth of the ∆J distribution. This will at least ensure that effect (2) is banished for distributions with tails weaker than than f (u) ∼ e −const.|u| x , where x is some constant smaller than 2. So effect (2) is certainly avoided even for distributions somewhat broader than a Gaussian. A more careful analysis would be required to establish whether this effect can ever play a role for broader distributions (with finite variance). 32 It remains to check that mechanism (1) does not affect the sign of J eff if the sign-rule-violating couplings are sufficiently weak. This is in fact straightforward to see in the path picture. Microscopically, a wrong-sign J 2 S i . S i+2 coupling in 1d, say, will allow the path to take steps of length 2 that incur a negative Boltzmann weight. However if the weight associated with these wrong-sign steps is sufficiently small, then after some slight coarse-graining the effective Boltzmann weight will be positive for all coarse-grained steps. (For example, the negative weight associated with the step of length 2 allowed by J 2 will be outweighed by the positive-weight trajectory between the same points given by two steps of length 1.) This is true also in higher dimensions, and is not affected by weak disorder, which modulates the local coarse-grained Boltzmann weights by an amount of order ∆ but does not change their sign if they are initially always positive and of order 1.
Above we noted that if the clean system had incommensurate spin correlations then random sign J eff is the natural expectation. This is clear in 1d but we note a caveat in 2d. Recall that in a given disorder realisation J eff is mediated by a coarse-grained 'optimal' path P between the defects which has a definite geometry on scales of order ξ. Note that in 2d the defects are also connected by VBS domain walls. The line tension of the optimal 32 A naive estimate of ξ * suggests that effect (2) may be absent for any distribution f with finite variance, but this estimate may be too crude. The initial disorder ∆ is small but grows under the RG. It will cause a crossover from random walk behaviour to disorder-dominated behaviour at a lengthscale ∆ −4 measured longitudinally, or equivalently a lengthscale ∆ −2 measured transverse to the path. [137] Naively one might think that in the present case this transverse lengthscale ∆ −2 should play the role of * in the rare-region argument of Ref. [136]. This would give ξ * ∼ p(∆) −const.∆ −2 , which is sufficiently large to ensure that ξ * ξ 2d , so that effect (2) cannot play a role at all. However this estimate may be invalid because for a sufficiently broad power law distribution it gives a ξ * that exceeds the previous result for the case ∆ = 0, p = 0. path P depends on the local VBS environment. Therefore it could be the case that the effective free energy of the path P is minimized when it is 'glued' to a VBS domain wall. That is, spin correlations could be mediated principally by one of the domain walls connecting the defects. Further, in this situation it is conceivable (although it sounds rather contrived) that the correlations mediated by such a domain wall could have a commensurate sign structure even when the correlations in the bulk of a VBS domain are incommensurate.
Finally we comment on the triangular lattice case. Recall that for columnar VBS order, within a superdomain, there is a correspondence with the columnar VBS and its defects on the square lattice. We noted in the text that in principle J eff for these defects could have the bipartite sign structure that is natural on this embedded square lattice, but that this cannot happen when VBS order is weak. For weak VBS order (and weak quenched disorder) the sign structure of the spin correlations on short scales should respect the symmetries of the triangular lattice. These symmetries are incompatible with the sublattice sign structure associated with the embedded square lattice. If this sign structure is broken on short scales it is unlikely that it will be restored on longer scales.

Vortex nucleation in the weak disorder limit
Here we argue that in the weak disorder limit VBS vortices are necessarily nucleated on the Imry-Ma lengthscale ξ 2d (∆) ∼ exp(J 2 /∆ 2 ), where ∆ is the strength of disorder. Introducing a pair of VBS vortices into a vortex-free state allows a rearrangement of the domain pattern on the lengthscale ξ 2d (∆). If the disorder configuration is favourable this rearrangement can reduce the energy by an amount of order ∆ × ξ 2d (∆). Since ξ 2d (∆) grows exponentially as ∆ → 0, this energy gain is large when ∆ is small, and overwhelms the O(1) core energy cost associated with the vortex core.
The exponential growth of ξ 2d for small ∆ also means that the typical energy scale for the defect-spin couplings J eff is doubly-exponentially weak, since J eff itself is exponentially small in defect separation. For the application to experiments one must consider stronger disorder, and the couplings between unpaired spins will not be exponentially small in any sense.

VBS order as O(n) vison condensate from Z2 spin liquid
In this appendix section we give more details for the discussion of the generalized VBS phases on the triangular lattice through the transition from the proximate Z 2 quantum spin liquid. This is a "starred"-Landau-Ginzburg theory [138], involving condensation of the Z 2 vison. The VBS phases we discuss include not only the columnar-VBS state of Fig. 1 but also the known [76][77][78] plaquette-VBS phases, whose resonances involve various 4-spin, 12-spin, or 16-spin clusters. It is more difficult to visualize a point defect in these complicated phases; nevertheless, we shall now show that these VBS phases all host point topological defects which are Z 2 vortices (i.e. the vortex is its own anti-vortex), and that these Z 2 vortices again necessarily host unpaired spin-1/2s.
The spin liquid variables enable a natural description of the VBS topological defects, in a unified language; the present section will be concerned with various aspects of this description. At the vison condensation transition, the various VBS phases are unified into a continuous manifold of parent VBS states. As described below (see also Refs.76 and 78), for the two simplest vison condensation transitions this manifold is a unit sphere in n-dimensions (n = 4, 6 respectively) but with antipodal points identified, i.e. a headless n-dimensional unit vector. The point defects are then easily understood by analogy to 120-degree magnetic order and nematic orders, whose order parameters are similarly headless: these defects are the Z 2 vortices. [139]. Here a VBS vortex necessarily requires a S = 1/2 spinon to be nucleated in the vortex core; the vortex S = 1/2 modes may be considered as a relic of the fractionalized Z 2 spinon.
Condensing the spinon leads to a magnetically-ordered phase, such as the 120 • ordered phase of the triangular lattice Heisenberg antiferromagnet. The vison does not carry magnetic quantum numbers, so its condensation does not lead to a magnetic (or any time-reversal-broken) phase. Instead, it leads to a lattice-symmetry-broken phase, such as a singlet or plaquette valence bond crystal.
To analyze the manifold of VBS order parameters that results from vison condensation, let us consider the vison condensation in more detail. The vison sees the spin-half spinon as a π flux. Each site of the triangular lattice carries a single spin-1/2 degree of freedom: the number of spinons (in total for both spin flavors) at each lattice site is constrained to be exactly one in the ground state. So the vison experiences a (−1) phase when hopping around any triangular lattice site. The vison hopping can thus be modeled as hopping on a fully frustrated honeycomb lattice, where the honeycomb is the dual lattice formed by triangular plaquettes, and the frustration implies an effective uniform magnetic field experienced by the vison, with flux π per hexagon plaquette. The fully frustrated honeycomb hopping problem is known [76,78] to produce different possible solutions, allowed by the spin liquid PSG, depending on the particular pattern of vison hopping. The simplest vison hopping produces four degenerate minima [real modes from wavevectors ±K/2, ±K /2, where K is a BZ corner], while a more complicated vison hopping (which need not necessarily be more complicated in terms of original spin variables) produces six degenerate minima [real modes from wavevectors ±M 1 /2, ±M 2 /2, ±M 3 /2, where M is the midpoint of a BZ edge]. The family of VBS orders in the latter scenario includes the columnar order [78] with wavevector M . This four-or six-component vector v of the vison Energy cost for fixed and optimized monomers on the square and triangular lattices. Here we compare the energy cost for both fixed-site and partiallyoptimized monomers between the square lattice (square symbols) and triangular lattice (triangle symbols), in a log-linear plot. For fixed defects, the growth of the standard deviation (dashed lines) ensures that optimized defects gain diverging energy, both for bipartite and non-bipartite lattices.
condensate describes the resulting VBS order parameter, with an overall sign redundancy. Though any given VBS order parameter is discrete it is thus natural to describe it as a discrete subset of this continuous manifold, where the manifold is a unit sphere in n dimensions with each pair of opposite points identified into a single point, written as S n−1 /Z 2 ≡ RP n−1 . The resulting manifold RP n−1 has point-like defects, characterized by the first homotopy group π 1 (RP n−1 ) = Z 2 , produced by paths from a point to its diametrically opposite point on the sphere. (Without modding by Z 2 the sphere homotopy is trivial for n > 2.) Thus there is a vortex point defect which is its own antivortex.
We next argue that the vortex core in the VBS orders carries a protected spin-half degree of freedom. To see the spin-half at the vortex core, consider condensing the vison into a VBS configuration with a vortex. The vortex represents a winding of the vison condensation vector v such that v winds from some initial value v 0 at angle θ = 0 near the vortex core, to the opposite value − v 0 at θ = 2π. (Recall that v 0 and − v 0 represent the same VBS order.) The vison field sees the vortex core as a π flux. But there is only one such object in the Z 2 spin liquid: the spinon, which necessarily also carries a spinhalf degree of freedom.

Numerical simulation of random-energy classical dimer model
Here we discuss details of the numerics on the random energy dimer model. We used L×L lattices with periodic boundary conditions in the x direction (e.g. identifying left and right edges of each lattice in Fig. 10), giving a cylinder. Note that this is a planar graph (identical results were found on a torus). We assigned random energies to the edges of the graph, taken from a uniform distribution on the integers {−∆, −∆ + 1, . . . , ∆}, with String of shifted dimers connecting monomers in the random dimer model. The optimal dimer coverings in the random-energy dimer model, computed for the full lattice and for a lattice with a pair of monomer defects, is the same everywhere except along a string connecting the two monomer sites: along the string all dimers are shifted by one. Here we show sample plots of these strings, for fixed defects [L = 80, left] and optimized defects [L = 64, right]. The string fractal dimension is computed to be d f = 1.28 (1) for the triangular lattice, and d f = 1.25(1) for the square lattice. ∆ = 5000. Plots of δE in the main text are given with δE measured in units of ∆. When looking at defects, we removed pairs of sites by removing both a random site from the equator of the lattice (x = x 0 , y = 0) as well as the site x = x 0 +L/2, y = 1 approximately opposite (the shift to y = 1 ensured defects were on opposite sublattices for the square lattice, as required for finding a dimer covering on the remaining graph; the same shift was implemented on the triangular lattice to enable direct comparisons of the results). We computed the complete dimer covering with minimum dimer energy using the famous blossom algorithm of Edmond [140,141] as implemented in the NetworkX python package (implementation by Joris van Rantwijk, "max weight matching" with "maxcardi-nality=True"; weights are interpreted as negative energies). At least 5 × 10 3 independent disorder realizations were computed for each parameter set for the case of optimized defects, with many more (10 5 ) disorder realizations for each parameter set for the case of fixed defects. For the calculation of partially optimized defects, we again restricted to opposite defect pairs on the equator and optimized over x 0 . This reduced the computational effort of the optimization and also reduced finite size effects from open boundary conditions. The change in system energy δE we found for these partially-optimized defects is an upper bound on δE for fully-optimized defects, so it suffices to show that partially optimized defects give a negative δE at large L.
Let us briefly give further technical details for the histograms shown in Fig. 2 of the main text. These are histograms of the numerically computed distribution of system energy change δE under introduction of a pair of monomers spaced L/2 apart, that are either fixed [top] or allowed to optimize their position along the equator of the cylindrical system [bottom]. Colors shown correspond to the various system sizes L = 8, ..., 64 with rainbow colors from blue to red. Histogram bin size is 0.05∆ and 0.15∆ for the fixed and optimized monomers respectively, where 2∆ is the width of the uniform distribution of graph edge weights (−∆, ∆); it is implied that δE is in units of ∆. Observe that for fixed monomers, the histograms are all numerically close to Gaussians, with mean and variance that converge quickly as L → ∞, and with a small tail for negative δE. For optimized monomers, the histograms (no longer Gaussian-like) narrow and shift towards negative δE as L increases towards the thermodynamic limit. The behavior of the peak with increasing L is consistent with the scaling ∼ −(log L) 1/2 expected for minima obtained by sampling a polynomial-in-L number of times from a fixed-defect distribution whose tail is ∼ exp −const. δE 2 .
For the sake of comparison with known results on the bipartite case, Fig. 9 shows the energy cost for fixed and optimized defects on the square and triangular lattices. On bipartite lattices, a fixed monomer costs energy log L since it is a long-ranged vortex, with log L interactions in the elastic medium which maps to this problem. This elastic-vorticity effect indeed does not occur for the nonbipartite case.
It is also interesting to study more subtle properties of the system with monomers. Introducing defect monomers modifies the optimal dimer configuration along a string which connects the two monomers; as shown in Fig. 10, this string is a fractal object. The energy gain due to this domain wall string excitation enables the monomers, on either bipartite or nonbipartite lattices, to be pulled down from arbitrarily high energies and nucleate on the lattice.

Defects in the classical dimer model on a non-bipartite lattice
Here we discuss monomer defects in the classical dimer model, focusing on the case of non-bipartite lattices.
A defect can be regarded as one endpoint of an excitation which is an open line instead of a closed loop. The typical energy cost or gain due to disorder can be parsed into contributions from sections of this line in a series of concentric annuli of radii k ∼ 2 k with k = 1, 2, . . .. The contribution from the kth lengthscale will be of random sign and of order θ = 2 −k|θ| . The sum of these contributions is convergent, giving a random variable with mean and variance of order one.
We can also construct a rare region which shows that the density for defects in the classical triangular lattice dimer model, at zero temperature, is at least exp[−c(K/∆) 2 ] for large K, where K is the core energy cost of a monomer defect. Consider a disc-like region D of size , where we have fixed the signs of all the bond energies to favour a particular dimer configuration in D with a defect at the center. The probability of such a region is exponentially small in the area, e −const. 2 . To remove the defect at the center, we must disturb the favoured configuration inside D along a line of linear extent ∼ connecting the center of D to its exterior.
This costs an energy of order × ∆ with ∆ the disorder strength. Therefore, for large , such regions correspond to places where the energy 'cost' for introducing the defect, relative to the ground state of the defect-free model, is δE ∼ K − const. ∆. By taking ∼ K/∆ we obtain a region where it is favourable to introduce a defect. This shows such regions have a density ρ at least as large as exp[−c(K/∆) 2 ] with c a numerical constant. Note that the rare regions necessary for this argument do not require the disorder distribution on an individual link to be unbounded.

Defining singlets with spin-orbit coupling
Here we consider the effects of spin-orbit coupling on microscopic definitions of two-spin singlet states; recall that VBS phases are generally defined as lattice-broken quantum paramagnets, which preserve time reversal symmetry but transform non-trivially under some lattice translations and rotation symmetries.
The strong spin-orbit coupling of Yb 3+ completely breaks the continuous SO(3) spin rotation symmetry, which naively suggests that the spherically-symmetric two-spin singlet state is no longer a well-defined eigenstate. However, as long as inversion symmetry is present, this is not the case. The singlet, which is odd under inversion, does not mix with the inversion-even triplet manifold. (Extensions to the case of spin orbit coupling without inversion centers will be discussed elsewhere [111].) To see this explicitly, consider a given bond (i, j) with its 3 × 3 matrix of spin interaction coefficients J µν ij . The matrix J µν is real and, if it is symmetric (e.g. if the bond has an inversion center), it can be diagonalized to yield realvalued eigenvectors, known as its principal axes. (The symmetry condition on J µν is equivalent to ignoring any Dzyaloshinskii-Moriya interactions, which may arise in particular disorder configurations but are forbidden if inversion symmetry is preserved on average.) The spin interaction is then written as a sum of exchanges through the principal axes x, y, z, i.e. J x S x S x +J y S y S y +J z S z S z with different coefficients J x , J y , J z . The singlet state (|↑↓ − |↓↑ ) is an eigenstate of each of the termsthis is clear for S z S z , and easy to see from the singlet's spherical symmetry also for S x S x and S y S y . Thus if the J x , J y , J z coeffients are on average sufficiently antiferromagnetic then the ground state of the two spins (i, j) will be the spin singlet state, despite the strong SOC.

Spin-pair approximation
Let n(E) be the number of singlet pairs with energy less than or equal to E, i.e. the function plotted in Fig. 6. The spin-pair approximation consists of treating each singlet pair as completely decoupled from its environment. The resulting formulas for susceptibility and specific heat are as follows, When the coupling between two spins has E < 0, corresponding to a ferromagnetic sign coupling and a triplet ground state, the integrands are given with energy −E.

Microscopic considerations for YbMgGaO4
As mentioned above, the experimentally observed neutron scattering peaks along the BZ boundary (also seen in the corresponding theoretical computation shown in Fig.5) require some small but nonzero amount of secondneighbor correlations. In particular, it is easy to see that pure first-neighbor correlations do not reproduce the correct peak locations even qualitatively (see Fig.11). Second neighbor correlations are weak but present.
They could be supported by a possible second neighbor exchange, which though would be unusual for f moments, could potentially arise via Yb-O-O-Yb exchange pathways via direct overlap of in-plane-oriented p-orbitals on adjacent oxygen ions. Even at J 2 = 0 strong second neighbor correlations can arise. Let us consider the nearest neighbor interactions in detail.
The spin Hamiltonian for YbMgGaO 4 in the clean limit involves multiple exchange terms whose form is determined by the crystal symmetries. The question of the exact form of the Hamiltonian has not been fully settled. [5,8,11,13,142] We note (see also Ref.32) that the so-called J ±± and J z± terms may be re-written through a pseudo-dipole exchange I and a triangular lattice Kitaev term K, H SOC ij = I( S i · r ij )( S j · r ij ) + K( S i · γ ij )( S j · γ ij ) (45) where r ij is the unit vector connecting sites i and j, and the Kitaev label γ ij ∝ r Yb-O1 × r Yb-O2 is defined as the unit vector perpendicular to both Yb i -O 1,2 -Yb j exchange paths, i.e. perpendicular to the plane spanned by Yb sites i, j and the two intervening oxygen ions. (On the triangular lattice, γ ij is fully determined by the bond orientation r ij , and { γ ij } form an orthonormal coordinate system whose (1, 1, 1) axis is normal to the lattice plane.) The terms can be related to the original formulation[2] as follows: Notation : J z , J ± , J ±± , J z± (46) Heisenberg J = 1, 1 2 , 0, 0 Pseudodipole I = 0, 1 4 , Kitaev K = 1 3 , For example, on a bond for which the pseudodipole exchange is alongx, the Kitaev exchange is then along 2/3ŷ + 1/3ẑ.
These interactions are together strongly frustrating, and it is quite conceivable that once some preferred nearest neighbor singlet states are formed, other nearby sites prefer to order into configurations that involve second neighbor singlet states; or, that resonances among multiple sites are formed, which entail both first and second neighbor correlations in some ratio. Regardless of the microscopic origin of the observed weak second-neighbor correlations, the point in the main text regarding interpreting the experimental neutron scattering is that it can be captured by a simple model of short ranged randomlyoriented singlets with a particular small ratio of second to first neighbor correlations.
Regarding promising regions in parameter space for stabilizing non-magnetic phases, we note the following fact about the Hamiltonian above: the parameter point K = −2, J z = J xy = +1 is known (via a transformation into an exactly solvable ferromagnet) to produce a stripy antiferromagnetic order which, at this parameter point, has zero quantum fluctuations [143,144]. Though the stripy phase should therefore extend over a large region of parameter space [32], tuning away from the fluctuationsfree point may be a useful route for destabilizing magnetic order. In this context it is also useful to observe that when upon adding Kitaev interactions to the Heisenberg antiferromagnet 120-degree order, the resulting classical ground state [145] involves a finite but low density of nucleated Z 2 vortices (the topological defects of the 120 degree order), with complicated long-wavelength spin configurations reminiscent of skyrmion crystals. The true quantum fate of such configurations upon further perturbations, though they may be difficult to study using unbiased quantum algorithms, may result in a quantum paramagnet phase.