Possible inversion symmetry breaking in the $S=1/2$ pyrochlore Heisenberg magnet

We address the ground state properties of the long-standing and much-studied three dimensional quantum spin liquid candidate, the $S=\frac 1 2$ pyrochlore Heisenberg antiferromagnet. By using $SU(2)$ DMRG, we are able to access cluster sizes of up to 108 spins. Our most striking finding is a robust spontaneous inversion symmetry breaking, reflected in an energy density difference between the two sublattices of tetrahedra, familiar as a starting point of earlier perturbative treatments. We also determine the ground state energy, $E_0/N_\text{sites} = -0.488(12) J$, by combining extrapolations of DMRG with those of a numerical linked cluster expansion. These findings suggest a scenario in which a finite-temperature spin liquid regime gives way to a symmetry-broken state at low temperatures.

Introduction.-Frustrated magnets, on account of exhibiting many competing low energy states, are a fertile ground for exotic physics. A celebrated example is the pyrochlore Heisenberg antiferromagnet, which resides on a lattice of corner sharing tetrahedra, depicted in the inset of Fig. 1. The classical Heisenberg model on this lattice has a highly degenerate ground state [1], forming a classical spin liquid [2] with an emergent gauge field [3].
In contrast, the ground state of the quantum pyrochlore antiferromagnet remains enigmatic.
While recent experimental evidence in the approximately isotropic S = 1 compound NaCaNi 2 F 7 shows a liquid like state down to low temperature [4], the S = 1/2 case is still open both in theory and experiment.
Theory work on this promiment quantum spin liquid candidate over the years has been formidable. Absent a systematically controlled method, various approaches have somewhat inevitably led to an array of possible scenarios. One strand of work has built on a perturbative approach, in which half the couplings (those on one tetrahedral sublattice) are switched on perturbatively. This has led to suggestions of a ground state which breaks translational and rotational symmetries [5][6][7], a valence bond crystal [8] or a spin liquid state [9]. On top of this, the contractor renormalization method [10] finds antiferromagnetic ordering in a space of supertetrahedral pseudospins, pointing to an even larger real-space unit cell. To render the problem more tractable, all these theories involve the derivation of an effective Hamiltonian, which is per se not exactly solvable and hence solved by some type of approximation, ranging from mean field theory to classical Monte Carlo numerics. On a different axis in theory space, parton-based theories yield an ordered state with a chiral order parameter [15] or a monopole flux state [17], while pseudofermion functional renormalization group suggests a spin liquid ground state [18] In view of this relatively wide range of ground state candidates, a controlled and unbiased treatment of the model is clearly desirable, if only to narrow the possible location of the goalposts somewhat. Unfortunately, most numerical approaches quickly reach their limits for frustrated magnets in d = 3. While exact diagonalization is currently limited to ∼ 48 sites [19], possible alternatives are series expansions such as the numerical linked cluster expansion (NLCE) [20][21][22][23][24][25][26][27][28][29][30][31][32][33] or high temperature expansions [34,35], which can be pushed down to low temperatures [33], although they do not provide access to the groundstate itself and are particularly challenged by many competing low energy states.
Here, we take DMRG one step further, by applying it to the pyrochlore lattice in d = 3, and present a study of periodic clusters with N sites = 32, 48, 64, 108. This demonstrates that DMRG can treat clusters with up to 108 sites reliably, significantly larger than previous exact diagonalization results of 36 sites [46]. Exploiting the SU(2)-symmetry of the model [47][48][49][50], we keep up to 20000 SU(2) states, (typically equivalent to 80000 U(1) states). We calculate the ground-state energy, the spin structure factor and low-energy excitations for these clusters, yielding an estimate for the ground-state energy per site in the thermodynamic limit of E 0 /N sites = −0.488 (12). The study of finite size clusters is complemented by a high order NLCE calculation, which excludes any scenario where E 0 /N sites > −0.471.
Our main finding is that the ground state of the larger (64 and 108-site) clusters we consider exhibits a breathing instability, rendering up and down tetrahedra (cf. inset of Fig. 1) inequivalent: one tetrahedral sublattice exhibits a lower energy than the other. Amusingly, our estimate for the ground state energy is compatible with that of the original perturbation theory with a simple mean field solution of the resulting effective Hamiltonian, where the inversion symmetry was maximally broken at the very outset of the calculation [5].
Model and methods.-We consider the pyrochlore antiferromagnetic Heisenberg model with S = 1/2: where the spins sit on the sites i, j of the 3d pyrochlore lattice and i, j denotes nearest neighbors. The lattice is a face centered cubic lattice with lattice vectors 2 a 3 , such that each lattice point can be expressed by R α,n1,n2,n3 = n 1 a 1 + n 2 a 2 + n 3 a 3 + b α , with integer n 1 , n 2 , n 3 and α ∈ {0, 1, 2, 3}. The model is obviously SU(2) symmetric. Our DMRG calculations are performed on finite size (N = 32, 48, 64, 108) clusters with periodic boundary conditions (cf. Tab. II of Appendix).
We apply the one-and two-site variants of SU(2) DMRG to reach high bond dimensions necessary to obtain reliable results in our three dimensional clusters. Since DMRG requires a one-dimensional topology, we  (2) bond dimensions χ (indicated by the labels) as a function of the two-site variance. Solid lines correspond to linear extrapolations to the error free limit, corresponding to infinite bond dimension and zero variance. We estimate the systematic extrapolation error as the half distance between the last point and the extrapolated value.
impose a one dimensional "snake" path on the threedimensional lattice, which defines the variational manifold. We use fully periodic clusters to reduce boundary effects and confirm that using a snake path which minimizes the bandwidth of the connectivity matrix improves convergence [33,51,52].
For small bond dimensions (χ 2000) we use the twosite version of the DMRG, and switch to the one-site variant to optimize the wave function for larger χ. Since the truncation error is not well-defined in the one-site variant case, we use the reliable two-site variance estimation to extrapolate towards the error-free case [53], because calculation of the full variance would be impractical due to its cost.
It turns out that even the calculation of the two-site variance becomes too costly for clusters with more than ∼ 100 sites, in which case we revert to the standard extrapolation as a function of the inverse bond dimension (cf. also [54]). We checked for several cases where both types of extrapolations were possible that they give the same results within the margin of error.
Ground-state energy.-Using DMRG, we calculate the variational ground state energy of finite clusters with high accuracy. By systematically increasing the bond dimension χ, we enlarge the variational manifold in a controlled way, such that we can extrapolate, χ → ∞, to the exact limit using a linear extrapolation as a function of the two site variance (cf. Fig. 2). We use an estimate of the systematic extrapolation error given by half the distance between the extrapolated value and the last DMRG point. Fig. 1 shows the extrapolated energies per lattice site of all finite clusters we considered in comparison with the available predicted ground-state energies in the literature. Our results show a monotonic growth of the ground-state energy as the number of sites is increased.
The periodic clusters we consider have either the full cubic (32,108) or an increased/reduced (48a, 48b, 48c, 48d, 64) symmetry of the pyrochlore lattice and represent the bulk due to the absence of a surface. The energies per site of different clusters as a function of inverse cluster size admit a fit to a quadratic polynomial, which we use to obtain an extrapolation to the thermodynamic limit. In order to get an estimate of the extrapolation error, we use gaussian resampling, using the systematic DMRG errorbars as standard deviation. This yields our best estimate for the ground state energy of E 0 /N sites = −0.488 (12). In this fit we considered only the cluster 48d among the 48-site clusters, which appears to be consistent with the other clusters, while other 48 site clusters have lower ground state energies.
Our extrapolated (χ → ∞) cluster energies and gaps are summarized in Table I. While the singlet gaps in the most symmetric clusters (32, 48d) are very small, the triplet gaps are sizable and roughly an order of magnitude larger. Since the 48d cluster does not obey all lattice symmetries, a reliable extrapolation is not possible, but our results are compatible with a scenario with a finite triplet gap, in which case all low energy excitations would be in the singlet sector as claimed in Refs. [6,10].
Our finite temperature NLCE [20][21][22][23][24][25][26][27][28][29][30][31][32][33] provides a complementary perspective. We have carried out this expansion in entire tetrahedra up to eight order (cf. [33] for details, as well as Appendix C), obtaining convergence for the energy per site in the thermodynamic limit as a function of temperature for temperatures T 0.2. Since the energy is a monotonic function of temperature, the converged part of E(T ) (cf. Fig. 7) provides an upper bound for the ground state energy E nlce ≈ −0.471J, which is consistent with the DMRG data and extrapolation. One can furthermore polynomially extrapolate the finite temperature NLCE energies to zero temperature (assuming an analytic behavior at low temperatures), see Fig. 7, and obtains −0.495 (15), which agrees remarkably well with the DMRG extrapolation and lies within its error bar, serving as a further corroboration of the DMRG energy. In light of these results we can confidently exclude a ground state energy per site larger than −0.47J.
Ground state symmetry-breaking.-To investigate the properties of the ground state in more detail, we calculate the total spin, and hence total energy, of up and down tetrahedra separately. This reveals an inequivalence of up and down tetrahedra (cf. Fig. 6 in the Appendix), suggesting a breaking of the inversion symmetry of the lattice. In our DMRG calculations, the snake path does not fully respect the symmetry between up and down tetrahedra, so that we need to verify that this symmetry breaking is intrinsic, and not due to a preference imposed by the snake path. We therefore introduce a small symmetry breaking 'breathing' perturbation, where we modify the couplings of up/down tetrahedra to be J = 1 ± ,  equivalent to the standard technique of including pinning fields. Fig. 3 shows the results for the total spin of up and down tetrahedra for opposite signs of the breathing perturbation in the 64 (108) site clusters as a function of the two-site variance (inverse bond dimension), admitting a linear extrapolation towards χ → ∞. The results reveal a clear selection of states with opposite symmetry breaking, as required for spontaneous symmetrybreaking. The order parameters for the larger, 108 site, cluster are slightly different for the two opposite pinning fields (Fig. 3, right panel), but that difference is much smaller than the extrapolated order parameter which differs only little between the two clusters. It is of course always possible in principle that the symmetry breaking vanishes when yet larger clusters are considered. Given the scaling of the computational effort with system size, the study of much larger clusters with the present method is, however, out of reach.
We next consider nearest neighbor spin correlations of the best (lowest-energy) wave functions |ψ 0 obtained in DMRG. For each pair of adjacent sites (i, j), we calculate the correlation function C ij = ψ 0 | S i · S j |ψ 0 . We plot the result for the clusters 64 and 108 in Fig. 4 (truncated  FIG. 4. Real space spin correlation Cij in the ground state (Sz = 0) for N = 64 (left) and N = 108 (right) shown in the cubic unit cell. The thickness of the red bonds corresponds to magnitude of the correlation between neighboring sites. The black lines are indicating bonds between sites with negligible correlations. to the cubic unit cell for ease of visualization), with the tube thickness proportional to the strength of the spin correlations.
The correlation pattern reveals that one sublattice (say, 'up') of tetrahedra contains more strongly correlated bonds than the other. These are found on opposite edges of 'up' tetrahedra. We note that the details of this pattern still depend strongly on the cluster geometry and we get opposite choices of correlated bonds in the two clusters, presumably due to different symmetry broken states picked by the different 'snake' paths in the two clusters.
Ground-state structure factor.-The static spin structure factor for different clusters, accessible in neutron scattering experiments, is obtained from the Fourier transform of the spin correlations (factor 4/3 from nor-malization 1/(S(S + 1)) for spin S = 1/2): where R i denote the real-space coordinates of sites and the index c denotes the connected part of the correlation matrix. The results for two cuts (Q x = Q y (top) and Qz = 0 (bottom)) in the three dimensional momentum space are shown in Fig. 5. One can readily recognise the bow-tie patterns, the hallmark of pyrochlore magnets [3-5, 9, 16, 18, 33]. Note that the 32-and 108-site clusters have full cubic symmetry, while the 64-site cluster does not, hence the structure factors looks slightly different in that case. The results for the spin structure factor and the absence of sharp Bragg peaks confirm that there is no long range magnetic ordering. The observed pattern is very close to what is found at finite temperature in the regime T 1 [33]. While the pinch points sharpen with increasing system size (and therefore momentum resolution), we are unable to extrapolate their width reliably to the thermodynamic limit to extract a correlation length. Interestingly, the breaking of lattice inversion symmetry is not straightforwardly visible in the magnetic structure factor.
Concluding discussion.-Our DMRG study has found the ground state of the SU (2) symmetric S = 1 2 Heisenberg antiferromagnet to discard lattice inversion symmetry in favour of a 'breathing' pattern of strong (weak) sublattices of up (down) tetrahedra. We extrapolate the energy per lattice site to −0.488 (12). The possibility of such spontaneous symmetry breaking has been a central question for this class of magnets, as several studies have used an explicit such symmetry breaking as a starting point of various perturbative schemes [5,6,10,55]. As the restoration of an explicitly broken symmetry in a perturbative scheme is generically not to be expected, a non-vanishing order parameter does not per se indicate spontaneous symmetry breaking.
Our results are thus important in that they provide largely unbiased evidence for the existence of this spontaneous symmetry-breaking, subject only to finite-size effects which are much reduced in comparison to previous studies. This also indicates that one of the prime Heisenberg quantum spin liquid candidates in three dimensions in fact exhibits at least one form of symmetry breaking.
In closing, we note that our extrapolated ground state energy lies close to the estimate obtained in the pioneering work by Harris et al. [5], in the abovementioned scheme of coupling the up tetrahedra perturbatively through the bonds of the down tetrahedra. These authors also found a long-range dimer ordering (cf. also [6]) compatible with the correlation pattern we observe in our calculations shown in Fig. 4. This first, simple and quite uncontrolled, approach to this difficult problem thus may turn out to have been already quite close to what will eventually be established as the final answer. The real space dimer correlation pattern shown in Fig.  4 suggests that the lattice inversion symmetry is broken in the ground state. In order to scrutinize this finding, we analyze the square of the total spin (morally the tetrahedron energy) of up and down tetrahedra in the lattice separately. Fig. 6 shows the extrapolation to the exact limit for the 64 site cluster. The extrapolation clearly suggests an imbalance between up and down tetrahedra, and confirms the finding from the real space dimer correlations. This is further corroborated by a high susceptibility towards inversion symmetry breaking perturbations, as discussed in the main text.