Random-Singlet Phase in Disordered Two-Dimensional Quantum Magnets

We study effects of disorder (randomness) in a 2D square-lattice $S=1/2$ quantum spin system, the $J$-$Q$ model with a 6-spin interaction $Q$ supplementing the Heisenberg exchange $J$. In the absence of disorder the system hosts antiferromagnetic (AFM) and columnar valence-bond-solid (VBS) ground states. The VBS breaks $Z_4$ symmetry, and in the presence of arbitrarily weak disorder it forms domains. Using QMC simulations, we demonstrate two kinds of such disordered VBS states. Upon dilution, a removed site leaves a localized spin in the opposite sublattice. These spins form AFM order. For random interactions, we find a different state, with no order but algebraically decaying mean correlations. We identify localized spinons at the nexus of domain walls between different VBS patterns. These spinons form correlated groups with the same number of spinons and antispinons. Within such a group, there is a strong tendency to singlet formation, because of spinon-spinon interactions mediated by the domain walls. Thus, no long-range AFM order forms. We propose that this state is a 2D analog of the well-known 1D random singlet (RS) state, though the dynamic exponent $z$ in 2D is finite. By studying the T-dependent magnetic susceptibility, we find that $z$ varies, from $z=2$ at the AFM--RS phase boundary and larger in the RS phase The RS state discovered here in a system without geometric frustration should correspond to the same fixed point as the RS state recently proposed for frustrated systems, and the ability to study it without Monte Carlo sign problems opens up opportunities for further detailed characterization of its static and dynamic properties. We also discuss experimental evidence of the RS phase in the quasi-two-dimensional square-lattice random-exchange quantum magnets Sr$_2$CuTe$_{1-x}$W$_x$O$_6$.


I. INTRODUCTION
In the quest to classify and characterize ground states and excitations of quantum many-body systems, disorder (quenched randomness) plays a central role. Beyond the fundamental scientific interest in understanding the interplay between quantum fluctuations and intrinsic randomness, there are also potential practical implications: In the same way as pure crystalline states of matter are often not optimal for achieving desired properties of materials, e.g., in the case of metals hardened by limiting the size of crystal grains, it is likely that quantum technologies will emerge that exploit disorder effects. For example, random spin chains have been proposed as key elements for memories [1,2] and state transfer channels [3] in quantum computing. Two-dimensional (2D) quantum spin systems, which we consider here, is another natural setting for exploring novel disorder-induced states. Recent experimental efforts have been devoted * yc.lin@nccu.edu.tw † waguo@bnu.edu.cn ‡ sandvik@bu.edu to searches for quantum spin liquids in quasi-2D insulators. Several candidate systems showing the qualitative signatures of spin liquids have been identified, e.g., in a series of organic salts where the spins reside on triangular lattices [4][5][6][7] and in the kagome-lattice herbertsmithite [8][9][10][11]. It has so far not been possible to unambiguously match the properties of these systems to theoretically proposed spin liquids, however, and it has been proposed that disorder effects are crucial for understanding the observed behaviors [12]. In a more extreme interpretation put forward recently [13][14][15][16][17][18][19], disorder is even responsible for realizing a certain spin liquid, the random singlet (RS) state, in some triangular, kagome, and frustrated honeycomb lattice systems, e.g., the triangular YbMgGaO 4 [20,21] where disorder is present in the form of random occupation of Mg and Ga ions on equivalent lattice sites between the magnetic layers. While such a state has not yet been observed in systems without geometric frustration, there is recent experimental evidence for an RS state in a square lattice system; the double perovskite Sr 2 Cu(Te 0.5 W 0.5 )O 6 . Here the disorder is in the form of random Te↔W substitutions relative to the isostructural compounds Sr 2 CuTeO 6 and Sr 2 CuWO 6 , which have dominant nearest-and next-nearest-neighbor spin interactions, respectively [22]. We will here show that disorder can induce a spinliquid-like state-a gapless state with algebraic correlation functions-in a 2D quantum spin system on the square lattice even without geometric frustration. To this end, we carry our large-scale quantum Monte Carlo (QMC) simulation studies of an S = 1/2 quantum spin model, the J-Q model, which in the absence of disorder hosts both a Néel antiferromagnetic (AFM) and a spontaneously dimerized valence-bond solid (VBS) ground state. The transition is driven by enhancing the formation of correlated singlets by increasing the multi-spin (here six-spin) interaction Q, which competes with the Heisenberg exchange J. We show that randomness in the coupling constants leads to the formation of domains in the four-fold degenerate VBS state, with different realizations of the bond order and with domain walls of the type expected [23] to lead to localized spinons at each nexus of four domain walls. These spinons form in correlated groups of even numbers, as a consequence of the domain-wall topology. We will show evidence for domain-wall mediated enhanced spinon-spinon interactions, which leads to singlet formation within the groups and no residual AFM ordering of the spinons. Due to rare events, in the form of singlet formation over large distances, the spin and bond correlations decay as power laws. We present detailed characteristics of this RS spin liquid, in the ground state as well as at elevated temperatures. As a contrast, we also consider a site-diluted system, in which the remnant local moments associated with removed sites are not spatially strongly correlated; thus residual AFM order forms and there is no RS phase.
We will argue that the RS state we identify here is the same one, in the sense of renormalization group (RG) fixed points, as the one proposed recently to arise out of a VBS on the triangular lattice in the presence of random couplings [17]. It may then also be a realization of the unusual magnetic states observed in YbMgGaO 4 and Sr 2 Cu(Te 0.5 W 0.5 )O 6 , and possibly in many other disordered spin liquid candidates as well. The possibility of creating this state with a "designer Hamiltonian" within the J-Q family of models is very significant, as this unfrustrated (in the geometric sense) system is amenable to large-scale QMC studies (without the "sign problems" plaguing simulations of models with frustration), thus allowing the state to be characterized essentially completely-far beyond the analytical calculations in Ref. 17 and the exact diagonalization (ED) numerics on small frustrated Heisenberg lattices in Refs. 13-16 and on slightly larger triangular lattices by density-matrix renormalization group (DMRG) calculations in Ref. 19. In particular, we are able to reliably study the AFM-RS quantum phase transition. This paper is devoted to a broad survey of the phase diagrams, quantum phase transitions, and basic ground state and temperature T > 0 properties of the 2D RS phase in different versions of the random J-Q model.
The paper is organized as follows: In Sec. II we discuss the broader context of our work and provide specifics of the models considered. In addition to the main focus on different kinds of disorder in the J-Q model, we will also discuss a simpler case as a point of reference: the statically columnar-dimerized Heisenberg model in which localized moments different from the VBS spinons form in the neighborhood of removed sites. In order to aid in the presentation and interpretation of the extensive QMC results of the later sections, in Sec. III we first discuss qualitatively the phenomena and mechanisms we have identified and quantify later. In Sec. IV we present results of ground-state projector QMC calculations of static properties of all the models considered, with the main focus on the order parameters and correlation functions in the RS phase in the cases where this state is attained. We demonstrate the existence of a universal continuous AFM-RS quantum phase transition. In Sec. V we discuss results at T > 0 which allow us to extract the dynamic exponent z of the RS phase. In Sec. VI we provide evidence for the mechanism underlying the formation of the RS state; spinon interactions mediated by VBS domain walls. We conclude in Sec. VII with a brief summary and further discussion of our results and their significance in the context of both theory and experiments.

II. BACKGROUND AND MODELS
A. Infinite-randomness fixed points and the random singlet phase Theoretically, when randomness is a relevant perturbation under RG transformations, fixed points corresponding to ground state phases and critical points appear beyond those realized in pure, translationally-invariant systems [24,25]. In some cases the RG flow converges to non-zero but finite disorder, e.g., at critical points in many quantum spin glasses [26][27][28][29], boson systems with random potentials [30] or random hopping [31], and Heisenberg antiferromagnets [32][33][34]. However, the randomness can also increase without bounds in the RG flow, leading to an infinite-randomness fixed point (IRFP). This broad class of fixed points has been extensively studied using strong-disorder RG (SDRG) methods in quantum systems in one [35][36][37][38][39][40][41][42] and higher dimensions [43][44][45][46][47] (in addition to many applications in classical statistical physics [48]). The most striking general property of the IRFPs is an infinite dynamic exponent z, i.e., the scaling relationship between energy (ǫ) and length (l) scales is exponential instead of the conventional power-law relation ǫ ∼ l −z . Moreover, rare instances of long-distance entangled spins (or particles) lead to different behaviors of the mean and typical correlation functions versus distance r, decaying, respectively, as a power law and exponentially.
An important example of an IRFP is the 1D RS phase, realized in the antiferromagnetic Heisenberg chain with random couplings [35,36,40,42], where long-distance entangled spins lead to the mean spin correlations decay-ing as r −2 [36,42] (while the typical correlations decay exponentially) and the entanglement entropy diverging logarithmically with the system size [39].
IRFPs have been identified also in 2D systems, primarily in transverse-field Ising models [45][46][47] but also in experiments on the superconductor-metal transition in Ga thin films [49]. However, no convincing case of such a phase or critical point has been reported in 2D quantum magnets with spin-isotropic interactions, such as the standard Heisenberg exchange, as far as we are aware. If an RS state exists in such systems, one would expect it to have algebraically decaying mean correlation functions, as in the 1D case, and if the state also corresponds to an IRFP the dynamic exponent would presumably be infinite as well. However, an RS state can also in principle exist which has finite z, although such a state corresponding to an RG fixed point at finite disorder strength does not exist in random Heisenberg chains. Finite-disorder fixed points have been obtained in SDRG calculations on the 2D Heisenberg model with various types of disorder [32,33], but it is not clear whether the SDRG method, by its construction and underlying assumptions, produces the correct fixed point when it does not flow to infinite disorder.
As mentioned in the Introduction, Sec. I, there are some experimental indications of 2D disorder-induced spin liquids with finite z in frustrated quantum magnets, according to interpretations supported by numerical studies of the S = 1/2 Heisenberg antiferromagnet with random couplings on the triangular and kagome lattices [13][14][15]19], and also on the honeycomb lattice with frustrated (J 1 -J 2 ) interactions [16]. These may very well be realizations of an RS state, as proposed. However, a full characterization of the putative RS ground state and its low-temperature thermodynamic properties (i.e., the form of the asymptotic long-distance correlations and the value of the dynamic exponent) was not possible, because of the limited lattice sizes accessible to ED [13][14][15][16] and DMRG [19] calculations. The recently developed theory of the RS state arising out of a VBS on the triangular lattice [17] contains ingredients not discussed in the context of the numerical works, in particular the role of VBS domains and localized spinons. It is still likely that the states discussed within these two approaches are actually the same finite-z RS state.
Here we consider a class of S = 1/2 quantum spin models on the 2D square lattice, with no geometric frustration but with interactions leading to weakened AFM order or nonmagnetic VBS states on uniform lattices. In systems with random couplings, the dynamic exponent is finite and varying throughout the RS phase, which is a clear indication of a class of finite disorder RG fixed points. Our results suggest a mechanism of pairing of localized spinons, which leads to the RS state instead of a weakly ordered AFM state (which has been regarded as the most likely state forming in the random VBS in the absence of frustrated interactions [17]). Importantly, this RS state in an unfrustrated, bipartite system can be induced also in cases where the pure host system is not yet in the VBS state (though not in the standard Heisenberg model with random nearest-neighbor couplings [50]), because local VBS domains are still created in response to the disorder. This observation, along with other considerations, also lends support to a universality scenario that connects our square-lattice RS state directly to the above mentioned states studied with ED and DMRG in frustrated lattices with various host states [13][14][15][16]19] and to that arising out of the triangular-lattice VBS state [17].

B. Random singlet state in the 2D J-Q model
We study a square-lattice Heisenberg antiferromagnet with nearest-neighbor exchange J augmented with certain multi-spin interactions of strength Q (the J-Q model). The unadulterated translationally invariant model is defined by the Hamiltonian [51,52] where P ij is the singlet projector for two S = 1/2 spins, In the sums in Eq. (1), ij indicates nearest-neighbor sites, and the index pairs ij, kl, and mn in ijklmn are neighbors forming a horizontal or vertical column, as illustrated in Fig. 1. The summations are over all pairs and columns, so that the Hamiltonian respects all the symmetries of the square lattice, including the 90 • rotation symmetry when J x = J y = J and Q x = Q y = Q as we have assumed in Eq. (1). We will introduce various forms of disorder in the model, including site dilution and random J and Q couplings drawn from suitable distributions; detailed definitions of the different cases are presented in Sec. IV. In the uniform system the Q interactions compete against the exchange terms J, disfavoring the strong AFM order present for Q = 0 (the standard 2D Heisenberg model [53]) by producing correlated local singlets. The interactions are not frustrated in the standard (geometric) sense, however, and the model is amenable to large-scale QMC simulations for all positive values of the ratio g = Q/J (with J ≥ 0, Q ≥ 0 being of primary interest) [54]. The ground state has AFM order for g < g c , with g c ≈ 1.50 [51], and is a spontaneously dimerized valence-bond solid (VBS) for g > g c . In the VBS phase the Z 4 symmetry of four degenerate columnar dimer patterns is broken when L → ∞.
A columnar VBS state and an AFM-VBS transition is also realized if the Q-interaction in Eq. (1) is replaced by a simpler one with only two singlet projectors [55]. The critical coupling ratio g c is then much larger, g ≈ 22, and the VBS order is much weaker throughout the phase. A large number of studies have been devoted to the issue of deconfined quantum criticality within this model [55][56][57][58][59][60][61][62]. Disorder effects on the VBS state are easier to study with the more extended Q term in Eq. (1), and we will here demonstrate RS behavior for a significant range of coupling ratios when either the J or the Q interactions are random. We expect these disorder effects to be generic for VBS phases on bipartite lattices, and, we will also argue that the RS phase we identify is even more generic, most likely being realized also in geometrically frustrated quantum spin systems.
Before the advent of the J-Q model, VBS physics was normally associated with geometric frustration, in models such as the J-J ′ Heisenberg model with nearest-(J) and next-nearest-neighbor (J ′ ) couplings. These systems are not amenable to large-scale QMC studies because of mixed-sign sampling weights (the sign problem), except at the variational level in sampling and optimizing wave functions [63,64]. While great progress has been made in the last several years on DMRG and methods based on tensor product state (TNS) for studying frustrated models (see e.g., the recent papers [65][66][67] for applications to the J-J ′ Heisenberg model), various convergence issues or limited system sizes still make it impossible to carry out calculations as reliable as QMC simulations of sign-problem free models.
The J-Q models exhibit many of the phenomena of long-standing interest in the context of frustrated quantum magnetism, in particular the AFM-VBS transition [68], which appears to realize the exotic deconfined quantum-critical point (DQC) scenario [69,70]. While it is presently not clear whether exactly this transition is also realized in the J-J ′ Heisenberg model [65][66][67], the phenomenon has attracted a great deal of interest as it is a prominent example of a quantum phase transition beyond the standard Landau-Ginzburg-Wilson framework. The J-Q models offer unique opportunities to study the emergent degrees of freedom-spinons and gauge fieldsthat are the ingredients of the field-theory description of the DQC point. A very interesting question is how these degrees of freedom respond to quenched disorder, and this is the topic of the present paper.
By the Imry-Ma argument [71], in the presence of even an infinitesimal degree of randomness in the local interactions, the VBS can no longer exist as a long-range ordered state, due to different columnar dimerization patterns being energetically favored in different parts of the lattice.
Thus, the uniform VBS breaks up into finite domains of different VBS patterns. An extreme case (in the sense of very small VBS domains) of such a disordered dimer state has been dubbed the valence-bond glass (VBG) [72]. It essentially consists of a random arrangement of short valence bonds and it has been discussed in the experimental context of the kagome-lattice material herbertsmithite [8,9], and also in 3D frustrated spin systems [73,74]. The kagome spin S = 1/2 lattice of the herbertsmithite is to some degree diluted with non-magnetic impurities, and these also liberate spinons from the singlet ground state [12]. It was argued that these spinons interact and form a gapless critical RS state. In this case the spinons can be regarded as a byproduct of the dilution, and in the original picture of the VBG without dilution [72] there were no such spinons.
In analogy with one dimensional spin chains with VBS ground states [42,75], and considering the nature of the elementary domain walls in 2D VBS states [23], one should expect a VBS broken up into domains to also have localized spinons at the nexus of domain walls. Therefore, interesting magnetic properties due to local moments can arise even without the explicit introduction of moments by dilution. Indeed, it was very recently argued [17] that a spin-liquid-like RS state arises in this way on the triangular lattice when the starting pristine system is a VBS. The RS state there is formed as a direct consequence of the randomly interacting localized spinons at the nexus of domain walls, while spinons do not appear in the scenario discussed in the context of the ED [13][14][15][16] and DMRG [19] studies. Localized spinons may still give rise to the physical properties observed in the numerical calculations, although spinons and VBS domains were not studied explicitly (which would also not be easy with the very small lattices considered). On the square lattice with bipartite interactions, this kind of state has not been previously expected, however, and it was argued that the most likely scenario for systems like the random J-Q model is that the liberated spinons form a subsystem with AFM order, instead of a fully disordered RS state [17]. This is at odds with our results to be presented here.
An example, illustrated in Fig. 2, of a well understood system in which residual AFM forms among impurity spins is the diluted columnar dimerized Heisenberg model, which we will later use as a bench-mark case for our numerical analysis techniques. For sufficiently large ratio j 2 = J 2 /J 1 of the intra-to inter-dimer couplings, in the quantum paramagnetic phase, the removed sites leave behind 'dangling' spins at the sites near the broken dimers, and these form a subsystem with AFM order due to effective bipartite interactions mediated by the inert spin-gapped dimer host [76]. Thus, the quantum phase transition out of the AFM ground state, at j 2 ≈ 1.91 in the intact system [54,77,78], is destroyed and replaced by a cross-over from strong to weak AFM order. In a VBG on the square lattice, one might imagine that the disorder induced spinons should be subject to a simi-  [54,77,78]), the ground state is approximately a product of singlets on the strong bonds, and when diluted the 'dangling spins' remaining at the 'broken dimer' adjacent to each removed spin constitutes a localized magnetic S = 1/2 moment. lar ordering mechanism [17]. However, our results and arguments suggest that the correlated nature of spinonantispinon pairs (and larger complexes of even numbers of spinons) in the randomized VBS was not taken fully into account previously. In particular, we argue that a key missing ingredient in the analysis of bipartite systems Kimchi et al. [17] is that the VBS domain walls act as channels of enhanced spinon-spinon interactions within the groups of even numbers of spinons, thus leading to stronger than expected tendency to local singlet formation and no residual AFM ordering.
Though it is not immediately clear whether the RS phase that we identify and characterize here corresponds to the same fixed point as the state identified on the triangular lattice by Kimchi et al. [17], this certainly is a strong possibility based on symmetry considerations. Moreover, similar to our results presented here, the random state on the triangular lattice does not have infinite dynamic exponent, but exhibits power-law correlations in both space and time. We further show that the RS state can also form in some cases even though the bipartite host system is not yet VBS ordered, still in the AFM state, as long as there are sufficient interactions (here Q terms) favoring the formation of some local VBS domains. This role could also be played by standard frustrated interactions, and it therefore appears most likely that the RS state in the disordered J-Q model actually is the same as those states discussed previously in the context of a variety of frustrated host systems, including the ED studies [13][14][15][16] and DMRG calculations [19]. In these numerical works, the physical picture presented for the nature of the RS state was different, however, with an emphasis in Refs. 13-16 put on the singlet pairs (Anderson localization of singlets) [15] and no reference to the localized spinons and VBS domains. These are actually the objects that form the key ingredients in the theory In Secs. IV and V we will present QMC results for the Hamiltonian Eq. (1) with random J and random Q couplings, as well as for a site site diluted system with no randomness in the remaining J and Q interactions. For reference we also present results for the diluted J 1 -J 2 Heisenberg model. To characterize the ground states of these systems in an unbiased way, we use a ground-state projector QMC method formulated in the valence-bond basis [79], and to obtain properties at temperature T > 0 we use the stochastic series expansion (SSE) method [80]. To make the results sections more accessible and concise, in Sec. III we first outline the physical scenario that arises out of the many different calculations reported in the subsequent sections.

III. DOMAIN WALLS AND SPINONS IN THE DISORDERED VALENCE-BOND SOLID
On the 2D square lattice and with the bipartite nature of a model such as the J-Q model, the main question regarding the disordered VBS state is whether the spinons localizing at each nexus of four domain walls [23] will form long-range AFM order or some other collective state with only short-range or algebraic spin-spin correlations. As already discussed in Sec. II B, one might suspect [17] that AFM order should exist for all values of g = J/Q, in analogy with the fate of the quantum paramagnet and Néel-paramagnetic quantum phase transition in Heisenberg models with static dimerization when spins are randomly diluted (Fig. 2). This picture neglects important spatial correlations among the localized spinons, however, as well as the nature of the VBS domain walls that connect the spinons.
To understand these spinon correlations, consider first an individual, localized spinon. As illustrated in Fig. 3, the four lattice bonds pointing out from the site of an unpaired spin (the core of the spinon) correspond to the four different VBS patterns. While a different pattern can also in principle form, with the bonds rotated by 90 • relative to those in the figure [23], our simulations of the J-Q model consistently show the 'star' pattern at the spinon (but this local arrangement should not change the properties of the domain walls discussed in Ref. [23]). The four bonds and the corresponding extended VBS domains can be associated with angles φ as indicated. Note that the energetically favored domain walls correspond to a π/2 phase twist [23], while walls with π phase change are unstable and break up into two π/2 walls (as shown explicitly in Ref. [81]). This is the origin of the proper classification of the symmetry of the VBS as a Z 4 , or 'clock' symmetry (as opposed to the full S n permutation symmetry if all domain walls were equivalent) [23,69]. Within a domain wall, the angle φ (properly defined by coarse graining) changes continuously.
Note that a spinon can be associated with either sublattice A or B, and the way in which the angle φ changes, increasing or decreasing, when going around the spinon in a given direction depends on the sublattice. Thus, we can also refer to the two cases as a spinon and an antispinon (but for convenience we will often just use the term spinon for both). This classification remains valid also in the presence of longer valence bonds, as long as only bonds connecting the two sublattices are allowed. This is exactly the case with bipartite interactions, where bonds connecting sites on the same sublattice are always eliminated when a state written in the valence-bond basis is time-evolved.
As also stressed in Ref. [17], when starting from the clean VBS, spinons always have to be introduced in pairs of spinons and antispinons. When separating the two members of a pair, domains form such that each spinon is connected to all four types of domains as in Fig. 3. As shown in Fig. 4, this leads to a four-stranded confining string, akin to the (more complicated) quark-confining strings in quantum chromodynamics [82]. Here we have not shown the details of the bonds within the domains, only the colors corresponding to the coding in Fig. 3. As already mentioned, in principle there will also be valence bonds of length greater than one lattice spacing, but the pictures remain valid as long as the probability of longer bonds decays sufficiently rapidly with the bond length. If we consider the total-spin singlet ground state, there will also be a bond connecting the spinon and the antispinon sites. Such a long bond corresponds to a small gap to the triplet; vanishing in the limit of large separation. In the non-random VBS, the spinons can not actually be far separated in this way, because other spinons can be excited from the vacuum (the VBS ground state) as the string energy becomes sufficiently high; thus the confining strings will break, which again is analogous to the case of quark confinement.
In a system with random couplings, different VBS angles φ ∈ {0, π/2, π, 3π/2} will be favored in different parts of the system and the domain size will be governed by the competition of the energy cost of the domain walls and the energy gains due to the disorder. In classical systems, according to the Imry-Ma argument [71], this always leads to domain formation at T = 0 in dimensionality D < 2, while for D > 2 the uniform state is stable in the presence of weak disorder. Considering entropy effects, the uniform state is also unstable at T > 0 in D = 2. Similarly one can expect quantum fluctuations to also always lead to domain formation in systems with two spatial dimensions at T = 0 [17]. At least for weak disorder, the domain walls should still be of the π/2twist type. In addition to single domains forming with this phase difference with respect to their surroundings over the whole length of the boundary, the domain-wall topology also allows for a different situation if localized spinons are allowed to form. As in the uniform VBS state discussed above, spinons forming in a VBS broken up into domains must also always appear in groups of an even number-half of the spinons and half of them antispinons. In Fig. 4, a quadruplet is shown along with the spinon pair already discussed. It is this inherent correlation among spinons and, importantly, the tendency to singlet formation within the groups, that we believe prohibit the formation of AFM order in the random VBS (which we will show is actually the RS phase) arising out of the VBS in the J-Q model. The effective interactions between the spinons should be mediated through the domain walls (and we will show explicit evidence for this), because they have much smaller local mass gaps than the bulk of the VBS domains (through which interactions between different spinon groups have to be mediated). We will also later comment on this picture in the context of SDRG theory.
According to our findings in Sec. IV, the above described disordered VBS is an RS state with mean spin correlations decaying with distance as r −2 . It arises out of the VBS state in the J-Q model with random cou- In the pure models, Λ = 0, there is a DQC point (red circles) separating the AFM and VBS phases. The VBS is destroyed, breaking up into domains, for any Λ > 0. In (a), which applies to the model with site dilution, there is no phase transition vs the coupling ratio Q/J when Λ > 0, only a cross-over (indicated by the wedge) between the standard AFM state and a VBG state with finite VBS domains in which weak AFM order forms among localized effective moments. In (b), which applies to the case of random coupling constants, there is a true continuous quantum phase transition between the AFM and RS phases for at least some range of Λ > 0.
plings (either random J or random Q, both of which we will study, or all random, which we have not considered). The form of the spin correlation function is, thus, the same as in the 1D RS phase, and the dimer (bond singlet) correlations decay with a higher power, likely r −4 , which again would be the same as in 1D [42]. Unlike the 1D RS state, but in agreement with the VBG state proposed on the diluted kagome lattice [12] and in the RS states proposed more recently [13][14][15][16][17], we do not find a divergent dynamic exponent. By investigating the temperature dependence of the uniform magnetic susceptibility we find z = 2 (T independent susceptibility) at the AFM-RS phase boundary and z > 2 (power-law divergent susceptibility) inside the RS phase.
In further support of a disordered VBS state with no AFM order, we also compare the model with random couplings with a site-diluted J-Q model. Here, like in the diluted J 1 -J 2 model in Fig. 2, there will be effective moments associated with the removed sites. Thus, while there may also be localized correlated spinons associated with the meeting points of four domain walls, now there are also moments at random locations and those are not subject to the strong singlet formation within groups of spinons. Indeed, in this case we find a VBS broken up into domains and weak AFM order, and no RS state exists in the phase diagram.
In Fig. 5 we summarize the two kinds of phase diagrams that we find for the J-Q model in the presence of the different types of disorder discussed in this paper. We expect these phase diagrams to be generic for disordered 2D quantum magnets that host AFM-VBS quantum phase transitions in the absence of disorder. Note the way the AFM-RS phase boundary has been drawn in Fig. 5 as tilted into the AFM phase, i.e., one can reach the RS state not only from the VBS phase of the pure system but also (for some types of disorder) from the AFM state even when it is quite far from the AFM-VBS transition. This is interpreted as the tendency to local VBS domain formation being favored by the disorder. On the square lattice the Heisenberg model with only nearest-neighbor couplings J, disorder in the form of random unfrustrated Js does not induce an RS phase [50], and a critical strength of frustrated interactions is presumably required, like in the other frustrated systems, to induce it [13][14][15][16][17][18][19]. The Q interactions of the J-Q model explicitly favor local correlated singlets and apparently mimic the effects of geometrically frustrated interactions in their ability to generate the RS state.

IV. GROUND STATE PROPERTIES
We here present QMC results for the J-Q model defined in Eq. (1) in the presence of disorder in the form of random J or random Q. In some cases we use a bimodal distribution of couplings J ij ∈ {J(1 − ∆), J(1 + ∆)} or Q ijklmn ∈ {Q(1 − ∆), Q(1 + ∆)}, with equal probability for the two values, and in other cases we consider uniform distributions with the couplings bounded by the above values. To contrast random couplings and site dilution, we also consider the J-Q model where a given fraction of the sites, randomly selected, are missing. All operators in Eq. (1) touching one or several missing sites are removed from the Hamiltonian. To bench-mark our calculations for the J-Q model against a case where it is known that site dilution induces AFM order in a quantum paramagnetic host, we also consider the diluted statically dimerized Heisenberg model illustrated in Fig. 2. In all cases, we average QMC results over a large number of independent realization of the disorder (hundreds to thousands) on square lattices with N = L × L sites and periodic boundary conditions. Below, in Sec. IV A we will first briefly describe the QMC algorithm used in the ground state calculations and also introduce the main observable we use to characterize the systems. In the following subsections, we present results for all the models; the diluted J 1 -J 2 model in Sec. IV B, the diluted J-Q model in Sec. IV C, and the random J and random Q systems in Sec. IV D and Sec. IV E, respectively.

A. Ground state projector method
The QMC method we use here projects out the ground state from a trial wave function |Ψ(0) written in the valence bond basis consisting of all possible tilings of the square lattice into bipartite singlet bonds. Acting with (−H) m on this state, we obtain an un-normalized state |Ψ(m) ; thus expectation values of operators A are evaluated in the form for sufficiently large m. The different propagation paths contributing to |Ψ(m) are sampled by expressing H as a sum over the J and Q terms in Eq. (1) and carrying out Monte Carlo updates on the corresponding strings (products) of m such operators acting on |Ψ(0) . In this process, the spin degrees of freedom are put back in by also sampling the ↑↓ and ↓↑ contributions to each valence bond (where one can show that the signs associated with the singlet always cancel out for systems with bipartite interactions) [54]. This way, the projector QMC method in practice becomes very similar to the finite-temperature SSE method [80], with the main difference being that the periodic imaginary-time boundary conditions in the SSE method are replaced by boundary conditions given by the trial state |Ψ(0) . The exact choise of this state is not critical, though a good variational state can improve the convergence rate in m significantly. The advantage of the projector approach relative to taking the limit T → 0 in SSE calculations is that the valence bonds restrict the system to the singlet sector (and other sectors can also be accessed by simple modifications). Thus, low-lyings S > 0 states that require very low temperatures to be filtered out in T > 0 calculations are excluded from the outset. For further technical details on the method we refer to the literature [54,79].
In the valence bond basis, expectation values are expressed using transition graphs [83,84] obtained by superimposing the bond configuration from the left and right projected states in (3). Spin-rotationally averaged quantities can be expressed using the loops of the transition graphs, e.g., the spin-spin correlation function between two sites i and j vanishes if the two sites are in different loops and is ±3/4 for sites in the same loop (with the plus and minus sign corresponding to sites on the same and different sublattices, respectively). Higherorder correlation functions involve more complicated expressions with the transition graph loops [85].

Order parameters and correlations
Here we will focus on the order parameters of the AFM and VBS phases. The former is the conventional sublat-tice (staggered) magnetization where the coordinates x i , y i ∈ {0, L − 1}. Since the simulations do not break the symmetry we evaluate the expectation of the squared order parameter, M 2 , which has a simple loop expression. The VBS order can form with horizontal or vertical bonds, and these are captured by the bond-order parameters where, for convenience, we have switched to a notation where the double subscripts on S x,y refer to the integer coordinates on the square lattice. In this case as well we need the squared order parameter, y , which has a reasonably simple direct transitiongraph loop estimator [85].
With the above order parameters we can also define the corresponding Binder cumulants. In the case of the O(3) symmetric AFM order the proper definition of the cumulant is where the coefficients are chosen such that, with increasing system size, U m → 1 in the AFM phase and U m → 0 if there is no AFM order. For M 4 as well there is a simple direct loop expression [85]. In the case of VBS order, the coefficients of the cumulant should be chosen as appropriate for a two-component U(1) symmetric vector order parameter, thus Here D 4 involves eight-spin correlation functions that in practice are too difficult to compute efficiently [85]. We therefore invoke an approximation in Eq. (7) that does not impact the scaling properties of the cumulant; we simply evaluate (D x , D y ) using the loop estimator for the two-point operators (5a) and (5b), and then use these classical numbers to evaluate D 2 and D 4 . While the expectation values entering in Eq. (7) are then not strictly the correct quantum-mechanical expectation values, they still reflect perfectly the absence or presence of VBS order in the system, and U D maintains the desired properties discussed above.
In addition to the squared order parameters M 2 and D 2 evaluated on the full lattice as described above, we will also consider the distance dependent spin and dimer Here open boundaries are used to avoid valence bonds crossing the boundaries. The thin red and blue arches correspond to valence bonds in the bra and ket state, respectively, while the thicker bonds represent the spinon strings that terminate at the unpaired spins (arrows up). Note that in each string (depicted with thicker lines) one of the unpaired spins is in the bra state and the other one is in the ket state.
correlation functions, where we spatially average over the reference coordinates (x, y) for each disorder sample. In the case of the spin correlations we will also consider the probability distribution of values without averaging over (x, y) or disorder realizations. The spin correlations have a staggered sign (−1) rx+ry , while the sign of the dimer correlator with x oriented bond as above is (−1) rx (and we take the proper average with the y-oriented ones). When presenting results we remove these signs. In C d (r) it is sometimes better to use the difference between even and odd distances instead of removing the squared mean value.

Spinon strings
In addition to the physical observables in the singlet sector discussed above, it is also useful to consider the lowest state with total spin S = 1, in which some aspects of spinons can be probed directly. In the valence-bond basis, an S = 1 state can be expressed with a "broken bond", e.g., with one bond replaced by two ↑ spins, one each on sublattice A and B (or with one bond treated as a triplet) [86][87][88]. These unpaired spins will propagate under the action of the Hamiltonian, and one can characterize their collective nature as bound or unbound, and, in the latter case, quantify the size of the bound state [68,88]. Here we will demonstrate another way to characterize an S = 1 state by simply using the number of sites involved in the spinon strings formed in the transition graphs. A string consists of the unpaired spin on a given sublattice in the projected bra and the ket state and a connected set of alternating bra and ket valence bonds, as illustrated in Fig. 6. As we will see in Sec. IV E, the mean number of sites in the strings scales very differently in the AFM and RS states, and this provides a way, along with other methods, to locate the phase transition between these two states. In addition, we will also use the difference in ground state energy between the S = 1 and S = 0 sectors to extract the spin gap. For technical details on how to carry out the simulations with broken valence bonds we refer to Refs. 81, 86-88.

B. Site Diluted J1-J2 static-dimer model
We begin our discussion of QMC results with a brief study of a statically dimerized system, where in the uniform system there is a quantum phase transition from an AFM to a trivial quantum paramagnet due to singlet formation at the stronger bonds. In the case of the columnar model illustrated in Fig. 2, the critical coupling ratio j 2c ≈ 1.91 [54,77,78]. For j 2 > j 2c , it is well known that effective S = 1/2 moments localize around diluted sites in such a system, and that these moments interact with each other by non-frustrated effective couplings mediated by the gapped host system [76], thus inducing AFM order also in the previously quantum-disordered phase. Here we use this system as a means of illustrating how this weak dilution induced AFM order is manifested in the quantities that we will later study in the more interesting models. For these illustrations we take the vacancy fraction p = 1/32, with a canonical ensemble such that exactly N/32 sites are removed, with equal numbers on the two sublattices. This density of vacancies is far below the classical percolation threshold, p c ≈ 0.407, beyond which no long-range order can exist. Fig. 7 shows results for both the squared sublattice magnetization and its Binder cumulant. The latter turns out to be a more sensitive quantity for detecting weak order. If there is a critical point separating the AFM phase from a non-AFM phase, the cumulants for two different system sizes, graphed versus the control parameter, should cross each other at a point that drifts toward the critical point with increasing L. However, as shown in Fig. 7(a), the crossing points in this case drift rapidly toward higher j 2 values and no convergence to a critical coupling can be found. In the inset of Fig. 7, the size dependence at two values of the coupling ratio deep inside the quantum paramagnet are shown. Here one can observe non-monotonic behaviors indicating asymptotic flows toward the value U M = 1 expected for long-range ordered AFM states. This behavior can be seen even though the order parameter itself, shown in Fig. 7(b), is very small. Here all the curves for different j 2 should extrapolate to M 2 > 0 when L → ∞, but for large j 2 the values are very small and not easy to extract precisely. With the behavior of the Binder cumulants, we can nevertheless confirm that there is long-range order at least up to j 2 = 5, and there is no reason to expect any other phase for still larger j 2 .
The reason for the decreasing AFM order with increasing coupling ratio j 2 deserves some discussion. This behavior can have more than one source and the most important should be: (i) The localized moments induce some AFM order in their vicinity and so each diluted site contributes effectively more than one unit of staggered magnetization. This effect decreases with increasing j 2 as the host becomes less susceptible to induced order. (ii) Some of the local moments will form singlets and do not contribute (or contribute very little) to the overall AFM ordering. This effect may also increase with increasing j 2 , as the effective interactions among moments at fixed distance becomes weaker and the distribution of couplings becomes broader. Therefore, some pairs of moments will become more specifically coupled to each other than to other more distant spins in their surroundings.

C. Site Diluted J-Q model
In the site diluted J-Q model, any J or Q term in Eq. (1) acting on one or more vacancies are excluded from the Hamiltonian. We consider small vacancy concentrations p. In the gapped VBS host, when Q > Q c , with Q c /J ≈ 1.50 [51], the vacancies are again expected to localize magnetic moments around them, and, by the same mechanism as in the J 1 -J 2 model studied above, one should expect weak AFM order to form among these. Thus, one would expect the sharp AFM-VBS transition to be ruined and the phase diagram to qualitatively be of the type in Fig. 5(a), with just a cross-over between strong and weak AFM order. The ground state for sufficiently large Q/J, in the local moment regime with weak AFM order, should still be different from the corresponding state in the J 1 -J 2 model, due to the formation of VBS domains with different angle φ in the notation of Fig. 3. This domain formation also should lead to localized spinons that are not induced directly due to sublattice imbalance at the vacancies, but form at the meeting points of four VBS domain walls. As we have discussed in Sec. III, these moments can be associated with quite different physics, but in the case of the diluted J-Q model this aspect of the problem should be masked at least partially by the more trivial moments forming at vacancies. We will here not study a possible interplay between the two different kinds of moments, but simply use the Binder cumulants to confirm the break-up of the VBS into domains and the formation of AFM order. Results for p = 1/32 at two different values of the coupling ratio are shown in Fig. 8. Here, in Fig. 8(a), we can again see, as we did in the case of the J 1 -J 2 model in Fig. 7, how the AFM Binder cumulant first decreases with increasing system size but then starts to grow when the number of moments becomes sufficient for AFM order to form. This cross-over occurs for larger sizes for the larger Q/J value, which is again similar to the behavior found for increasing coupling ratio J 2 /J 1 in the J 1 -J 2 model.

D. Random Q model
We next consider randomness in the Q interactions, with an extreme case of bimodal coupling distribution where each Q term in Eq. (1) is either absent or present (with equal probability). Here we take the strength the present six-spin couplings as 2Q, so that the parameter Q is the average six-spin coupling. As Q increases, the effective value of the disorder strength, Λ in Fig. 5(b), also increases when defined in relation to the constant J coupling. We will demonstrate a quantum phase transition between the AFM phase and the phase that we characterize as an RS phase as the coupling ratio Q/J increases. We will argue that the phase diagram is of the type schematically illustrated in Fig. 5(b), though we will not consider the full phase boundary versus Λ. We will demonstrate the existence of a quantum critical point separating the two phases along one path in parameter space and also characterize the ground state properties of the RS phase in various ways.

VBS domains and apparent lack of AFM order
First, in Fig. 9, we visualize the VBS domains forming in this kind of system for large Q/J, where the pure system is deep inside the VBS phase. Here we observe several instances of meeting points of four domain walls, where spinons are expected to be localized. Note that the static dimer pattern shown in Fig. 9, which is just a representation of the nearest-neighbor spin correlations, can be misleading due to the fact that it does not convey completely the quantum fluctuations. A thin line or the absence of any line on a given site implies large fluctuations of the associated spins, as further explained in the caption of Fig. 9, but the nature of those fluctuations is not apparent. Later, in Sec. VI, we will also visualize the local spin fluctuations and demonstrate that they are small within the bulk of VBS domains and large at regions corresponding to spinons and domain walls. Despite possible shortcomings of this type of visualization, it nevertheless makes clear the typical domain size and the manner in which domains meet. A notable feature is that there are mainly domain walls of the type where the angle φ (Fig. 3) changes by π/2, as would be expected according to the discussion in Sec. III. Some very short segments of π domain walls can also be seen, with a line of bonds oriented perpendicularly to those of the adjacent domains located in the gap between those domains. The π domain walls in a pure system with a two-fold degenerate VBS are gapless with deconfined spinons [82], and in a disordered system with a pinned π domain wall one can expect localized spinons to form pairwise as well. These spinons can also be regarded as meeting points of four domains, with two of the domains being extremely narrow (chain-like). Examples of local VBS patterns indicative of such spinons can also be seen in Fig. 9, in the form of π phase shifts between the VBS patterns of chain segments between two domains.
The main question now is whether AFM order is induced among the localized spinons that presumably exist in the random VBS environment. We again study the AFM Binder cumulant, Eq. (6), as a function of the Q interaction. For convenience, to span the full range of interactions, we graph U M versus Q/(J + Q) in Fig. 10(a). Interestingly, unlike in the diluted models (Figs. 7 and 8), in this case it appears that the cumulants for different system sizes develop a common crossing point as L increases; the standard signal of a quantum phase transition of the AFM state. Furthermore, as shown in Fig. 10, for values of Q/J larger than the apparent asymptotic crossing point, the cumulants decrease steadily toward zero and there are no indications of any upturn expected if the state has weak AFM order. One could of course wonder whether the turning point might occur only for even larger system sizes, but the very different behaviors of the crossing points between the diluted models, where they drift strongly as the system size increases (as shown in Fig. 7 in the case of the J 1 -J 2 model) suggests that the phase diagrams really are different.

Existence of a phase transition
The possibility of AFM order for large Q/J in the random Q model can be be excluded if we can convincingly establish the existence of a quantum critical point where the AFM order parameter and related quantities exhibit scaling. To this end, we will analyze the drift of the cumulant crossing points, and also consider an alternative way of locating the critical point.
As discussed in Sec. IV A, QMC simulations in the valence-bond basis allow also for studies of the lowest triplet state, which is associated with strings representing spinons in the sampled transition graphs (see Fig. 6). In an AFM state one can expect the spinon strings to cover a finite fraction of the system (and then the spinons are not well-defined particles [88]). We therefore define the string fraction λ as the mean fraction of sites covered by one of the spinon strings. In Fig. 11 we demonstrate that indeed λ approaches a constant when L increases inside the AFM phase, while in the RS phase λ ∝ L −1 . We do not have a rigorous explanation for the latter behavior, but it appears to be a very robust feature of the RS phase. Superficially, it would seem to indicate that the spinons are not completely localized but involve of the order of L spins. However, it should be noted that many spinons can be involved in forming the lowest total spin triplet, and the spinon strings will migrate during the simulations between all of them. The mean string frac- tion therefore is not really probing an individual localized spinon, and its physical meaning should be further investigated. Here we just exploit its apparent utility in locating the AFM-RS transition. Interestingly, as shown in Fig. 12(a), when graphed versus the coupling ratio, Lλ for different system sizes exhibits crossing points. This would not necessarily be expected when the behavior throughout the RS phase is λ ∼ L −1 , but is still possible due to scaling corrections; indeed, the fact that the crossings occur at decreasing angles when L increases and all the curves are close to each other for large coupling ratios suggest that corrections to the dominant power law are responsible. While the crossing point is still quite well defined and suggestive of a critical point, the weak size dependence inside the putative RS phase makes it hard to accurately extract the crossing points between curves for, e.g., system sizes L and 2L when L is large. Nevertheless, we have extracted a few crossing points and compare them with the crossing points extracted from Binder cumulant data such as those in Fig. 10. As shown in Fig. 12(b), in both cases the size dependence is consistent with a flow to a common point as L → ∞, with power law correction in 1/L. We do not have enough Lλ cross points to be able to do an accurate independent fit, but since we expect that both data sets scale to the same critical point we impose this condition in the fit shown in Fig. 12(b). The critical point so extracted is Q c /J ≈ 1.23. We take this analysis as strong evidence for a quantum critical point separating the AFM phase and a non-magnetic phase that we argue is an RS phase.

Correlation functions
Next, we consider the mean spin and dimer correlation functions. Fig. 13(a) shows the mean spin correlations, Eq. (8a), at the largest distance on the periodic lattices, r = L √ 2, versus the system size L. For three different coupling ratios inside the RS phase, we find the same behavior; a power-law decay corresponding to the distance dependence C s (r) ∝ r −α with α = 2. Instead of carrying out line fits to find α, we here just show comparisons with the form with α = 2, but individual fits in all cases are also consistent with this value. Interestingly, C(r) ∝ r −2 is also the form at the RS fixed point in 1D [36], though in that case there are apparently also multiplicative logarithmic corrections [42] that we do not find here in 2D. In the case of the dimer correlations defined in Eq. (8b), Fig. 13(b) shows results at the longest distance where we have extracted the relevant connected piece of C d (r) as the difference between even and odd distances r, which produces less noisy results than the method of subtracting the mean value in Eq. (8b). Here the relative error bars are still rather large for the larger systems, and we only show consistency with the form C d (r) ∝ r −4 , which again is the same form as in 1D (up to the log corrections found in 1D) [42].
It is also interesting to investigate the probability distribution of the values of the correlation functions in the spatially non-uniform system. Here we again consider the longest distance r ij = L √ 2 on the periodic square lattice and accumulate in histograms all the individual spin correlations C ij = C(r ij ) for spins at sites i, j separated by this distance, with a large number of disorder realizations used to produce reasonably smooth distributions. In this case it is important to run rather long simulations for each individual disorder realization, so that the statistical errors do not influence the distributions significantly for the smaller instances of C(r ij ) (in contrast to the mean disorder-averaged values, where one only has to make sure that the individual simulations are equilibrated and the final statistical error is dictated by the number of disorder instances). There will always be some problems with large relative errors for the smallest correlations, and therefore we expect the distributions presented below to be most reliable at the upper end of the distribution.  (9) has been set to a = 1/3, close to its optimal value for collapse of the data for the larger systems. The blue fitted curve on the left side of the distribution corresponds to the power-law behavior P (v) ∝ v n with n = 11. In (b) the scaling variable x = ln(Cij L 2 ) is used.
To investigate scaling of the distributions, we first attempt a scaling variable similar to one applicable to endto-end spin correlations of the random transverse-field Ising chain, which realizes an IRFP [89], and transform the histograms to the distribution P (v). In Ref. 89 the exponent a = 1/2, but here this does not work, and we therefore consider a as a fitting parameter. This indeed works quite well for the larger system sizes if a ≈ 1/3, as shown in Fig. 14(a). We also need the resulting data-collapsed distribution to be consistent with the mean correlation function, for which we previously found C ij ∝ L −2 . We can obtain a power law if the behavior of the probability distribution P (v) for the scaled variable v follows a power law close to 0; P (v) ∝ v n . It is easy to see that the contribution to the mean value from small v then decays as C ij ∝ L −a(n+1) , and with a = 1/3 we therefore need n = 5. The behavior in Fig. 14(a) is not consistent with this value of n, instead giving an exponent n more than twice as large (corresponding to C ij ∝ L −4 ), as shown with a fitted curve in the figure. However, the part of the distribution away from the region where the power law applies still changes the scaling of the mean value to the observed L −2 form for the rather small systems we have access to, for which e −vL a in Eq. (10) is not yet very small when v ≈ 2 ∼ 3. For large system sizes, the power law region would always dominate the integral and with the fitted form we would then obtain an L −4 decay. Since our data do not extend very close to v = 0 we can not exclude that the distribution still changes and evolves into the v 5 form as v → 0 and C ij ∝ L −2 .
Considering the apparent inconsistencies arising with the scaling variable v above, we exploring an alternative form of the distribution. Fig. 14(b) shows distributions P (x) with the scaling variable x defined as In this case any P (x) gives the desired L −2 decay of the mean. Though the data collapse is not as good as in Fig. 14(a), the behavior does seem to improve with increasing L, especially at the high end of x.
At the IRFP, a scaling variable of the form (9) and P (v) ∝ v n for small v implies different behaviors of the typical and mean correlations (exponentially versus power-law decaying) is a consequence of the divergent dynamic exponent [89]. As we will show in Sec. V, the RS state in the random J-Q model has finite dynamic exponent, and the scaling with the variable in Eq. (11), which implies the same power-law decay of the mean and typical values, may appear more plausible from this perspective. However, the scaling with the logarithmic variable in Fig. 14(a) works noticeable better and we cannot exclude that mean and typical values will scale differently even though z is finite. It would clearly be useful to study larger system sizes and further test the two scenario for the distributions. The inverse-square distance dependence of the mean correlations already appears to be well-established by the good scaling for a wide range of system sizes and three different Q/J values in Fig. 13.

E. Random J model
In the random J model, all Q couplings are included and the J couplings are drawn from a distribution. We have consider bimodal as well as continuous distributions and find qualitatively the same kind of behaviors as above in the random Q model. We therefore only provide a few illustrative results showing these similarities. Fig. 15 shows results for the order parameters and Binder cumulants at Q/J = 2 for the extreme bimodal case where half of the J couplings are set to 0 and the rest to 1 (which we here take as the value of J in the ratio Q/J). For reference we compare the size dependence of these quantities with the corresponding pure system (all J = 1). The results indicate that both order parameters vanish when L → ∞, with the VBS Binder cumulant again showing a non-monotonic behavior with a drop toward zero starting when L is of the order of the typical VBS domain size. For Q/J = 2 we conclude that the system is in the RS phase.
To confirm the existence of a critical point separating the AFM and RS phases, Fig. 16(a) shows scans for several system sizes of the Binder cumulants versus Q/J for the same bimodal J distribution. For U M we again see crossing points apparently converging toward a critical point, similar to the behavior in the random Q case in Fig. 10. The (L, 2L) crossing points are graphed versus the inverse system sizes in Fig. 16(b), along with the crossing points of the scaled string fraction Lλ. These two finite-size estimates of the critical point again approach Q c from different directions. The data here are of better statistical quality than the random Q data in Fig. 10, and with both data sets we can fit the corrections to the infinite-size critical point to the expected forms ∝ L −ω . Requiring the fits to have the same value of Q c /J but allowing for different values of the exponent ω, we obtain Q c /J ≈ 0.72 and the exponents ω ≈ 1.5 (for the cumulant crossings) and 2.3 (for the string quantity). Given the rather small number of points and not very large system sizes, the exponents should be regarded as "effective exponents" that are still influenced by neglected higher-order corrections.  Fig. 12) are also shown. Fits (the curves shown) to the latter data set and that for the UM crossing points were carried out using power-law corrections, ∝ L −ω (with ω ≈ 1.5 and 2.3 for the UM and Lλ set, respectively), with the constraint of the same value of the crossing point Q * /J when L → ∞. Fig. 16(a) also shows the behavior of the VBS cumulants. It is clear that the crossing points here do not converge but flow to larger Q/J as the system size increases, as would be expected when arbitrary weak disorder destroys the VBS phase. The (L, 2L) crossing points are graphed versus 1/L in Fig. 16(b).
Overall, with the results presented above and many other cases, we find very similar behaviors for the random Q and random J models, indicating that the RS phase induced by these types of disorder is the same one. One notable aspect of the specific random J model for which we have presented results here is that the RS phase can arise not only out of the VBS phase of the pure model but also from the AFM state. The critical coupling extracted in Fig. 16 is at Q/J ≈ 0.72, where the pure model with all J = 1 Heisenberg couplings is still well inside the AFM phase (the AFM-VBS transition of the pure system taking place at Q/J ≈ 1.50). With the way we have defined the bimodal coupling strengths with J = 0 and J = 1 at random locations, we can reach the RS from the AFM phase simply by removing some fraction of the J interactions when Q is between 0.72 and 1.50. This random removal of J couplings enhances the ability of the Q terms to cause VBS formation, which in the ran- dom system only can take the form of a domain-forming VBS. Thus, it seems very plausible that the same RS state will also be generated if the host system includes some frustrated interactions that weaken the AFM order and favor local formation of VBS domains in a disordered system, instead of the Q terms considered for that purpose here. Such frustrated disordered systems can include the Heisenberg model on the triangular lattice, which is equivalent to the square lattice with half of the diagonal couplings activated. It would then appear quite plausible that RS state we have identified here on the square lattice is actually the same state as that discussed previously for frustrated systems.

F. Universality of the AFM-RS transition
Given our results presented above, it appears most likely that the AFM-RS transition is universal and that the RS phase itself has universal properties, such as the 1/r 2 power-law decay of the mean spin correlations (but we will show in Sec. V that the dynamic exponent is not universal inside the RS phase but varies continuouslythough it also is universal at the AFM-RS transition). An often used characteristic of a critical point is the value of the Binder cumulant. This quantity is universal, in the sense that it is independent on microscopic details, but, unlike many other universal quantities, such as critical exponents, it depends on boundary conditions and aspect ratios of the system [90][91][92]. In the projector QMC method we effectively take the limit of the time-space aspect ratio β/L → ∞ and the system geometry is also the same for both the random Q and random J models. Thus, we have identical boundary conditions and aspect ratios, and would expect the same value of the Binder cumulant at the AFM-RS transition point.
In Fig. 17 we show results for three disorder types for which we have sufficient data to carry out careful studies of the scaling of the AFM cumulant at the (L, 2L) cross-ing points; in addition to the bimodal Q and J cases we also show results for a continuous distributions of J, with values drawn uniformly from the range [0, 2]. Remarkably, the cumulants for all cases not only appear to flow to the same point in the limit of infinite size, but even the leading correction in 1/L seems to be the same (even as regards the prefactor of the power-law correction). This correction appears to be almost linear, and we analyze the data under this assumption, though it is more likely that the form is L −ω with ω just close to 1. For the two bimodal distributions all the data fall on the line as closely as would be statistically expected (with excellent goodness of fit), while for the continuous distribution we see that the data for the smaller sizes deviate more significantly, indicating that the higher-order corrections do depend on the kind of disorder distribution. These results clearly lend further support to the existence of a universal AFM-RS critical point, and, therefore, to the existence of the RS phase.

V. FINITE TEMPERATURE PROPERTIES AND THE DYNAMIC EXPONENT
Finite temperature properties are useful for extracting the dynamic exponent z and may be the most direct route to connect to experiments. We will here consider the uniform magnetic susceptibility, and the local susceptibility defined by the Kubo integral where S z r (τ ) is the standard imaginary time-dependent spin accessible in QMC simulations. We here use the SSE method and refer to the literature, e.g., Ref. 54, for further technical information. In this section, we average the local susceptibility over all the sites r of the system (as well as over disorder realizations) and call this averaged quantity χ loc . In Sec. VI we will show an example of the spatial dependence of χ loc (r) for a fixed disorder realization.

A. Power-law behaviors
At a quantum critical point of a system such as those considered in this work, where the magnetization is a conserved quantity, the susceptibility should scale with the temperature as [93] where d = 2 in our case. In contrast, the local susceptibility is sensitive to the fluctuations of the non-conserved Here β/ν should be equal to κ/2, where κ = 2 is the exponent we have found for the decay of the spin correlations; C s (r) ∝ r −κ . Thus, we expect the asymptotic form χ loc ∝ T 1/z−1 , which diverges faster than the uniform susceptibility Eq. (14). We also note that, in the alternative (less likely) scenario where κ = 4 (Sec. IV D 3), we would have χ loc ∝ χ u . For the above forms to be valid, we not only have to reach sufficiently low in T , but also the system size has to reach the range where there is no longer any size dependence left. This requirement limits the temperatures we can reach, as demonstrated in Fig. 18 in the case of the uniform susceptibility of the random Q model close to the critical point and inside the RS phase. We can still see critical behaviors emerging for a range of temperatures for the largest system sizes. In Fig. 18(a), at Q/J = 1.25, which should be close to the AFM-RS transition according to the results in Fig. 12, we find very little size dependence, indicating, by Eq. (14), that z = d = 2 at the transition. The small increase seen at low T before the finite-size form is most likely due to Q/J not being exactly at the AFM-RS transition but slightly inside the RS phase. Well inside the RS phase, at Q/J = 2 as shown in Fig. 18(b), we find a clearly divergent low-T behavior of χ u . Since the overall magnitude of the susceptibility originating from the localized spinons is small, when fitting to the expected power-law form we also include a constant, as a natural leading correction to the asymptotic divergent form. This works well and the value of the exponent given by the fit corresponds to z ≈ 5.6. Thus, we find that z increases rapidly as the RS phase is entered. Figure 19 shows results even further inside the RS phase, along with fits such as those discussed above. Here we find z ≈ 6.2 at Q/J = 4 and z ≈ 8.3 when Q/J → ∞ (J = 0). In the latter case it should be noted the bimodal disorder distribution, where half the Q couplings are set to zero, can lead to isolated spins that contribute ∝ 1/T to the susceptibility. However, we avoid this issue by "patching" such rare isolated spins by adding a randomly chosen Q interaction for each of those spins to connect them to the rest of the system. Figure 19 also shows results for the local susceptibility. Here we do not observe the expected faster divergence than in χ u , given by Eq. (15), and it appears that the asymptotic temperature regime has not yet been reached. This could be because the local susceptibility only contains a small fraction of the dominant staggered response, q = (π, π) in momentum space, and therefore one may expect large corrections from all the other momenta at which the response is weaker. We conclude that lower temperatures would be required in order to see the behavior predicted by Eq. (15). An easier way to detect the dominant dynamic response, but that we have not yet pursued, would be to compute the susceptibility in momentum space at q = (π, π).

B. Griffiths-McCoy singularities
To properly classify the proposed RS state, we need to consider the fact that disordered systems generically have regions in parameter space called Griffiths, or Griffiths-McCoy, phases. These phases or regions are characterized by spatial 'commingling' of two phases [94,95]. Fluctuations in the quenched disorder can favor a phase B within a limited part of a system that is overall in a phase A. Griffiths phases, which do not always have well-understood RG fixed-point analogues (but sometimes they do [96]), appear close to critical points and are normally associated with weaker singularities than the actual critical points (for reviews, see Refs. 24 and 25). The singularities arise from exponentially rare regions (e.g., large domains of phase B inside phase A) and have the most profound effects on dynamical properties.
In quantum systems, Griffiths phases typically have large but finite dynamic exponents, with associated divergent susceptibilities if z > d. The large z values (long time scales) motivates the often used term "glass" for these phases, though a Griffiths phase is not normally associated with the multitude of thermodynamic states (by replica symmetry breaking and related phenomena) of classical and quantum spin glasses (and it was also claimed that the VBG state does undergo replica symmetry breaking [72] but this may be a consequence of a classical treatment). Examples of Griffiths phases include the Bose glass in the Bose-Hubbard model with random potentials [30] and the Mott glass in particle-hole symmetric boson systems where randomness is introduced in the hopping constants (and there are indications that this state can also form with random potentials due to emergent particle-hole symmetry [97]). The spin analogue of particle-hole symmetry is also present in 2D random exchange Heisenberg antiferromagnets, where Mott-glass phases have been identified [98,99].
An important question is whether the RS state we identify in the random J-Q model is also a Griffiths phase. We argue that it is not, because equal-time correlations in Griffiths phases normally decay exponentially with distance (a fundamental consequence of the rareregion mechanism), while we find strong evidence for power-law correlations.
There is a further strong argument against the RS phase being a Griffiths phase: If, in the language above, we consider the AFM as phase A, there is no obvious phase B with which A can commingle to form the RS phase as a Griffiths phase. The RS phase is then actually that phase B, and in principle Griffiths singularities could appear due to comingling of the AFM and RS phases close to the phase boundary. However, since the AFM and RS phases are both gapless, the Griffiths singularities would be very hard to detect and would very unlikely be responsible for the power laws we have identified here. Most likely, they would only cause scaling corrections and no separately identifiable Griffiths phase in addition to the AFM and RS phases.

VI. SPINON INTERACTION MECHANISM
The way the localized spinons interact with each other is a crucial ingredient in the formation of the RS state. In order for singlets to be gradually "frozen out" as the energy scale is reduced (as in the Ma-Dasgupta strongcoupling RG procedure [35,43]), and for AFM order not to form on large length scales, it is necessary that the interactions are not completely random. The fact that spinons are created in pairs (spinons and antispinons) when VBS domains are formed already implies a correlation that favors closer typical distance between a spinon and the nearest antispinon, because the domain walls will provide an effective potential due to domain-wall energy between a spinon and anti-spinon site connected by a wall. There is, however, potentially also another effect, namely, the effective magnetic interactions between the spinons are likely mediated mainly by the domain walls. The putative role of domain walls as mediators of spin correlations was mentioned in Ref. 17 but was not developed into an actual mechanism suppressing AFM order and causing the singlet formation in the RS state. Here we will provide evidence for this mechanism within our models on the square lattice. We note that the effective interactions should have the same bipartite nature as the microscopic interactions in the pure system, as was discussed in a generic situation in Ref. [17].
First, let us consider a domain wall in the pure J-Q model in its VBS state. According to the DQC theory [69], the thickness of a domain wall between VBS domains, across which the angle φ defined in Fig. 3 changes by π/2, is not governed by the standard correlation length ξ, but by a longer length scale ξ ′ (i.e., diverging faster than ξ as the DQC point is approached). This affects the scaling of the energy density of the domain wall as the critical point is approached [70,81], which may also have a counterpart at the AFM-RS transition. We here only mention this and do not explore the domain wall thickness further. Instead we discuss the spin gap of a domain wall, i.e., the energy difference between the S = 0 ground state and the lowest S = 1 state in a system with a domain wall imposed by boundary conditions. Figure 20 shows an example of a domain wall, where the bond thickness on a 32×32 lattice corresponds to the magnitude of the spin correlation on that bond, and the colors of the bonds are coded as in Fig. 3. Here the boundary conditions are periodic in the vertical direction but in the horizontal direction the interactions have been modified (see Ref. 81) so that the edges are locked into VBS realizations differing by the angle φ = π/2 of an elementary domain wall. Here it should be noted that the length scale over which the angle φ changes in Fig. 20 is not the intrinsic domain wall width, because the location of the wall also has quantum fluctuations that smear it out when expectation values are computed. The spin gap of the wall is still a completely well defined quantity, as long as the S = 1 excitation (observed, e.g., with the spinon strings illustrated in Fig. 6) is not repelled from the wall. We have confirmed that the excited spin is attracted to the domain wall.
The spin gap is obtained by simply taking the difference between total ground state energies computed in the two spin sectors. Fig. 21 shows results for the uniform system without domain wall (obtained with fully periodic L×L lattices) and with domain walls on lattices with two different aspect ratios, as a check for the expected independence on the lattice geometry when L → ∞. For small systems with a wall, the gap is strongly influenced by the boundary modifications, which here extend three rows into the system on each side, and one should not draw any conclusions on the differences between the system with and without the domain wall until L is much larger and the actual domain wall has converged to its intrinsic thickness. For large L, it is clear that the gap on the domain wall is significantly smaller than in the bulk, as one might have expected just from the fact that the domain wall has weaker order, i.e., more fluctuations, than the bulk VBS. Thus, in a non-random system, a domain wall will be a more effective mediator of correlations, and thereby of effective interactions between impurity spins, than the bulk VBS. The above results for a pure infinitely long domain wall should only be taken as suggestive of enhanced spinon interactions along domain walls in the disordered system. We can obtain further evidence by examining the spatial variations of the local susceptibility, Eq. (13), for individual disorder realizations (see Ref. [87] for similar calculations for a diluted Heisenberg system). A large susceptibility can be taken as a sign of a small local gap, through the sum rule (here at T = 0) where S loc (r, ω) is the local dynamic spin structure factor, which satisfies the sum rule For any finite system, the spectral weight in S loc (r, ω) does not extend all the way down to ω 0 , and in a singlemode approximation, where there is only a single δfunction at ω = ∆, we can extract the 'local gap' as ∆(r) = 2S loc (r)/χ loc (r) = 1/2χ loc (r). In the realistic case where there is a broader distribution of spectral weight, χ loc (r) can still be regarded as a proxy for the typical local low-energy scale, and it should then also be a measure of the local ability of a region of the system to mediate spin-spin interactions.
In Fig. 22 we show the spatial dependence of the local susceptibility for the same Q disorder realization as in the illustration of VBS domains in Fig. 9. Several bright spots on the susceptibility map can be observed, and some of them can be matched with meeting points of four  Fig. 4. The values of the susceptibility defined in Eq. (13) have been rescaled so that the maximum is 1, and the color coding is shown on the bar on the right side.
VBS domain walls, where spinons should localize. Naturally the sites on which the spinons reside should have enhanced susceptibility (and note that a single spinon will be spread out over several sites due to quantum fluctuations). There are also bright regions in Fig. 22 that cannot be specifically identified as likely spinon locations in Fig. 9, showing that also other VBS defects can be associated with small local gaps. It is not possible in this kind of picture to accurately mark out all the domain walls based on the susceptibility map-for this to work well we should go to a currently intractable limit where the domains are much larger (thus, requiring large lattices). However, it is clear from comparing Figs. 9 and 22, that the susceptibility is low in the interiors of large domains, and that in turn means that the susceptibility of the domain walls is enhanced over that of the bulk. This supports the notion that the domain walls act as mediators of spinon-spinon interactions, which should play an important role in the formation of the RS state.
As mentioned in Ref. 17, although the SDRG procedure [35,36,43] on the 2D square lattice normally flows away from infinite randomness, one cannot exclude a flow to a finite-randomness fixed point. Such fixed points have been obtained in SDRG calculations on various bipartite and frustrated 2D Heisenberg models with disorder [32,33], but the physical properties of those fixed points do not appear to correspond to the RS phase discovered here. The key physical ingredients underlying the RS phase-VBS domains and localized spinons-are unlikely generated in the SDRG procedure applied to bipartite Heisenberg models, and the 'cluster states' generated in the presence of frustration also appear to be quite different. It is furthermore very difficult to apply the SDRG approach to more complicated interactions like the sixspin Q terms used here (which are very difficult to deal with even in 1D systems, though that has been done [42]) since many kinds of effective couplings can be generated during the decimation process.
It may be more fruitful to consider an SDRG process carried out only on the subsystem of the localized spinons. Though we have not actually carried out such an RG procedure and it is not clear how to actually formulate the spinon subsystem quantitatively, an intuitive picture follows from our results and arguments. Because of the formation of groups of even number of spinons (half of which are anti-spinons), an explicit spatial correlation exists between spinons and antispinons that will automatically lead to a reduction of statistical sublattice imbalance upon course graining, i.e., within a region of finite length l the relative difference between the number of spinons and anti-spinons should decrease much faster as l increases when compared to the case of randomly located spinons with Gaussian fluctuations in the sublattice imbalance. This statistical effect in combination with the domain-wall mediated spin interactions, should lead to singlets gradually "freezing out" one-by-one without first forming larger moments due to sublattice imbalance. This picture is very similar to the flow of the SDRG in the 1D random Heisenberg chain [35,36], which, however, leads to z = ∞ instead of the continuously varying finite dynamic exponent found here.

VII. CONCLUSIONS AND DISCUSSION
Using the J-Q model, we have demonstrated that an RS phase can be induced by disorder in a quantum spin system even though all microscopic interactions are bipartite, lacking the geometric frustration that so far was believed to be a necessary ingredient for this type of 2D state. It appears most likely that the state is the same one, in the RG sense, as those previously identified in frustrated systems [13][14][15][16][17][18][19], though the lack of definitive quantitative results of the previous works for, e.g., various exponents governing power-law behaviors, makes it difficult to definitely ascertain this at the moment. For example, it was argued that the low-T susceptibility follows a Curie form in the frustrated honeycomb Heisenberg model in the RS phase [16], while we have here demonstrated a T −a behavior with varying a < 1 in the random J-Q model (and a → 0 as the AFM phase is approached). Clearly ED studies of lattices with only up to ≈ 20 sites cannot be used to reliably address the detailed form of the divergence, as we have seen even with much larger systems here. Also in the theory of Kimchi et al. [17] it was not possible to obtain quantitative values of most of the exponents pertaining to the RS phase in the triangular lattice, though we note a more recent work in which scaling forms for the heat capacity (which we have not yet investigated) were obtained under various conditions and compared with experiments [18]. We note in particular that the previous works have not discussed any details of the AFM-RS phase transitions, for which we have obtained specific results here on power laws both at T = 0 and T > 0.
Given some of they key results that we have obtained here, such as the r −2 form of the decay of the mean spin correlation functions and the temperature independent magnetic susceptibility at the AFM-RS transition (and the divergent behavior with varying exponent inside the RS phase), targeted calculations aiming at these specific universal characteristics can hopefully soon be carried out also for the frustrated models. One promising calculational route here is tensor network states tailored specifically to disordered spin models [100,101]. Though such calculations are certainly challenging, it may still be possible to reach larger system sizes than in the previous exact ED and DMRG studies.
If indeed the universality of the RS phase encompasses the J-Q model as well as the multitude of frustrated quantum magnets, the ability to study the former with large-scale unbiased QMC simulations has significant consequences in the context of experiments. It will then be possible to relate observed power laws directly to unbiased calculations, e.g., to test relationships between the power laws for different physical observables. Although the J-Q model does not represent the correct microscopic interactions of specific materials, its phases can still contain the experimentally relevant low-energy physics. This is in the spirit of "designer Hamiltonians" [102], which are tailored to realize collective quantum states and quantum phase transitions while at the same time being amenable to numerical calculations, especially sign free QMC simulations, on large scales without approximations.
A specific experimental question would be to study the AFM-RS transition. A recent candidate system in which this may be possible is the effectively square-lattice (quasi-2D) material Sr 2 Cu(Te x W 1−x )O 6 , which so far has been synthesized at x = 0.5 [22]. The corresponding isostructural compounds Sr 2 CuTeO 6 and Sr 2 CuWO 6 have dominant nearest-and next-nearest-neighbor spin interactions, respectively, due to the different orbital properties of the plaquette centered Te and W ions. With random distribution of these ions, it was argued that an RS type state may form. By tuning x it should then be possible to study the AFM-RS transition as well, and test, e.g., our prediction of the T independent susceptibility at the transition and the divergent behavior as the RS phase is entered.
Many interesting QMC calculations are called for as extensions of the initial study of the random J-Q model presented here. For example, the evolution of the RS phase as a function of an external magnetic field (which was recently studied in the case of the triangular lattice in Ref. 18) is very interesting theoretically and also from the experimental perspective. The field can also can be included with the SSE method (as recently done in the case of the J-Q chain [103]). The dynamical signatures, e.g., the dynamic spin structure factor, can also be studied using SSE supplemented by analytic continuation techniques [104], and it will be interesting to compare the 2D RS phase with the case of the RS state in the random exchange Heisenberg chain, which was also recently studied with the above mentioned techniques [105].