Magnetization process and ordering of the $S=1/2$ pyrochlore magnet in a field

We study the $S=1/2$ pyrochlore Heisenberg antiferromagnet in a magnetic field. Using large scale density-matrix renormalization group (DMRG) calculations for clusters up to $128$ spins, we find indications for a finite triplet gap, causing a threshold field to nonzero magnetization in the magnetization curve. We obtain a robust saturation field consistent with a magnon crystal, although the corresponding $5/6$ magnetization plateau is very slim and possibly unstable. Most remarkably, there is a pronounced and apparently robust $\frac 1 2$ magnetization plateau where the groundstate breaks (real-space) rotational symmetry, exhibiting oppositely polarised spins on alternating kagome and triangular planes. Reminiscent of the kagome ice plateau of the pyrochlore Ising magnet known as spin ice, it arises via a much more subtle `quantum order by disorder' mechanism.

Introduction.-The Heisenberg antiferromagnet on the pyrochlore lattice, is one of the most frustrated three dimensional magnets, and as such a prime candidate for exotic, specifically quantum spin liquid, behavior. Indeed, its classical variants are topological magnets (classical spin liquids), exhibiting a Coulomb phase [1][2][3] in the limit of low temperature both for the Heisenberg [4] and the Ising (spin ice) variants, with the latter also hosting deconfined magnetic monopoles as quasiparticles [5][6][7][8].
The classical variants exhibit either no magnetization plateau for the Heisenberg model (at least in the absence of magnetoelastic coupling [38][39][40]), or, in the Ising case, the very rich physics of kagomé ice [41][42][43]. The latter arises for a field applied in a [111] direction, which has a large projection onto the local Ising axes on one quarter of the spins, which therefore quickly get polarised. The remaining three quarters of the spins reside on kagomé * hagymasi@pks.mpg.de † schaefer@pks.mpg.de ‡ moessner@pks.mpg.de § dluitz@pks.mpg.de Normalized ground state magnetization m = ψ0 |S z tot |ψ0 of different pyrochlore clusters as a function of the external magnetic field h. The saturation magnetization is given by the fully polarized state of N spins: msat = N/2. For the 64 site cluster the upper half, while for the largest clusters 108 and 128, only the part of the curve at the strongest fields can be reliably calculated. Fig. 2) which enter a stable partially polarised plateau at intermediate field strengths.
Here, we study the magnetization process of the S = 1/2 Heisenberg antiferromagnet from zero field to saturation with DMRG, which has recently turned out to be FIG. 2. Onsite magnetization for the m/ms = 1/2 plateau of the 108 site and 128 site clusters. We show only the cubic unit cell cut of the clusters for better comparison. Green balls correspond to polarized spins in one direction, red balls represent a polarization in the opposite direction. The size of the balls is proportional to the onsite magnetization. Gray shaded triangles are guiding the eye to indicate the polarized kagomé planes.
very useful regarding the zero-field properties [19]. Most saliently, we find an incompressible magnetic phase with a 3:1 spin polarization ratio, signaled by a robust plateau at half saturation.
The ground state corresponding to this plateau breaks the rotational symmetry of the lattice. It exhibits polarized kagomé planes along the field direction and antipolarized interplane sites. Unlike in kagomé ice, this pattern of 3:1 disproportionation arises spontaneously, being selected from an exponentially large number of possible 3:1 disproportionations in what may be termed a quantum order by disorder mechanism. Also, the minority sublattice has negative, rather than enhanced, Zeeman energy -a spontaneous instance of ferrimagnetism.
Furthermore, we obtain a value of the saturation field which is consistent with what one obtains for the magnon crystal state, an exact eigenstate of a range of frustrated Heisenberg magnets [25].
Methods.-We determine the field-dependence of the ground state magnetization for periodic clusters ranging from N = 32 to N = 128 spins, using SU(2) and U(1) DMRG [52][53][54][55]. Although DMRG is by nature a one dimensional method [56][57][58][59][60], it has been successfully used in two [61] and recently three dimensions [19,62], by 'linearizing' the system along a snake path at the price of non-local interactions within the snake. Conservation of the total spin, S tot and its z-component, S z tot allows us to target and optimize the ground states within the concomitant symmetry sectors. We observe that despite the fact that the SU(2) representation is very efficient in higher spin sectors, the convergence while optimizing the energy is sometimes better when explicitly enforcing the U(1) spin symmetry only. Since the wave function is represented as a matrix-product state with finite bond dimension, extrapolation to infinite bond dimension is necessary. This is usually done by either extrapolating as a function of the truncation error or the variance [59].
(2, 0, 0) T (0, 2, 0) T (0, 0, 2) T 8 We optimize the wave function for small bond dimensions ( 2000) using the two-site DMRG algorithm, but for larger bond dimensions we switch to the single-site variant with subspace expansion [54]. For the bond dimensions we used (up to ∼ 20000 SU(2) or U(1) states) the calculation of the full variance is impractical, and we therefore extrapolate these energies to the error free limit as a function of the two-site variance. This quantity was shown to be equally good compared to the truncation error [63], see Appendix A for further details. Magnetization curve.-We consider the ground state of the S = 1/2 pyrochlore Heisenberg antiferromagnet in a finite magnetic field h > 0. The spins are located on a pyrochlore lattice defined by the fcc lattice vectors a 1 = 1 2 (1, 1, 0) T , a 2 = 1 2 (1, 0, 1) T , a 3 = 1 2 (0, 1, 1) T , together with four tetrahedral basis vectors b 0 = 0, b i = 1 2 a i . This model conserves the total magnetization S z tot = i S z i , since [H, S z tot ] = 0, and hence the Hamiltonian decomposes into symmetry sectors with fixed total magnetization m = −N/2, −N/2+1 . . . N/2. This means that each eigenstate |n m of H is also an eigenstate of S z tot : S z tot |n m = m|n m and therefore the eigenstates of the Hamiltonian are independent of the field h, and their In the absence of a field, h = 0, the ground state of the Hamiltonian is in the m = 0 sector. For finite values of h > 0, the energies of states exhibiting a finite magnetization m = 0 will change by −hm and potentially become the overall ground state of the system. This leads to the characteristic jumps of the magnetization in Fig. 1. The field strength at which a transition of the ground state magnetization from m to m + 1 occurs is entirely determined by the minimal energy of all states in the sectors m and m + 1 in the absence of the field, E 0 . Multiple transitions between m, m + 1, m + 2 . . . can coincide, leading to a larger change at a given field strength. To obtain the full magnetization curve for a given cluster, we therefore have to calculate the lowest energies in all magnetization sectors at zero field using DMRG, respecting the U (1) symmetry associated with the conservation of S z tot . In fact, since the total spin S 2 tot also commutes with both H and S z tot , we can also use the full SU(2) symmetry as mentioned above.
If the ground states of adjacent S z tot sectors are separated by a gap in the thermodynamic limit, the magnetization will remain locked into the lower magnetization sector for a range of fields proportional to the gap, leading to a plateau in the magnetization curve, and hence an incompressible magnetic state.
We display the resulting magnetization curves in Fig. 1 for different clusters ranging from 32 to 128 sites, noting that for large clusters it is only possible to determine the large field part of the magnetization curve.
A prominent feature is the large jump of the magnetization from its maximum to 5/6 m sat near the saturation field h sat = 4J. Position and height of this jump can be understood using the same reasoning as in the case of kagomé [28,32]. The exact ground state of symmetry sectors with very large S z tot is a crystal of localized magnon modes, while the ground state of the sector with maximal S z tot is the trivial fully polarized state. One possible densest packing of independent magnon modes on the pyrochlore lattice is given by densely arranging magnons localized on hexagonal motifs in the kagomé planes, leading to a magnetization plateau at m/m sat = 5/6 (see Appendix B for details). Each magnon mode reduces the energy by 4J, so that h sat = 4J (Eq. (B4)). Each independent magnon excitation requires three unit cells (twelve sites), which fixes the densest packing. Due to the requirement of commensurability of the densest packing with the cluster geometry, we find the 5/6 plateau only for the 108 site cluster in Fig. 1. On other clusters we can obtain even larger jumps at h = 4J, as modes localized on shorter loops winding across periodic boundary conditions can yield a denser magnon mode filling as a finite-size effect. This yields, for the 32 and 48 site clusters, a broad plateau at m/m sat = 3/4 which is not representative for the infinite lattice. The narrowness of the 5/6 plateau on the 108 site cluster in turn calls into question its stability in the thermodynamic limit.
We attempt to extrapolate the widths of the magnetization plateau observed in finite size clusters ( Fig.  1) to the thermodynamic limit (Fig. 3), using linear fits in 1/N . There is little indication that the 3/4 and 5/8 plateau retain a finite width. In contrast, we have strong evidence for a finite width of the 1/2 plateau in the thermodynamic limit (red in Fig. 3), which is located in the field range h  (11) agrees with that of the recent variational Monte Carlo estimate, 0.40(4), in the thermodynamic limit [21] and the data shown in Fig. 3 is consistent with a finite triplet gap.
It is worth commenting on the actual values of the magnetic field to make the comparison with experiments easier. Since a material, which realizes the S = 1/2 Heisenberg model is still lacking, the closest relative we can consider is the S = 1 NaCaNi 2 F 7 [51]. Assuming the same value of J ∼ 3.2 meV in a S = 1/2 material (and g factor g ∼ 2), the saturation field would correspond to B sat ∼ 110 T and the 1/2 plateau is expected to start at ∼ 68 T. Such high values of magnetic fields are accessible in pulsed field experiments.
Properties of the 1/2 plateau.-We turn to the correlations of the 1/2 plateau. We focus on the largest cluster, where finite-size effects due to short loops winding across periodic boundaries are suppressed. Nonetheless, we include results for smaller clusters for finite-size extrapolations, as in the preceding analysis. Indeed, while the 32 site cluster develops a uniform magnetization throughout the lattice with ψ 0 |S z i |ψ 0 ≡ 0.25, larger clusters exhibit a more complex pattern, stabilizing inequivalent spins with differing polarization. These are arranged with respect to one of the four rotational axes defined within each tetrahedron. The planes perpendicular to this axis are alternating triangular and kagomé planes, containing 1/4 and 3/4 of the spins, shown in red and green in Fig. 2, respectively: each tetrahedron contributes a 'base triangle' to the kagomé, and its apex spin to the triangular, plane.
The kagomé spins (A sites in Tab. II) acquire a polarization along the field, while the triangular spins (B sites) are polarised in the opposite direction in our clusters of size 64, 108 and 128. Fig. 2 shows the onsite magnetization pattern ψ 0 |S z i |ψ 0 in a cubic unit cell for the clusters of size 108 and 128, along with a list for different clusters in Tab. II. While the largest cluster with 128 sites develops this pattern perfectly, smaller clusters can exhibit defects in the form of lines with small onsite magnetization passing through kagomé planes, which we attribute to the existence of short resonant loops across the periodic boundaries in these clusters. The number of such defective sites is listed in the caption of Tab. II and vanishes for the 128 cluster. We emphasize that this symmetry breaking is distinct from the one we found for the zero-field problem, where the lattice inversion symmetry appears to be broken, as the two sublattices of tetrahedrons have different energy densities [19]. On the 1/2 magnetization plateau, the inversion symmetry is preserved, while the rotational symmetry of the lattice is broken by the emergence of a preferred [111] axis. The rotational symmetry around this axis is not broken, but those within the three other kagomé planes are.
The symmetry breaking is naturally reflected in the spin structure factor where R i denote the real-space coordinates of sites and the index c denotes the connected part of the correlation matrix (the factor 4/3 comes from normalization 1/(S(S + 1)) for spin S = 1/2). This is plotted in Fig. 4 for several clusters. The larger clusters exhibit clearly discernible bright lines in the structure factor, discarding the rotational symmetry present for the 32-site cluster; in contrast to the m/m s = 0 state, the inversion symmetry remains intact. Earlier work [44] predicts another possible pattern, the R-state, for the S = 3/2 case. We discuss the competition of the R-state with the best variational wave func- tion in Appendix A and conclude that the R-state has higher energy. Discussion -The pyrochlore Heisenberg antiferromagnet in a field, like its zero-field cousin, illustrates the capacity of frustrated magnets to exhibit a plethora of instabilities, discarding the rotational symmetry at halfmagnetization and forming an incompressible state. It is interesting to embed this finding in a broader context.
Firstly, the idea that the half-magnetization plateau goes along with a 3:1 disproportionation of sites seems entirely natural: indeed, for collinear spins, such a 3:1 ratio is the only way to obtain half-magnetized tetrahedrons. Note, however, that such 3:1 tetrahedrons can be tiled in exponentially numerous ways on the pyrochlore lattice, indeed mapping onto dimer coverings of the diamond lattice, which have a finite entropy of S ≈ 0.13k B [64][65][66]. The selection through fluctuations of a specific one (or subset) of such tiling with lowered symmetry is known as order-by-disorder. In this sense, our magnetization plateau exhibits a form of quantum order-bydisorder, although the assignment of what term of the Hamiltonian contributes the fluctuations is to some degree a matter of choice. At any rate, the emergence of ordered ferrimagnetism in the plateau is a striking instability of a highly frustrated quantum magnet.
We close by contrasting the 1/2 plateau to the kagomé ice plateau mentioned in the introduction. In a real pyrochlore material with spin-orbit coupling, an Ising anisotropy needs to go along with non-collinear easy axes: the local easy axis is the [111] direction joining a site with the centers of its tetrahedrons. This turns a uniform applied magnetic field into a staggered Zeeman field according to the sublattice [67]; as mentioned above, when applied along a [111] direction, it polarizes the triangular planes more strongly than the kagomé ones, as the easy axes projections differ by a factor of 3. We find it most intriguing that this general setting arises for the S = 1/2 Heisenberg plateau by spontaneous rather than explicit symmetry breaking, and the question of interpolating between the two immediately poses itself. Note that the two cases differ considerably in (non-symmetry) 'details': the triangular layers are strongly positively polarised in kagomé ice, in contrast to their weak negative polarization in the Heisenberg S = 1/2 plateau.
There clearly remains much further scope for studying the highly frustrated quantum magnets in d = 3 with and without applied fields, and for the foreseeable future, recent technological progress promises to yield previously inaccessible interesting data such as those underpinning the present article.

ACKNOWLEDGMENTS
We acknowledge support from the Deutsche Forschungsgemeinschaft through SFB 1143 (project-id 247310070) and cluster of excellence ct.qmat (EXC 2147, project-id 390858490). I.H. was supported in part by the Hungarian National Research, Development and Innovation Office (NKFIH) through Grants No. K120569 and No. K134983. Some of the data presented here was produced using the SyTen toolkit [52,53].

Appendix A: Further details of the DMRG simulation
Matrix-product operators (MPOs) can be constructed by hand for the simplest 1D models with nearestneighbor interactions. However, this task becomes difficult when long-range interactions are present and needs to be done in an automatic way. While the corresponding MPO of a single product, e.g. S z i S z j , can be represented by an MPO of small bond dimension, the bond dimension grows exponentially by summing up multiple terms. Luckily, the MPO can be compressed effectively using the deparallelisation method [68] without any information loss in many cases. We start by optimizing a random matrix-product state (MPS) with the corresponding size of the cluster using DMRG where the long-range correlations are captured automatically up to a given bond dimension. The choice of the 3D→1D mapping can have an influence to the overall convergence as well as the final bond dimension of the MPO, which can also have an impact on the speedup of the calculation. However, we did not observe a significant difference in terms of computation time and convergence properties for the different paths in the three dimensional clusters we investigated.
The convergence problems in three dimensions are even more severe than in two dimensions, since there are more periodic bonds yielding to increasing long range interaction. While we do not face convergence problems for the 32 and 48 site clusters, DMRG often exhibits bad convergence and gets stuck in local minima for the larger clusters due to its local-update nature. In these cases, the pattern in Fig. 2 is only partially generated if we initialize the algorithm with random states. This can be monitored either by examining the truncation error or the two-site variance as a function of the bond dimension. Metastable states induce a non-monotonic behavior of these quantities, that is, the energy decreases but the truncation error or two-site variance increases. When the convergence is smooth, the energy follows typically a linear behavior as a function of the two-site variance [63].
To stabilize one of the magnetically ordered states in our simulation, we therefore use the pinning-field technique [61]. We apply a magnetic field (at low bond dimensions) to the same set of sites for each tetrahedral unit cell which compatible with the polarized kagomé planes such that the ordered state is stabilized. We then remove this pinning field and continue the optimization procedure while further increasing the bond dimension. With this approach the overall convergence becomes much better and smoother as it is indicated by Fig. 5.
To assess the quality of the variational ground state, we compare the total energies at a fixed bond dimension χ = 8000 for the 64 and 108 site clusters using different initial conditions: starting from random initial states or using the pinning-field technique to start from ordered patterns. The ordered pattern produces lower energies (∼ 0.1J) than starting from random initial states. We also consider another possible pattern, the R-state [44], which was predicted for the S = 3/2 case. This state also yields higher energies than our best variational state, although only with a small difference ∼ 0.05J and is therefore clearly in the low energy manifold. The error of the extrapolated energies are defined as one quarter of the distance between the lowest energy DMRG datapoint and the extrapolated value, which is commonly used in the community [69]. The same definition of the error bars is used in the main text as well and should be understood as an estimate of the systemic extrapolation error.

Appendix B: Localized magnons in high fields
Constructing analytic solutions for quantum many body systems is a challenging discipline and succeeds only in special cases. Therefore, solutions of the ground state in the form of quasi particles in high fields were a huge success for spin systems in certain frustrated lattices [23][24][25][26][27][28][29]. Probably the most famous example for these quasi particles are independent and localized magnon excitations in the two dimensional kagomé lattice.
Kagomé lattice Localized magnon excitations describe the ground state of the Heisenberg model in a large external field on the kagomé lattice. Theses are confined to non-touching hexagons such that the description can be limited to a single star of David. Starting from the fully polarized state | ↑ ... ↑ , a single magnon state is given by (up to normalization) where j runs over the hexagon. The localization can be easily verified since each corner spin of the star is attached to two spins of the inner hexagon. The contributions of flipped spins propagating to corner sites are canceled due the alternating sign. Hence, the magnetic excitation remains within the hexagon and preserves the alternating sign structure such that the hopping contribution is diagonal. Not only are the magnons localized but also they are ground states of the Heisenberg antiferromagnet for high fields. For simplification we set h = 0 and focus on the hopping H ± and interaction term H z individually: Within the hexagon, the sign of each term in Eq. (B1) is inverted by H ± . Hence, H ± |m = −J|m and the hopping term reduces the energy by J. In the kagomé lattice, each site is attached to 4 other sites. The contribution to the energy by the interaction H z due to single spin flip is a reduction of 2J in contrast to the fully polarized state. In total, a single magnon reduces the energy by 3J. Due to the localization, multiple independent states can be placed within the kagomé lattice as long as they are separated by at least one site. In this manner, each hexagon is uniquely assigned to 9 sites in the kagomé lattice to achieve the densest filling. The complete tiling of hexagons describes the ground state and corresponds to a magnetization of m/m sat = 7/9. The energy per site is reduced to E 7/9 = 1/6J from the fully polarized state E 1 = 1/2J. Pyrochlore lattice The tetrahedral unit cell of the pyrochlore lattice consist of four sites. Each face of this tetrahedron defines one out of four orientations of parallel kagomé planes in the lattice. The apex spins form a separating triangular plane between the kagomé planes. Equivalently to the two dimensional case, the same formalism can be used to generate localized magnons in the pyrochlore lattice. As visualized in Fig. 6 by the red dotted triangles, 9 sites are uniquely assigned to each localized hexagon in the kagomé plane (black sites), such that no supercells are sharing the same site. The corresponding magnon excitation is illustrated by the red circles. Red and blue sites correspond to the upper and lower separating triangular layer respectively. In addition to the 9 sites laying inside the plane, we include the