Broken-Symmetry Ground States of the Heisenberg model on the Pyrochlore Lattice

The spin-1/2 Heisenberg model on the pyrochlore lattice is an iconic frustrated three-dimensional spin system with a rich phase diagram. Besides hosting several ordered phases, the model is debated to possess a spin-liquid ground state when only nearest-neighbor antiferromagnetic interactions are present. Here, we contest this hypothesis with an extensive numerical investigation using both exact diagonalization and complementary variational techniques. Specifically, we employ a RVB-like many-variable Monte Carlo ansatz and convolutional neural network quantum states for (variational) calculations with up to $4\times 4^3$ and $4 \times 3^3$ spins, respectively. We demonstrate that these techniques yield consistent results, allowing for reliable extrapolations to the thermodynamic limit. Our main results are (1) the determination of the phase transition between the putative spin-liquid phase and the neighboring magnetically ordered phase and (2) a careful characterization of the ground state in terms of symmetry-breaking tendencies. We find clear indications of spontaneously broken inversion and rotational symmetry, calling the scenario of a featureless quantum spin-liquid into question. Our work showcases how many-variable variational techniques can be used to make progress in answering challenging questions about three-dimensional frustrated quantum magnets.


I. INTRODUCTION
Featureless ground states of interacting quantum spins with exotic properties and emergent excitations are highly sought after. Identifying such quantum spin liquids (QSL) is particularly challenging in threedimensional (3D) systems due to fast scaling of the Hilbert space with linear system size. A key ingredient that favors liquid over ordered phases is geometrical frustration, with the pyrochlore lattice a prominent frustrated 3D lattice. On the pyrochlore lattice, even the (classical) antiferromagnetic spin-1/2 Ising model has a ground-state manifold with an extensive degeneracy, governed by the "two-in, two-out" spin ice rule, instead of magnetic order [1][2][3][4]. Upon inclusion of a small exchange term, excitations with fractionalized "magnetic" charges and gauge photons emerge, the hallmarks of a U (1) QSL [5][6][7][8][9][10].
At the same time, the nature of the ground state away from perturbative limits around the Ising point on the pyrochlore lattice is subject to long-standing debates. A prominent example is the ground state of the SU (2) Heisenberg nearest-neighbor antiferromagnet: Perturbative or mean-field treatments suggest the ground state to be dimerized [11][12][13][14][15][16][17][18][19][20] or to possess chiral magnetic order [21]. Alternatively, the aforementioned U (1) QSL phase might continue from the small exchange interaction region to the SU (2)-symmetric point without a phase * nikita.astrakhantsev@physik.uzh.ch † tom.westerhout@ru.nl transition [8,10,[22][23][24]. Resolving this debate is a challenging theoretical problem with direct relevance to real compounds, such as rare earth molybdenum oxynitride pyrochlores R 2 Mo 2 O 5 N 2 , which are expected to admit an antiferromagnetic isotropic spin-1/2 Heisenberg model description [25]. Remarkably, the recently synthesized Lu 2 Mo 2 O 5 N 2 experimentally shows neither sign of magnetic order nor spin freezing [26,27]. A QSL phase, if present, is expected to be close in parameter space to several symmetry-breaking ordered phases [28][29][30][31][32]. These ordered phases can be induced by including next-to-nearest neighbor couplings [22,28] or large transverse exchange interactions, with the latter stabilizing spin-nematic order [8,9,23,24]. This suggests a close competition between nearly degenerate QSL and ordered phases.
In this work, we assess whether the spin-1/2 pyrochlore quantum antiferromagnet hosts a QSL or a symmetrybroken ground state [38] in the vicinity of the SU (2)symmetric point and further identify the adjacent phases. For this we use state-of-the-art variational Monte Carlo (VMC) methods. We use two complementary and highly flexible ansatz wave functions, thus controlling the inherent parametrical bias. On the one hand, we use the many-variable variational Monte Carlo (mVMC) method inspired by the resonating-valence-bond (RVB) wave function [39], and recently adopted for two-dimensional frustrated models [40]. The second ansatz uses neural network quantum states (NQS) [41], which have proven to be a powerful addition to the numerical toolbox for many-body quantum systems [40,42]. Finally, ED calculations (we employ a novel package SpinED [74]) are used to benchmark both approaches and guide our analysis.
Our methodology allows us to obtain reliable wave functions on clusters as large as 256 = 4×4 3 sites and perform extrapolations to the thermodynamic limit, which is notoriously hard for a 3D frustrated system. Figure 1 summarizes our main results in the form of (a) the phase diagram and (b) our energy extrapolation to the thermodynamic limit. In the putative QSL phase, we observe clear signatures suggestive of spontaneous breaking of both inversion and rotation symmetry. Furthermore, we extrapolate the phase transition of this nonmagnetic symmetry broken phase with the adjacent magnetic phase upon addition of next-nearest neighbor coupling j 2 . We find the phase transition at 7 times smaller j 2 /j 1 as compared to previous FRG results [22], which expands the magnetically ordered phase significantly.
The rest of this paper is organized as follows: In Sec. II, we introduce the model and our methodology. In Sec. III, we discuss our results, starting with a validation of the method by combining results from ED with the two variational approaches and a characterization of the magnetically ordered phase. This is followed by a thorough analysis of the symmetry-breaking signatures in the nonmagnetic phase. Finally, Sec. IV concludes with a discussion.

II. MODEL AND METHODOLOGY
A. Model and system geometry We consider the interacting spin Hamiltonian H(λ, j 1 , j 2 ) = j 1 i,j ĥ λ (i, j) + j 2 i,j ĥ λ (i, j), (1)  The non-magnetic phase (green) in the vicinity of (1, 0) is separated by a phase transition (yellow) from the magnetically ordered k = 0 phase (red). The U (1)π QSL phase is shown in grey. Further indicated are the symmetry-breaking properties of the respective phases with respect to inversion (Î) and rotation (R) symmetry. The ED phase boundary was obtained for a 4 × 2 3 system. Blue data points with errorbars show our infinite-volume phase-boundary extrapolations at λ = 0.8 and 1.0. Note that our data does not allow us to reliably draw conclusion about the shape of the phase boundary in the region λ 1. The shape of the yellow region, shown in Fig. 1 (a) at λ 0.8 is a guide to the eye consistent with perturbative results for small λ. The inset in the upper left shows the pyrochlore structure with sublattice indices. (b) The best variational energies obtained within this work using mVMC on equilateral (EC) and non-equilateral (NEC) clusters at λ = 1, j2/j1 = 0. The dashed line represents the infinite-volume extrapolation. For comparison, the figure also shows the DMRG result of Ref. [38] without bonddimension extrapolation. The errorbars at Ω −1 = 0 represent the infinite-volume extrapolations obtained within this work and DMRG bond-dimension extrapolation of Ref. [38], see main text for details.
where . . . denotes summation over nearest neighbors, . . . over next-to-nearest neighbors on the pyrochlore lattice and the local quantization axis of spin points towards the center of a tetrahedron. In this work, we are mainly interested in determining the (λ, j 2 /j 1 ) phase di-agram. Two limiting cases of the phase diagram are readily identified. On the one hand, the point (λ = 0, j 2 /j 1 = 0) corresponds to the classical Ising model which hosts the spin-ice ground state manifold. On the other hand, λ = 1 corresponds to the SU (2) symmetric Heisenberg model. Fragments of this phase diagram, which is schematically depicted in Fig. 1 (a), were previously studied. The classical Ising model orders magnetically upon inclusion of a small j 2 /j 1 term [22,[29][30][31][32]. The so-called k = 0 magnetic order does not break translation symmetry, while the magnetic moments within each of the four sublattices are ordered ferromagnetically, such that the total magnetization within each tetrahedron vanishes. On the j 2 /j 1 = 0 axis, perturbation theory in λ > 0 reveals the U (1) π QSL phase [6,43,44].
In order to better clarify the symmetries of the problem, we start by discussing the geometry of the pyrochlore lattice. The inset in Fig. 1 (a) shows the pyrochlore crystal structure. The pyrochlore lattice Λ is the union of four sublattices Λ α , α = 0, 1, 2, 3, which in Cartesian coordinates are spanned by the unit vectors e 1 = (1, 1, 0) T , e 2 = (0, 1, 1) T , e 3 = (1, 0, 1) T , Here, e 0 = (0, 0, 0) T and L 1 , L 2 , L 3 are the number of unit cells in the three respective directions that span the lattice. In each unit cell, the resulting four sites form a tetrahedron. Note that each site of this lattice is shared between exactly two tetrahedra. In what follows, we refer to "up-tetrahedra" and "down-tetrahedra" as indicated in blue and orange, respectively in Fig. 1 (a). The total number of sites is denoted by Ω throughout the manuscript. We consider pyrochlore clusters with periodic boundary conditions. Their number of bonds is 12 times the number of unit cells. We associate with one unit cell the four lattice sites labelled 0, 1, 2, 3 in Fig. 1 (a) as well as the 12 bonds comprising the up-and down-tetrahedra that have their corners labelled. The space group of the pyrochlore lattice is F d3m with a point group isomophic to O h . In the following, to study symmetry breaking tendencies, we will mainly focus on equilateral pyrochlore clusters (L 1 = L 2 = L 3 ), which have a point group D 3d . It is generated by C 3 rotations around the "easy-axis" (1, 1, 1) T , which cyclically exchange the sites (1 → 2 → 3 → 1) [see Fig. 1 (a)], inversion symmetry r → −r, which also exchanges down-and up-tetrahedra and a mirror symmetry, which reflects the cluster with respect to the plane passing through the bond 23 and the middle of the bond 01 [45]. Importantly, such equilateral clusters avoid geometric bias for the symmetries studied in our calculations.

B. Many-variable wave function
The many-variable variational Monte Carlo (mVMC) method has proven successful in studies of strongly correlated phases [46,47], including the QSL and valence bond solid (VBS) in the j 1 -j 2 Heisenberg model on the square lattice [40,48]. Here, we use the highly-optimized realization of mVMC from Ref. [39,49]. At the heart of the mVMC method is a mapping of spin operators to fermionic bilinearŝ where i labels a lattice site, a = x, y, z, and σ x , σ y , σ z are the three Pauli matrices. The RVB-like pairing state has the form where single occupation is ensured by the P ∞ G singleoccupation Gutzwiller projector. Note that the P ∞ Gprojected fermionic Hilbert space can be mapped to the original Hilbert space of spin operators. The wavefunction value σ|φ pair of a specific spin configuration |σ is evaluated using the Slater determinant of the matrix with elements f i,j . Here, σ represents a string of ±1, which, for each lattice site, stands for the respective spin eigenstate in the S z basis. The parameters f i,j are optimised using the stochastic reconfiguration optimization technique [50], which can also be seen as a way of performing stochastic imaginary-time dynamics in the variational manifold [41,51].
The general idea of the NQS method is to use the spin configuration σ as input for a neural network Ψ, interpreting the result Ψ(σ) as the (not normalized) wave function component, corresponding to the basis vector σ. To cope with the large system size and possible overfitting [52], we employ the convolutional neural network architecture (CNN) with real parameters and an elaborate alternating training technique (see Appendix A for details). This choice of architecture allows for better optimization and avoids certain instabilities [64]. We also benchmark it against the Restricted Boltzmann Machine (RBM) network, comprising one fully-connected dense layer, which is the classic NQS architecture [41,55]. The details of the NQS architectures and the variational parameter training can be found in Appendix A.

D. Symmetry-projected wave functions
To obtain highly accurate variational wave functions and energies, we employ quantum-number projections, i.e., we impose the ansatz state to transform in a chosen irreducible representation of the symmetry group.
Any point-group symmetryĜ is projected by applying it until the symmetry orbit is exhausted where ξ is the desired projection quantum number and |Ψ ξ the symmetrized state. Similarly, within mVMC, the projection onto the total spin S is performed by superposing the SU (2)-rotated wave functions [49]. At λ = 1, when only U z (1) spin rotation symmetry is present, the symmetry is enforced by only working in the space of total-zero magnetization M z |ψ = 0. In the NQS method, the spin-parity projector is applied: where Ψ(±σ) is the wave function evaluated at spin configurations σ and the spin configuration flipped along the z-axis, −σ. Such a projector selects wave functions of either even or odd total spin. For the mVMC ansatz, momentum and point-groupsymmetry projection can partially be performed by directly constraining the variational parameters f i, j . Otherwise, taking all Ω 2 parameters f i, j independent and symmetrizing the wave function at the end leads to a prohibitively large number of terms in the projector. In addition to computational cost, this makes the optimization procedure prone to false minima convergence. As a compromise, for clusters larger than 48 sites, we impose translational symmetry on the variational parameters f i, j and project the other symmetries using Eq. (6). Similarly, within the NQS method, the CNN architecture automatically imposes translational symmetry on the network parameters, while other projections are done using Eq. (6). To avoid the possible false minima convergence, in both methods we employ long Monte Carlo samples and try dozens of random initial approximations to select the best energy (see Appendix A).

E. Symmetry-breaking susceptibilities
One of the most direct characterizations of an ordered phase is through its symmetry breaking pattern. Consider an operatorÔ = Ω −1 iô i withô i acting locally, which measures symmetry breaking. In other words, the operator transforms non-trivially under the actions of the symmetry group. In a finite volume, data indicates that the pyrochlore ground state belongs to a trivial representation of all symmetries (point-group and spin), which forbids direct observation of the Ô = 0 condensate. Instead, one measures the (equal-time) susceptibility χÔ = Ô †Ô , which vanishes or survives in the thermodynamic limit if the phase has a symmetric or symmetry-broken ground state, respectively. Note that a non-vanishing susceptibility can only arise due to the establishment of long-range order.
To probe the SU (2)-symmetry breaking through longrange magnetic order, we calculate χM k , introducing the k-dependent operator where r i is the physical position of the site i and k takes values in the extended Brillouin zone. Since the Hamiltonian is SU (2)-symmetric, we restrict ourselves to the z-spin component only.
To probe the point-group symmetry-breaking tendency in the absence of magnetic order, we construct dimertype operatorŝ where with = 1, 2, 3.
Here, ξ ∈ {±1} and ω ∈ {1, exp (±2πi/3)} are the eigenvalues of inversion and rotation, respectively. Condensation ofÔ with non-trivial ξ or ω signals symmetry breaking. Note that the C 3 rotation does not mix the bond groups {(0, i), 1 i 3} and {(j, i), 1 i, j 3}, so the specific eigenvalue of C 3 does not fix the relative phase between these groups. The phases given in (10) are equal on opposite bonds within a tetrahedron, which is suggested by the dimer-dimer correlation pattern obtained within ED. Practically, such choice improves the signal-to-noise ratio in Monte-Carlo measurements of the susceptibility.

III. RESULTS
We present our results in the following order: First, we demonstrate that the variational energies we obtain compare favorably to the previous studies and then show finite size extrapolations based on computations with system sizes beyond those available in the literature. We show that our variational wave functions correctly capture the magnetic order in the k = 0 phase and the absence of magnetic order at the nearest-neighbor Heisenberg point. Second, we present our finite size extrapolation of the transition point between the two phases along the j 2 /j 1 axis for λ = 1. Third, we discuss the results of order-parameter and susceptibility calculations, elucidating the symmetry-breaking characteristics of both phases.

A. Accuracy of wave functions
The accuracy of the wave functions we obtained can be compared with recent DMRG data on clusters up to 4 × 3 3 [38] at the most frustrated λ = 1, j 2 /j 1 = 0 point. In Fig. 1 (b) we show the best energy values obtained within DMRG and within our work. The mVMC errorbars arise from statistical uncertainty and the errorbar at Ω −1 = 0 is estimated as 1/2 absolute difference between the value obtained on the largest 4 × 4 3 cluster and the extrapolation result [38]. The energies are listed in Table II of Appendix H. On all clusters for which DMRG data is available, our variational energies agree with or are lower than the ones obtained by DMRG 1 . Our result −0.477(3) for the infinite volume extrapolation is in agreement with the extrapolation from the DMRG calculations −0.488 (12) and improves upon the previous Variational Gutzwiller-projected mean field result [65]. This independent benchmark with DMRG confirms the ability of the variational Monte Carlo method to express the frustrated ground state even in the large volume.
Frustrated systems typically host many competing phases in a small energy window [40], thus a favorable comparison of energies is not always a guarantee of accuracy for a given variational ground state. However, we have further verified that physical observables obtained within the two variational parametrizations -mVMC and NQS -are also in striking agreement. To this end, we have studied the physical observables along the phase transition between the putative QSL phase at j 2 /j 1 = 0 and magnetically ordered k = 0 phase at large j 2 /j 1 (details about the k = 0 phase can be found in Appendix C).
The ground-state energy behavior along the phase transition qualitatively agrees across all our approaches. Figure 2 (a) shows the ground state energies at λ = 1 as a function of j 2 /j 1 obtained within ED and both variational methods on the equilateral clusters. The energy maximum signals the phase transition between the frustrated and the magnetic k = 0 phases. Note that a maximum is present for all approaches and clusters, but forms a more pronounced kink with larger cluster size, indicating a phase transition of first order.
The NQS CNN results look slightly off in the frustrated phase. However, at L = 2 the CNN energy falls between first and second excited states 2 , while overlap with the ground state is | ψ ED |ψ NQS |= 0.86, which suggests that observables computed with |ψ NQS will be dominated by the system ground state.
Further support that the wave functions are consistent and physically correct can be obtained by comparison of magnetic susceptibilities towards SU (2) spin-rotation symmetry breaking. In Fig. 3 (a-b) we show the magnetic susceptibility χM k plotted along the high-symmetry line in the 3D extended Brillouin zone. Both variational methods not only agree with ED and among themselves, but also clearly distinguish between the non-magnetic phase in Fig. 3 (a) and the magnetically ordered k = 0 phase in Fig. 3 (b) in accordance with the features described in Appendix C, such as peaks at k = X, L and ratio between their heights.

B. Magnetic phase boundary
Based on these consistent results, we are able to determine the shape of the phase boundary in the (j 2 /j 1 , λ) plane and extrapolate its location to the thermodynamic limit. As the position of the phase boundary, we take the maximum j 2 /j 1 such that the susceptibility on the whole XW line is higher than the peak's half maximum (see Fig. 3) 3 . Applying this criterion to the ED data for L = 2 yields the phase diagram shown in Fig. 1 (a). Within ED, the phase transition is located at j 2 /j 1 ∼ 0.06 on the SU (2)-symmetric axis. Further, in the interval 0.09 j 2 /j 1 0.12 there occurs a level crossing of excitations belonging to S = 1 and S = 0 spin sectors, indicating a transition to a magnetic phase [40,66,67]. With decreasing λ, the non-magnetic phase width decreases within ED and vanishes at the classical spin ice point, meaning that even an infinitesimal positive j 2 /j 1 term breaks the spin-ice degeneracy and establishes magnetic order. The ED phase boundary is well described by a (j 2 /j 1 ) * = 0.06λ 2 fit, shown as dashed line in Fig. 1 (a). The FWHM-analysis details, the raw ED phase diagram, and level spectroscopy can be found in Appendix D.
A complementary, direct signature of the phase transition comes from the non-monotonous behavior ("kink") of the energy as a function of j 2 /j 1 , shown in Fig. 2 (a). In the L = 2 ED-accessible case, the energy maximum is located at j 2 /j 1 ∼ 0.06, which coincides with the FWHM analysis. Larger clusters with L = 3 and L = 4 show similar behavior. To systematically obtain the kink position at large clusters, we perform a procedure comprising (1) obtaining the best wave function at metastable points belonging surely to the frustrated (j 2 /j 1 = 0.0) or the ordered (j 2 /j 1 = 0.2) phases, (2) adiabatically varying j 2 /j 1 , adjusting the wave function, and lastly, (3) locating the energy level crossing as an indication for the phase transition. An example of this hysteresis optimization is shown in the right panel of Fig. 2 (b). As the right panel shows, the energy level crossing is accompanied by an abrupt change in the ground state magnetic susceptibility. We find that the hysteresis optimization also significantly improves variational energies at the intermediate Using the hysteresis optimization, we perform the phase boundary infinite-volume extrapolation at λ = 1 and λ = 0.8. We consider the geometries 4×2 3 , 4×2 2 ×4, 64,96,128,196,256 and 108 sites, respectively, but perform the extrapolation using equilateral clusters of the shape 4 × L 3 only. The transition points for all methods and geometries are shown in Fig. 2 (c). Strikingly, the extrapolated range of the nonmagnetic phase at λ = 1 shrinks further as compared to the 4×2 3 system and the resulting critical value (j 2 /j 1 ) * = 0.0295 (30) is an order of magnitude smaller than the result obtained from FRG, (j 2 /j 1 ) * = 0.22(3) [22]. Note also that the phase-transition locations at the 4 × 3 3 cluster roughy agree within the two variational approaches, suggesting Spin/space-group spectroscopy at the SU (2)symmetric point. (a) Triplet gap as a function of inverse linear system size L. Data points labelled ED0.0, NQS 0.0 (difference between even and odd spin sectors projected using Eq. (7)), and mVMC0.0 represent the spin gap at j2/j1 = 0.0 obtained within ED, NQS, and mVMC, respectively, while mVMC0.2 shows the vanishing spin gap at j2/j1 = 0.2. Dashed lines on the right represent values obtained at non-equilateral clusters within mVMC at j2/j1 = 0.0. The errorbar at L −1 = 0 represents 1/4 of the difference between extrapolation and the value obtained on the largest possible cluster. The data can be found in Table II  that the conclusions are not subject to strong variational bias. We use the difference between the mVMC data infinite volume extrapolation and the result obtained on the largest 4 × 4 3 cluster as the estimate of extrapolation uncertainty, which is shown as an errorbar in Fig. 1 (a) and Fig. 2 (c) at Ω −1 = 0. Similarly, the λ = 0.8 point extrapolates to (j 2 /j 1 ) * = 0.0180 (35). The extrapolation results at λ = 0.8 and λ = 1 are shown with error bars in Fig. 1 (a). Taken together, our results strongly suggest that the phase boundary in the thermodynamic limit is located at substantially smaller j 2 /j 1 than the one extracted from L = 2 ED and FRG, but remains at finite j 2 /j 1 in the region λ 1. Since the initial ED phase boundary was well described by the parabolic fit, in Fig. 1 (a) we also include the phase boundary uncertainty region as bounded by the two parabolic curves in the thermodynamic limit as a guide to the eye 4 . We also point out that due to the shrinking width of the non-magnetic phase and its decreasing energy gain (as compared to the ordered k = 0 4 Note that our data does not allow us to reliably draw conclusion about the shape of the phase boundary in the region λ 1. The yellow region, shown in Fig. 1 (a) at λ 0.8 is a mere assumption. phase) 5 , variational approaches employed in this study converge to the ordered regime at λ 0.6 upon performing the hysteresis technique on the j 2 /j 1 = 0 axis. Thus, we infer the non-magnetic region shape at λ < 0.6 only from the ED data presented in Appendix D.

C. Symmetry breaking in the non-magnetic phase
We characterize the non-magnetic phase through its symmetry-breaking tendencies and contrast it with the behavior in the magnetic phase. We begin with the study of magnetic correlations and the spin gap, i.e., the energy difference between the lowest states in the S = 0 and S = 1 sectors. A vanishing spin gap allows for spontaneous spin-rotation symmetry breaking, characteristic of spin-nematic or magnetic phases. Numerous studies have already assessed the magnitude of the spin gap at λ = 1, j 2 /j 1 = 0 [12,14,23,34,38]. However, as seen in Table. II of Appendix H, no definite conclusion can be made, since the results are sensitive to the specific cluster geometry [34]. Here, we obtain the spin gap on all available clusters, but, most importantly, on three equilateral clusters 4 × L 3 , which retain the pyrochlore point-group symmetries under consideration. We thus believe that our data provide the most reliable spin gap extrapolation to date. As is seen from Fig. 4 (a), spin-gap extrapolations dramatically differ in the magnetic and the non-magnetic phase, represented by parameter choices j 2 /j 1 = 0.2 and j 2 /j 1 = 0, respectively. The gap extrapolates to zero and to a finite value ∆ S /j 1 = 0.40(4) in the magnetic and non-magnetic phases, respectively. The latter agrees with the recent DMRG result 0.36(3) [38].
A vanishing spin gap allows to establish an SU (2)breaking order parameter. We substantiate this via an extrapolation of the magnetic susceptibility χM X = M † XMX shown in Fig. 5 (a). Since the phase-transition point j 2 /j 1 is size-dependent, at intermediate j 2 /j 1 no such extrapolation is possible for the system sizes available to us and we show only j 2 /j 1 , which fall into one phase for all lattice volumes. For j 2 /j 1 0.025, χM X extrapolates to zero showing the absence of long-range magnetic correlations with this wave vector. Here we employ the Ω −1 extrapolation as we expect spin-spin correlations to decay exponentially with distance in the nonmagnetic phase [68]. Note that scaling Ω −r with r < 1 in the non-magnetic phase (e. g. due to a large correlation length) could only strengthen the conclusion of vanishing magnetic order. On the other hand, correlations at j 2 /j 1 0.12 extrapolate to finite values. To show that, we use the Ω −2/3 scaling arising from gapless Goldstone mode contributions to long-range spin-spin correlation in 3D [69]. The finite extrapolation values imply that the magnetic operatorM X will acquire a finite expectation value, M X = 0, in the thermodynamic limit in the ordered k = 0 phase, while the expectation value vanishes in the non-magnetic phase. Similarly to magnetic order and the corresponding SU (2)-symmetry breaking, we study the point-group symmetry breaking, which may be associated with dimertype order, for instance [70]. We compute the energy difference between the ground state and the lowest state in sectors with a different inversion and rotation eigenvalue. Only the equilateral clusters with 4 × L 3 obey rotation symmetry. In Fig. 4 (b) we show these gaps in the two phases at j 2 /j 1 = 0.0, 0.2 plotted against the inverse linear system size L −1 [40,66,67].
The inversion-symmetry gap extrapolates to a finite value in the k = 0 phase, which is in agreement with the fixed-point wave function for this phase, since the inversion operator does not mix the pyrochlore sublattices (see Appendix E for detailed description of magnetic phase and frustrated phase fixed-point wave functions). In contrast, the inversion-symmetry gap vanishes in the non-magnetic phase, opening the prospect of spontaneous inversion symmetry breaking.
The rotation-symmetry gap between the ground state (eigenvalue 1) and the lowest rotational doublet (eigenvalues e ±2πi/3 ) is found to be very small on the 4 × 2 3 cluster in both phases and falls below the error bar of our calculations in larger volumes, allowing for a spontaneous breaking of the rotation symmetry. In the k = 0 phase, this is in accordance with the sublattice order, which spontaneously breaks rotational symmetry. However, the rotation gap closing in the non-magnetic phase is a novel result.
To further support these results, we compute the sus-ceptibilities χÎ and χR to condensation of the operatorŝ corresponding to the inversion pseudoscalar and rotation E irreducible representation operators, whereÔ(ξ, ω) is defined in Eq. (9). In Fig. 5 (b-c) we show the corresponding susceptibilities for both phases, corroborating the gap analysis. In the non-magnetic phase the infinitevolume extrapolations of χÎ and χR are non-vanishing and, notably, χR is an order of magnitude larger than χÎ , which agrees with the relative gap magnitudes (see Appendix D). Likewise, in the k = 0 phase, only χR extrapolates to a nonzero value. For extrapolation, we employ the Ω −1 scaling based on our expectation that, if the dimerized phase indeed stabilizes, it breaks no continuous symmetry and has no gapless modes that would lead to non-exponential long-range dimer-dimer correlation saturation (see Fig. 14 in Ref. 71), as it happens with spin-spin correlations. If, on the contrary, no dimer order is truly established, dimer-dimer correlations decay algebraically and may cause scaling other than Ω −1 . Luckily, such change will lead to a scaling of the form Ω −r , with r < 1, which increases the slope against the one presented in the plots and thus further solidifies the conclusion. Real-space-resolved dimer-dimer correlations supporting this scenario are shown in Appendix G. We emphasize that if, despite this argument, the scaling dramatically differs from Ω −1 , e. g., L −1 , this would not allow us to draw a conclusion about rotation susceptibility extrapolation in the non-magnetic phase shown in Fig. 5 (c). The insets of Fig. 5 (b-c) show χÎ and χR obtained within ED. They show trends consistent with the larger-volume results.
A summary of operator expectation values based on these results is shown in the frames in Fig. 1 (a). Specifically, the rotation operator acquires a nonzero expectation value R = 0 in both phases, while the inversion operator is nonzero Î = 0 only in the non-magnetic phase.

D. General dimerization pattern analysis
Having observed non-vanishing susceptibilities towards establishment of rotation and inversion breaking dimerization patterns in the non-magnetic phase, we classify and measure (within ED/mVMC) all symmetry-allowed dimer observables, and compare them to those of analytically known model wave functions.
All the 12 × L x × L y × L z bonds of the pyrochlore lattice can be labeled with the unit cell index and an index within the unit cell. We will use the following convention: We choose each site of the 0-sublattice as the origin of a unit cell, where up-and down-tetrahedra meet, and assign the 12 bonds making up these two tetrahedra to the unit cell with the origin at site 0, as shown in Fig. 1 (a). Given a quantum state, we then compute the dimerization tensor χ D ij,µν = D µ iD ν j − D µ i D ν j with the indices 0 i, j < L x × L y × L z = Ω/4 running over unit cells in the lattice and 0 µ, ν < 12 enumerating bonds within unit cells. The eigenvectors that correspond to the largest eigenvalues of χ D , represent the dominant correlation patterns developed within the state of interest.
The ground state obtained within the 4 × 2 3 ED study belongs to the k = 0 momentum sector and the trivial irreducible representation of the point-group symmetry, while on larger symmetric clusters 4 × 3 3 and 4 × 4 3 this property is enforced by the procedure described in Section II D. Thus, χ D ij,µν shares the space group symmetries of the lattice (or finite cluster) and all its eigenstates transform according to an irreducible representation thereof. For all our numerically obtained variational wave functions, the Fourier transform of χ D ij,µν with respect to i, j indices has its dominant eigenvalues in the q = 0 sector. Hence, the dominant dimer-dimer correlation patterns are translationally-invariant, unlike for instance for the j 2 /j 1 Heisenberg model on the square lattice [40], where dimerization order breaks translational symmetry.
Restricting our consideration to translationallyinvariant (q = 0) eigenvectors of χ D µν (q = 0), we classify them by the irreducible representations of the O h point group acting on the remaining 12-dimensional linear space, namely the A 1g , A 1u , E g , E u , T 1g and T 1u representations. The details about those representations and the action of symmetry operations can be found in Appendix F. From the definition given in Eq. (12), we readily associate theÎ andR dimer operators, previously introduced in Eq. (12), with the A 1u and E g irreducible representations, respectively.
In Fig. 6 (a) we show the dominant eigenvalues of χ D µν (q = 0) obtained from ED (4 × 2 3 sites) and mVMC (4×3 3 and 4×4 3 sites) states. Notably, the E g irreducible representation is dominant. At the same time, we find as subleading eigenvector E u , which carries non-trivial quantum numbers of both rotation and inversion.The appearance of E g and E u as the irreducible representations of the largest eigenvectors reaffirms our conclusions from Sec. III C that both rotation and inversion symmetry are likely broken in this phase.
Next, we want to connect the observed χ D µν (q = 0) eigenvalue hierarchy to a specific dimerized state. To this end, we consider the only two short-range translationallyinvariant fully-dimerized patterns shown in Fig. 6 (c): two dimerized bonds in the same tetrahedron (denoted by |α ) or in opposite tetrahedra (|β ). In general, the state may be a superposition of the two |ψ(ϑ) = cos(ϑ)P |α + sin(ϑ)P |β , which we parametrize by the angle ϑ. To obtain a state that transforms trivially under all space group operations, as do our numerically obtained wave functions, we include the projectorP discussed in Eq. (6).
In the infinite volume, the eigenvalues of χ D µν (q = 0) are computed exactly for this family of model wave functions. We repeat this analysis on the small 4 × 2 3 cluster where one can store the wave function directly. In Fig. 6 (b) we show eigenvalues of χ D µν (q = 0) as a function of ϑ for the finite and infinite volume. The eigenvalues (scaled to match the simulation data) are shown in Fig. 6 (a) for tan(ϑ) = 1/2 (this choice leads to a good qualitative agreement between models and numerical results and was chosen as intermediate value between 0, where degeneracy between E g and E u is still present and 1, where crossing with low-lying eigenstates happens). We point out good agreement in the dominant eigenvalues hierarchy and suggest |ψ(arctan(1/2)) as a wave function qualitatively reproducing the numerical data. In this scenario, in the thermodynamic limit, when the point group symmetry is spontaneously broken, the wave function would be a superposition of the dimerization patterns |α and |β .
In addition, we discuss the pseudospin-ferromagnet wave function introduced, for instance, in Ref. [15]. Each isolated tetrahedron on the pyrochlore lattice has two orthogonal singlet states, which can be treated as states of the local pseudospin-1/2. We consider the pseudospinordered wave function introduced in Ref. [15], where the tetrahedra are grouped into 4 meta-sublattices as shown in Fig. 6 (d). We consider ferromagnetic pseudospin order in each of the sublattices. The pseudospin of the three sublattices (B, C, D) is obtained from the mean-field theory and results in a short-range dimerization pattern. The polarisation of A, expressed in terms of chiral superspin superpositions |± as (|+ + e iθ |− )/ √ 2 is obtained by accounting for quantum corrections and depends on the phase θ. In Fig. 6 (a)   projecting it withP to the trivial symmetry sector. Notably, the eigenstate structure is in sharp qualitative disagreement with our numerical data, for instance, due to appearance of q = 0 representation of χ D ij,µν among the leading contributions, which disfavours the pseudospin wave function scenario.

IV. DISCUSSION AND CONCLUSION
The spin-1/2 Heisenberg model on the pyrochlore lattice is an iconic candidate system for realizing a 3D QSL. This expectation is supported by the proximity of a U (1) QSL in the vicinity of the Ising point and the absence of magnetic correlations at the Heisenberg point [22][23][24]65]. To date, no consensus on the existence of the QSL phase has been reached, because the 3D many-body spin problem challenges all numerical techniques available. Confirming the QSL nature of a ground state is intrinsically harder than establishing an ordered phase: besides proving the absence of long-range order, emergent fractionalized excitations should be found [72].
We have performed a comprehensive numerical search for spontaneous symmetry breaking in the putative QSL phase and quantified its extent in phase space. Combining two complementary variational approaches and exact diagonalization, we are able to gather strong evidence for a symmetry broken rather than a featureless QSL nature of the phase, based on extrapolations to the thermodynamic limit. We characterize the phase as not magnetically ordered, but spontaneously breaking rotation and inversion symmetry by establishing long-range dimer order. This contradicts the QSL realization at the Heisenberg spin-1/2 pyrochlore as discussed in Ref. [72]. Rather, the symmetry breaking pattern for the nonmagnetic phase is consistent with the fixed-point dimerized wave function shown in Fig. 6 (c). Namely, simultaneous breaking of rotation and inversion symmetries by the dimer order would also be observed in the fully-dimerized state. In support of this, we found direct numerical evidence for dimer correlations in the ground state.
Its strong geometrical frustration makes pyrochlore materials promising candidates for experimental observation of non-magnetic spin systems. For instance, spin-1 pyrochlore NaCaNi 2 F 7 with j 2 /j 1 ∼ −7 × 10 −3 was recently claimed to show spin-liquid-like behavior in neutron scattering experiments down to low temperatures [73]. However, a spin-1/2 compound with nearly SU (2)-symmetric and dominantly nearest neighbor AFM interactions is still to be identified.
We find that already a moderate admixture of next-tonearest-neighbor interactions establishes long-range magnetic order, emphasizing that only materials in a narrow parameter regime are expected to show non-magnetic ground states. Concretely, we showed that the width of this non-magnetic phase is (j 2 /j 1 ) * = 0.0295 (30), which is substantially smaller than predicted in a recent FRG study [22].
The seemingly small extend of parameter range is comparable to what is seen in other models with candidate QSL phases as well, in particular since it also extends to negative j 2 . For instance, for the j 2 /j 1 square Heisenberg model, the width of a putative QSL phase was found to be ∼ 0.05 [40].
Due to the shrinking of the non-magnetic phase as λ is tuned from 1 to 0, the energy competition with the magnetic phase becomes even tighter and the local minima problem prevents us from obtaining a reliable non-magnetic variational wave function at small λ and j 2 /j 1 = 0 for the U (1) π QSL. Severe finite size effects existing in this regime should also be noted: the 4 × 2 3 cluster hosts no U (1) π QSL state within perturbation theory on the spin-ice manifold as the λ 3 -order hexagon flips are dominated by λ 2 terms specific to this cluster size [7]. As a result, we were unable to locate the transition point between the observed symmetry-broken phase at λ = 1 and the well-known U (1) π QSL phase at λ = 0. This transition, if it indeed exists, represents a crucial element of the pyrochlore puzzle that is yet to be resolved within future studies.
In summary, using many-variable Monte Carlo methods, we were able to refine the phase diagram of the spin-1/2 pyrochlore Heisenberg model both qualitatively and quantitatively, assembling strong evidence against a featureless QSL ground state. This still leaves open the possibility of a symmetry broken state with concomitant topological order. The NQS method casts each the spin configuration σ in S z basis to bit representation σ i ∈ {−1, +1} and uses it as an argument of a Neural Network Ψ. The result Ψ(σ) is interpreted as the (not normalized) wave function component corresponding to the basis vector σ. In case of a feed-forward neural network, which is most commonly applied in the NQS method, the input v 0 = σ undergoes a sequence of transformations whereŴ i is the weight matrix and b i is the bias vector. Parameters (Ŵ i , b i ) are the variational parameters of the NQS ansatz. The linear transformationŴ i v i + b i is followed by the application of nonlinearity f , which is an essential ingredient of the procedure, since otherwise only linear functions could be encoded. In this work, we apply the ReLU (rectified linear unit) nonlinearity [75] which is commonly used in artificial intelligence applications. In NQS, fixing quantum numbers of the point symmetry group is done similarly to mVMC, using Eq. (6). Imposing translational symmetry on the level of variational parameters requires construction of manifestly translation-invariant neural network. This "hard-coded" translational invariance turns out to be crucial for obtaining significant overlap with the QSL phase ground state [52], as it effectively increases the number of Monte Carlo samples L x × L y × L z times.
Translational invariance is achieved by using convolutional layers with periodic padding. Consider a spin configuration v 0 l,x,y,z = σ l,x,y,z ≡ σ l,r to be rearranged in a 4D-tensor of shape (C 0 = 4, L x , L y , L z ). Then, application of one convolutional layer reads Here, the kernel matrix K ll (∆r) depends only on the relative distance ∆r between elements v i+1 l,r and v i l ,r+∆r , while bias b i l depends only on the sublattice index. If r + ∆r extends beyond the boundary, periodic boundary conditions are applied, which makes this architecture manifestly translation invariant. Note also that the output vector v i+1 l,r has shape (C i+1 , L x , L y , L z ), so the number of channels C i+1 can vary while the spatial dimensions remain untouched.
As output the network should produce wave function in the form Φ(σ) = (log A, e iφ ), where log A is the amplitude logarithm and e iφ is the complex phase of the wave function element. In this work, we train two separate neural networks, one of them producing amplitude, and the other one producing phase, which was shown to increase precision in a number of cases [76,77]. So, the network should output one number obeying translational l,r is first rearranged into the (4, Lx, Ly, Lz) tensor, and then undergoes a sequence of convolutions in accordance with Eq. (A2). After the last convolution, the tensor (C, Lx, Ly, Lz) is spatially averaged (in this example C = 16) and is fed to the fully connected layer, given by Eq. (A1). The resulting number is then interpreted either as log|Ψ(σ)| or argΨ(σ).
invariance. To construct one number from the hidden tensor of shape (C, L x , L y , L z ), we first take the mean value over all spatial dimensions (C, L x , L y , L z ) → C, and then apply a standard dense layer that maps Cdimensional vector to phase or amplitude (Eq. (A1)).
Recently it has been shown that on moderately large systems, absence of "hard-coded" translational invariance might lead to outstanding wave function approximations [55]. Instead of CNN, one uses a standard shallow fully connected network with just one layer inspired by Restricted Boltzmann Machine (RBM) architecture used in early NQS studies [41]. In this way, no symmetries are encoded in the NN parameters, but rather all symmetry projectors, including momentum, are applied at the end. In this work, we also employ this architecture to improve the CNN results on the 4 × 2 3 and 4 × 3 3 systems.
Obtaining significant overlap with the QSL phase ground state is known to be a nontrivial task within the NQS method. The problem can be effectively reformulated in terms of the wave function sign structure (note that since the Hamiltonian Eq. (2) is real, phase of any element can be chosen φ = 0, π). In an ordered phase, wave function element's sign can be easily inferred from the spin configuration σ, e.g. via the Marshall-Peierls sign rule [78]. In this case, the network usually shows extreme accuracy, even outperforming existing approaches [41,42]. However, as the system is moved away from the ordered phase towards the QSL phase (for instance, by increasing j 2 /j 1 in the spin-1/2 Heisenberg model on square lattice), the overlap may drop to almost zero [42,52,77]. It was shown that in the QSL phase the wave function phase structure has a much larger complexity [77], which might be difficult to catch if a wrong training method or network architecture is applied. In case of 2D frustrated magnets it was shown that the right choice of network architecture can increase overlap with the ground state from 0 to 0.9 at the maximally frustrated point [52].
In case of 3D frustrated magnets, the architecture alone turns out to be not enough to grasp the correct QSL ground state properties. To deal with frustrations, we introduce a novel algorithm of alternating learning.
It was initially shown in [77] that one can improve the final training result by performing training in two stages. During the first stage, one sets amplitudes of all spin configurations equal log|Φ(σ)|= 0 and trains only phases. During the second stage, one trains both phase and amplitude networks simultaneously. We extend this idea and add the intermediate alternating phase, consisting of many tiny alternating trainings of phases with amplitudes fixed and amplitudes with phases fixed. Empirically, this extension can be motivated as follows. Note that within the NQS method the phase lives on the complex circle, while the true ground state phases can be defined purely real ±1. As the result, a converged network phase distribution usually has two major clusters at φ and φ + π, where φ is an arbitrarily found physically irrelevant wave function gauge phase. After full convergence, the two phases are separated by a "potential barrier" and phase of a single spin configuration can not "tunnel" between these minima. However, if the amplitude landscape is totally flat, log|Φ(σ)|= 0, then the tunneling can happen with a higher probability. During the first training stage, we train only phases to find the best energy available with the constraint log|Φ(σ)|= 0. After this constraint is removed, amplitudes quickly converge and the potential barrier between φ and φ + π sets in. If, however, the amplitudes are trained only for a small time after which the phases are trained to adjust for slightly changed amplitudes, we allow some signs to "tunnel" between the two minima since the barrier is not yet too high. the vicinity of the phase transition, however, the two phases have a tight energy competition. Variational methods, being usually prone to ordered solutions, might get trapped in false minima and require a large number of random initial approximations to resolve the correct phase near the critical point. It turns out, however, that the search for the best variational parameters can be efficiently replaced with so-called "hysteresis optimization".
The first step is to optimize the wave function at two values of j 2 /j 1 lying deep within the adjacent phases. The two sets of optimized wave function variational parameters are then used as the initial approximation for the hysteresis procedure. Namely, starting with the wave function trained deep within the k = 0 phase, we gradually decrease j 2 /j 1 , optimize the wave function for the new value of j 2 /j 1 , and obtain the energy E k=0 (j 2 /j 1 ). During the evolution, we observe the susceptibility pattern χM k to verify that the wave function still has the initial phase features. Repeating the same procedure in the opposite direction by starting with the QSL phase wave function, we obtain E QSL (j 2 /j 1 ) that intersects E k=0 (j 2 /j 1 ) at some j 2 /j 1 = (j 2 /j 1 ) * , which provides us with an accurate estimate for the phase transition position. An example of the hysteresis optimization with mVMC for the largest 4×4 3 lattice is shown in Fig. 2 (b). As the right panel shows, adiabatically evolved energies cross at j 2 /j 1 ∼ 0.032. At this point, wave function of the magnetic phase starts to outperform the one of nonmagnetic phase in terms of energy. Importantly, at any j 2 /j 1 the best energy obtained with the hysteresis optimization was never greater than the ones obtained with numerous trials starting from random approximations. Thus, the procedure provides us with a better energy in addition to saving the computational time. At every point we select the best available wave function and plot the magnetic order, which abruptly establishes if j 2 /j 1 is tuned above the phase transition point. This is consistent with the notion of metastable wave functions defined within magnetic and non-magnetic phases.
Appendix C: Properties of the magnetic k = 0 phase Momentum resolved spin-spin correlation functions are the most pronounced ordered phase fingerprints and thus are a powerful tool to track phase transition between ordered and disordered phases [22,39,40,68]. In terms of χM k = M † (k)M (k) , the frustrated phase and the neighboring k = 0 ordered phase have several distinctive features. The phase, known also as the "sublattice ferromagnet", has spins within each sublattice ordered ferromagnetically, while still maintaining the zero-sum spinice rule within each tetrahedron. In case of perfect classical ordering, there are C 2 4 = 6 possible ground states satisfying this requirement. Each of these states can be tracked by two peaks of χM k=X at a pair of inversionrelated X-points (see Fig. 3). These states also have equally high peaks at all L-points (middle points of the large BZ faces) with the corresponding susceptibilities ratio χM k=X /χM k=L = 4. In Fig. 8 (a-b) we show correlations magnitude over the 3D extended Brillouin zone and the k x = k y cut called [hhl].
The frustrated phase, as compared to the k = 0 ordered phase, lacks any long-range magnetic order, χM k = 0, ∀ k ∈ BZ. The susceptibility multiplied by the system volume ΩχM k , however, is finite and governed by the local spin-spin correlations. These correlations show the "diffusive" behavior, i.e. distributed over the BZ boundary (see Fig. 8 (c-d)). On the [hhl] cut, the bow-tie pattern in the vicinity of k z = 4π (vertical axis) is the attribute of the local spin-ice rule [22].
Appendix D: ED data: spectroscopy and phase diagram A qualitatively correct result for the shape of the phase diagram in terms of frustrated and magnetic phases can be obtained with ED on a 4 × 2 3 cluster. Beforehand, note that the k = 0 phase, though having qualitatively similar features on both λ = 0 and λ = 1 axes, such as magnetic susceptibility peaks at X and L points, has different peak heights. In the classical λ = 0 case, even a tiny positive j 2 /j 1 perfectly orders the system, giving rise to extreme χM k=X = 1/12. Contrary, in the quantum λ = 1 case, quantum corrections and SU (2) symmetry do not allow perfect classical sublattice ferromagnetic ordering, which significantly reduces the peak. This is illustrated in Fig. 9 (c), where we consider fixed-j 2 /j 1 cuts of the phase diagram. Note that indeed if j 2 /j 1 > 0, the peak height χM k=X would reach 1/12, provided λ was sufficiently small. Note also that as λ > 1, magnetic correlations are further reduced.
To unite these behaviors in the vicinity of λ = 0 and λ = 1 points, we consider the full width at half maximum (FWHM) characteristic defined as where FWHM is a vector collinear to X−W. In Fig. 9 (d) we show the peak width |FWHM|/2π in the (λ, j 2 /j 1 ) plane. If there is no such momentum k along the XW line to satisfy the FWHM criterion, we put |FWHM|/2π = ∞.
Exact diagonalization provides another way to pinpoint phase transition between magnetically-ordered and frustrated phase by locating positions of level crossings between different total spin S excitations [40]. In Fig. 9 (a) we show the energy gaps of excitations with distinct momenta and spin. Notably, the k = 0 excitation with S = 1 gradually softens and in the region 0.09 j 2 /j 1 0.125 it crosses non-zero momentum excitations of S = 0. This region may be taken for the spectroscopy estimation of the phase transition point, which is in rough agreement with the FWHM analysis. Note that softening of S = 0 mode is a signature of tendency to spontaneous SU (2)-symmetry breaking and stabilization of the magnetic phase.
Finally, in Fig. 9 (b) we show inversion and rotation gaps as a function of j 2 /j 1 at λ = 1 within ED, that complements the susceptibility analysis in Section III C. Note that the inversion gap is significantly larger than the rotation gap, which explains why the susceptibility to rotation breaking was found to be of order of magnitude larger than the inversion breaking susceptibility. Also, the inversion gap decreases towards the non-magnetic phase, which hints developing tendency to inversion symmetry breaking, observed in Section III C within susceptibility analysis. Appendix E: Spontaneous symmetry breaking from ED data In Fig. 9 (b) we show that even on the ED-accessible 4 × 2 3 cluster the rotation gap between the ground state s-orbital and first excited p x /p y -orbitals is vanishingly small. As shown in Section III C this gap vanishes in the thermodynamic limit.
To get more insight we apply the procedure suggested in [79]. Assuming that these three low-energy states are "already degenerate", we can linearly superpose them at no energy cost. To probe the symmetry breaking, we select a bond dimer operatorD ij =Ŝ i ·Ŝ j , residing on the fixed bond AB (see Fig. 10 (c)), and study how its expectation may change as the ground states are superposed.
Within ED, we measure the 3×3 pair-wise expectation value matrixÔ where 0 α, β 2 index the approximately-degenerate ground states and the two excited states. Eigenvalues λ Dij of O αβ stand for theD ij expectation values on superpositions whereD ij is diagonal. We remind here that the doublet gap only becomes smaller as j 2 /j 1 is increased to positive values, as seen in Fig. 9 (b). In Fig. 10 (a) we show λ Dij eigenvalues as a function of j 2 /j 1 . Note that a nearly-degenerate pair of λ Dij splits in the limit of large j 2 /j 1 and the intermediate eigenvalue switches to become near-degenerate with the upper eigenvalue. This dramatic rearrangement of λ Dij is a yet another signature of a phase transition. To get more insight, we describe this behavior in terms of fixed-point classical wave functions. At large j 2 /j 1 , when the system is in the k = 0 phase, the classical ground state has each sublattice ordered ferromagnetically, maintaining the zero spin per tetrahedron. In total, there are 4! /(2! 2! ) = 6 ways to construct such ground state. An example of such ground state is shown in Fig. 10 (b). Since our case is not classical and the SU (2)-symmetry is present, we enforce the symmetry s z → −s z by turning 6 such ground states into 3 equal-weight superpositions symmetric under s z → −s z spin flip. Each state of this basis breaks rotation symmetry. TheD ij bond operator is diagonal in this basis and its expectation equals ±1/4 if the adjacent spins are (non-)collinear. Thus, the eigenvalues are 4λ Dij = {+1/4, −1/4, −1/4}, which is in a striking agreement with the right side (magnetic part) of Fig. 10 (a).
In the other limit, j 2 /j 1 → 0, where magnetic ordering is absent, the levels picture can be well described by the dimerized pattern breaking rotation symmetry. An example of such pattern is shown in Fig. 10 (c). In a conventional dimerization pattern, each spin belongs only to one dimer. In our case, since the inversion symmetry is preserved (we only superpose inversion-symmetric states), a ground state must be taken as an equal-weight superposition of the red pattern and blue patterns, dimerizing up-tetrahedra or down-tetrahedra respectively and related by the inversion symmetry. This pattern can be rotated, producing three distinct ground states, each of them breaking rotational symmetry. The dimer operator is diagonal in this basis as well, and has the eigenvalues read λ Dij = {−3/8, 0, 0}, which is again in qualitative and even quantitative agreement with the left side (dimer part) of Fig. 10 (a).
Note, that we by no means claim that these tentative magnetized and dimer states are the ground states of the true system. For instance, the dimerized pattern has a much lower energy and only a 30 % overlap with the ground state. However, we note that the observed eigenvalues λ Dij evolution picture is in acute agreement with this fixed point rotation-breaking wave functions treatment. We thus conclude that, since the rotation doublet gap closes in thermodynamic limit, the ground state may easily form a superposition with a strongly broken rotation symmetry, where the word "strongly" means that the magnitude of the symmetry-breaking operatorD ij is in striking agreement with the very typical magnetic and dimer symmetry-breaking states. Here ω = exp (2πi/3). The inversion eigenvalue λ = +1 corresponds to g-representations and λ = −1 corresponds to u-representations. The E g/u ± basis states span the 2-dimensional representation and T g/u a/b/c basis states span the 3-dimensional representation. we expect that all excitations are gapped. As the result, dimer-dimer correlations saturate with distance exponentially fast D iDj ∼ A + B exp(−r/ξ). This fast saturation leads to the inverse-volume finite-size correction scaling 6 . Following Ref. [71], in Fig. 11 we show D 01 0D α j -D 01 0 D α j , where α labels the 12 bonds that belong to a unit cell as introduced in Sec. III D and in Fig. 10, as the function of distance |r 0 − r j | between the 6 Strictly speaking, following this logic, the Ω −1 extrapolation used for inversion-breaking susceptibility in the magnetic phase at j 2 /j 1 = 0.2 in Fig. 5 should not be used in this case, as we assume no dimer order establishment. In absence of dimer order, dimer-dimer correlations will saturate polynomially and will lead to L −α scaling with α < 3. This, however, will not change the conclusion since a steeper slope will only lead to an even smaller extrapolation result.
unit cells 0 and j, measured in the non-magnetic phase at the 4 × 4 3 cluster within mVMC. Notably, the correlations quickly saturate to non-vanishing values, which are dependent on α. In another words, we expect small correlation length ξ L, which paves way to finite dimer order parameter susceptibility and justifies the scaling employed in Fig. 5.

Appendix H: Ground state energies and spin gap
The energy-per-spin is the foremost way to assess the quality of the variational wave function. The spin-1/2 pyrochlore ground state energy at j 2 /j 1 = 0 was previously estimated within various approaches [13,21,24,65,81,82]. In Table II we summarize the ground state energies from this work as well as other references which employ DMRG or exact diagonalization. Note that the quasi-planar geometry used in exact diagonalization of 28-and 36-spin clusters in [34] has a significant effect on the ground state energy, which emphasizes strong geometric dependence of the results and the importance of choosing clusters with the maximum possible number of spatial symmetries.
Similarly, the non-vanishing spin gap is important not only for proving the absence of magnetic order, but is also essential for the perturbative dimer model analysis of the pyrochlore [14,15]. In Table II  clusters. Note again the strong dependence of spin gap on the cluster geometry.