Sublattice Pairing in Pyrochlore Heisenberg Antiferromagnets

We argue that classical pyrochlore Heisenberg antiferromagnets with small further-neighbor couplings can order in a state where pairs of sublattices form antiparallel spirals. The spiral ordering wave vectors of the two pairs are in general different from each other, and are constrained by which sublattices are being paired. This sublattice pairing state generally breaks inversion and most rotation symmetries. Its existence depends on the antiferromagnetic nearest-neighbor coupling which favors the spins on each tetrahedron to sum to zero. To substantiate our argument, we extend the nematic bond theory; a diagrammatic large-$N_s$ method, to non-Bravais lattices, and we demonstrate that the predicted state is indeed realized at low temperatures in a large region of exchange coupling space. We also carry out a spin wave calculation which suggests that the sublattice pairing state is coplanar.


I. INTRODUCTION
The Heisenberg antiferromagnet on the pyrochlore lattice has gotten much attention as it is a spin liquid candidate. This is mainly motivated by the antiferromagnetic (AF) nearest-neighbor classical model, which is predicted to be disordered at all temperatures [1][2][3][4]. However, real pyrochlore magnetic materials are seldom described by the nearest-neighbor model alone. It is therefore important to understand the effects of further-neighbor couplings, and when and what magnetic order they may cause.
In this article we propose a new kind of ordered state for pyrochlore Heisenberg antiferromagnets: a sublattice pairing (SLP) state, where sublattices pair up, and each pair form antiparallel spirals.
Ordering transitions as a result of adding furtherneighbor couplings has been studied in mean-field theory [2], and it is known that further-neighbor interactions induce symmetry breaking in the purely classical J 1 -J 2 model [5][6][7][8].
The third nearest-neighbor couplings are known to be important for several pyrochlore materials, and in many cases more important than the second nearestneighbor couplings [9][10][11]. There are two inequivalent third nearest-neighbor couplings on the pyrochlore lattice: J 3a which goes in the direction of J 1 , and J 3b , which goes across the hexagons, see Fig. 1. Existing theoretical works including third nearest-neighbor couplings either treat the two as equal or set J 3b to zero when studying ordering transitions [5,12,13].
In this article we treat the Heisenberg antiferromagnets classically. We focus on the hierarchy of magnetic scales J 1 > J 3b ≥ J 2 , J 3a , which might be important for the Gd 2 B 2 O 7 pyrochlores [9] and the AFe 2 O 4 spinels [11].

II. SUBLATTICE PAIRING
An AF nearest-neighbor coupling favors the spins on each tetrahedron to sum to zero [1,2], i.e. 3 i=0 S R,i = 0 for the up-tetrahedra and 3 i=0 S R− ai,i = 0 for the down-tetrahedra. If each sublattice orders in a singleq spiral state S R,i = u i cos Q i · R + v i sin Q i · R, the Fouriertransformed condition for the up-tetrahedra gives where the δ's should be understood modulo a reciprocal lattice vector. This is satisfied by what we refer to as an SLP state. In an SLP state the sublattices form pairs, such that each sublattice pair (i, j) shares the same ordering wave vector Q (i,j) and has antiparallel spins: where u (i,j) and v (i,j) are orthonormal vectors. If the two ordering wave vectors are different, this state satisfies also the condition for the down-tetrahedra if for both pairs of sublattices. Figure 2 shows the planes in momentum space where Q (0,1) and Q (2,3) satisfy this equation.
The ordering wave vectors Q (i,j) are generally found by minimizing the energy. As the SLP state minimizes the J 1 terms of the energy, it is sufficient to minimize the further-neighbor energy terms subject to the condition Eq. (5). As a first example, we consider the pure J 1 -J 3b model. The third nearest neighbors on the pyrochlore lattice couple sites from the same fcc sublattices, and J 3b alone effectively reduces each of the four fcc sublattices to a set of decoupled parallel triangular planes. An AF J 3b will then favor 120 • order in each triangular plane. For a single plane, there are two such chiral ordering vectors, given by clockwise and counterclockwise 120 • rotations. Since the triangular planes in a set are decoupled, the addition of any wave vector orthogonal to the triangular planes will still give 120 • order in each plane, but with an additional inter-plane rotation. This gives rise to a set of lines in momentum space for each of the four sublattices, along which the J 3b -part of the energy is minimal. This is illustrated in Fig. 2. These lines would correspond to rods of scattering if J 1 = 0. The lines intersect at points where the J 3b energy of two sublattices is minimized by the same Q vector. These wave vectors are given by (0, 4π/3, 4π/3) and symmetry-related vectors and lie also on the planes satisfying the tetrahedron condition Eq. (5). The example shown in Fig. 2 is Q (0,1) = ±(0, 4π/3, −4π/3) and Q (2,3) = ±(0, 4π/3, 4π/3), and the corresponding SLP 120 • configuration is illustrated in Fig. 1.
To study the ordering wave vectors for the more general J 1 -J 2 -J 3a -J 3b model, we make use of the fact that this model can be recast into aJ 1 -J 3a -J 3b model when the tetrahedra conditions are satisfied, withJ 1 = J 1 − J 2 andJ 3a = J 3a − J 2 [7]. Q (i,j) is then found as the wave vector that satisfies Eq. (5) and minimizes theJ 3a -J 3b energy sum of the pairing fcc sublattices i and j. For −3 ≤J 3a /J 3b ≤ 1 with AF J 3b , we find that Q (i,j) is given by the vectors symmetry related to with h=2 arccos −(1 +J 3a /J 3b )/2 , that satisfy Eq. (5).

III. NEMATIC BOND THEORY
In order to investigate the occurrence of SLP states in the pyrochlore Heisenberg model, Eq. (1), we employ the nematic bond theory (NBT) [15]. The NBT is a large-N s approximation, where N s is the number of spin components, leading to a set of self-consistent equations for classical Heisenberg magnets. It has previously been employed to the square, cubic and triangular lattices [15][16][17][18]. In this article, we extend the NBT to non-Bravais lattices with m sublattices. Consequently, quantities like the exchange coupling J q and the self-energy Σ q become m × m matrices in sublattice space.
In momentum space, the Hamiltonian is where the q-sum goes over the first Brillouin zone, and i and j are sublattice indices. In the NBT, the classical spins are integrated out by introducing a constraint field λ R,i ensuring unit length of the spins, | S R,i | = 1, through The remaining integrals over the constraint field are treated separately for the average constraints ∆ i and the fluctuations around the average. The average constraints are treated by the saddle-point approximation, and the fluctuations are treated through diagrammatic perturbation theory (large-N s approximation). As in Ref. 15, the diagrams consist of solid and wavy lines, representing spin and constraint propagators, respectively. In addition to sublattice indices, we have here also included directions to the spin propagators to allow for breaking of inversion symmetry. The spin propagator K −1 q,ij is then to be understood as carrying momentum q from i to j.
The spin and constraint propagators are renormalized by the self-energy Σ q and the polarization Π q , respectively, through the Dyson equations in Fig. 3. The large-N s approximation is performed through a set of selfconsistent equations for the self-energy and the polarization, which are shown diagramatically in tions are [15] where K −1 q is the renormalized spin propagator and D q is the renormalized constraint propagator.
As the constraint field now has a sublattice index, we get a separate saddle-point equation for each sublattice where V = L 3 is the number of unit cells. These m saddle-point equations give the temperature T . They must all give the same temperature for the solution to be physical. Following the derivation in Ref. 17, the free energy per unit cell (excluding vertex corrections) is The self-consistent equations are solved by iteration starting from a random self-energy and equal values of the ∆ i s. Each iteration gives an overall negative contribution to the self-energy. To avoid the general increase in temperature associated with this, the ∆ i s are renormalized in every iteration by subtracting from them the minimum eigenvalue among all Σ q . In addition, each ∆ i is adjusted very slightly so that Eq. (12) gives the same value of the temperature for all sublattices. We iterate until the temperature has converged, and then employ K −1 q , Σ q and D q to calculate the free energy. For each initial value of the ∆ i s, we thereby obtain f and K −1 q with a corresponding T .
For a random initial self-energy the NBT might not converge to the lowest temperatures. In those cases we initialize the iterations using a guessed form of K −1 q with peaks at suitable momenta. If different initial conditions converge to different states, we pick the state with the lowest free energy.
To get information about the spin correlations we calculate the quantity A q is periodic with twice the reciprocal lattice vectors and the associated extended Brillouin zone (EBZ) is a truncated octahedron with dimensions twice those of the first Brillouin zone (1BZ) of the fcc lattice.
A q is closely related to the spin structure factor is manifestly inversion symmetric, A q is not as it reflects the symmetries of the self-energy Σ q,ij . We take lack of inversion symmetry in A q to indicate that the spin state breaks inversion symmetry.

IV. RESULTS
For the pure AF nearest-neighbor Hamiltonian, NBT gives no symmetry breaking down to the lowest temperature studied (T 10 −9 ), and A q shows O h symmetric extended maxima on the square surfaces of the EBZ with pinch points at X EBZ [19].
For the J 1 -J 3b model with J 3b = 0.2, the maxima of A q at high temperature occurs at W EBZ . As the temperature is lowered, these maxima move into the hexagonal EBZ surfaces keeping the full O h symmetry. Then at T c = 0.195 (for L = 36) the NBT free energy reveals a first-order phase transition into a low-temperature phase with a total of eight peaks in A q in the EBZ [20]: at Q (0,1) + b 1 + n 2 b 2 + n 3 b 3 for n 2 , n 3 ∈ {0, 1}, and at Q (2,3) + b 2 + n 1 b 1 , Q (2,3) + b 3 + n 1 b 1 for n 1 ∈ {0, 1}, with Q (0,1) = (0, −4π/3, 4π/3) and Q (2,3) = (0, 4π/3, 4π/3). b i denote the reciprocal lattice vectors for the fcc Bravais lattice. The symmetry of A q is thus reduced from O h to C 2v . In particular, inversion symmetry and all three-and four-fold symmetries are broken. When adding also A − q to obtain the spin structure factor, the C 2v symmetry is increased to D 2h . The peaks in A q can be explained as originating from a 120 • SLP state where antiparallel spirals on sublattices 0 and 1 order at Q (0,1) , and antiparallel spirals on sublattices 2 and 3 order at Q (2,3) . We find that such first-order transitions into the C 2v symmetric 120 • SLP phase occur for all positive values of J 3b , see Fig. 5. We next check the stability of the SLP phase when adding J 2 . The finite-temperature phase diagram obtained using NBT for J 3b = 0.2 is shown in Fig. 6. The SLP phase is stable in the region −0.2 ≤ J 2 < 0.372, and the ordering wave vectors follow Eq. (6). For FM J 2 < −0.2, the SLP state becomes unstable to a doubleq state, reminiscent of the multiq states investigated in Refs. 6-8, and 21, where two ordering vectors are present on all sublattices.
For J 2 = −0.2, we find a special case of the SLP state, labeled SLP-X in Fig. 6, where all sublattices have the same ordering vector X 1BZ . A q has maxima at opposite corners of four of the EBZ square surfaces; (2π, 4π, 0), (2π, 0, 4π) , and subleading peaks with half maximum intensity at the four points (0, ±2π, ±2π), consistent with the SLP ordering vector Q (0,1) = Q (2,3) = (2π, 0, 0). The peak locations transform into each other by the D 4h subgroup of O h . Thus, this phase is inversion symmetric, as opposed to the general SLP phase. This SLP-X phase extends both into the doubleq and the general SLP regions at finite temperatures, and will therefore generally cause two ordering transitions as the temperature is lowered from the disordered phase, first one into the SLP-X phase and then another into the doubleq or general SLP phase at a lower temperature, see Fig. 6.
For J 2 = J 3b at low temperatures, when biased into it, NBT converges to an SLP-(0, π, π) state where A q also has D 4h symmetry. In this state the two spirals each reduce to a collinear configuration, one with ordering wave vectors Q (0,1) = ±(0, −π, π) and the other with Q (2,3) = ±(0, π, π). This state exists up to a finite temperature, but not all the way up to the disordered phase. We have not found the proper state at the intermediate temperatures, as we have not been able to get NBT to converge in the shaded region in Fig. 6. Nevertheless, we have indicated phase boundaries around the SLP-(0, π, π) phase as it necessarily must be separated by phase transitions from the inversion symmetry-breaking SLP phase surrounding it.
For sufficiently strong AF J 2 , the SLP phase gives way to a singleq state with an ordering wave vector along the ΓX 1BZ line. In this phase all sublattices order at the same wave vector, and the tetrahedra conditions are no longer satisfied. The phase boundary between the SLP phase and the ΓX 1BZ phase is estimated to lie between J 2 = 0.366 and J 2 = 0.374 from our NBT calculations. The vertical line at J 2 = 0.372 is chosen from where the minimum of J q changes character.
In Fig. 7 we map out the low-temperature phase diagram in the J 3b -J 2 coupling space using NBT. It is seen that the SLP phase exists in a large region.
On the AF J 2 side, the SLP phase ceases to exist when it becomes energetically favorable to violate the tetrahedra conditions. For small and intermediate values of J 3b we find the singleq ΓX 1BZ phase. For large J 3b the SLP phase is stable up to J 2 = 1, which is the limit set by the mapping from J 1 -J 2 toJ 1 -J 3a .
On the FM J 2 side, the SLP phase borders the doubleq phase where each sublattice has two ordering vectors. The two ordering vectors are (2π, l, l) and (2π, l, −l), where l increases from 0 (X 1BZ ) to π/2 (U 1BZ ) as J 3b decreases from −J 2 . It reaches π/2 for J 3b equal to a small positive J 2 -dependent value. For J 3b less than this, the two ordering vectors shift to (0, k, k) and (0, −k, k) with k decreasing slowly from 3π/2 (K 1BZ ) as J 3b decreases further. Spins in such doubleq spirals will only obey the length constraint when the number of spin components N s > 3. The NBT is derived from the large-N s limit without vertex corrections, and is extrapolated down to N s = 3. We believe that the doubleq state produced by NBT is a remnant of the large-N s limit and that the extrapolation down to N s = 3 does not take the length constraint sufficiently into account.
We have also performed a similar stability analysis in J 3b -J 3a coupling space with J 2 = 0. There, for AF J 3b , the SLP state is stable in an even wider region; whenever J 3a ≤ J 3b .
To get information about the relative orientation of the spiral plane vectors u (i,j) and v (i,j) of the two SLP spirals, we have performed a spin wave calculation to compute the entropy. We find that entropy favors the SLP state to be coplanar for our model, Eq. (1), with the two SLP spirals sharing spiral plane vectors. Consequently, entropy favors collinear states in the special cases of SLP-Γ and SLP-X.

V. DISCUSSION
We have shown how the classical AF Heisenberg model on the pyrochlore lattice with small further-neighbor couplings orders with coplanar SLP. In SLP, pairs of sublattices form antiparallel spirals. The ordering wave vectors are in general different for the two sublattice pairs, and are found to be the wave vectors that minimize the total J 3a -J 3b energy of the paired fcc sublattices subject to the tetrahedra conditions.
For the pure J 1 -J 3b model, we find a 120 • SLP state. This state simultaneously satisfies both AF J 1 and AF J 3b and is thus realized for all J 3b > 0. It is separated from the disordered phase by a first-order phase transition with a T c and latent heat that goes to zero as J 3b → 0, Fig. 5. This is consistent with the AF nearestneighbor model not ordering [1][2][3][4]. Note that the critical temperatures are likely to be overestimated as NBT excludes vertex corrections [16].
The SLP state generally breaks inversion symmetry, except when four times the ordering wave vectors are reciprocal lattice vectors. We note that recent numerical results indicate that quantum fluctuations of the purely AF spin-1/2 and spin-1 models also induce inversion symmetry-breaking [22][23][24].
As special cases of the SLP state, where all sublattices order at the same wave vector, we find SLP-Γ (Néel) and SLP-X for our model. SLP-Γ (Néel) has previously been identified as the ground state for the J 1 -J 2 model with small AF J 2 [7,21,25]. SLP-X is realized along the lineJ 3a = J 3b > 0 (and close to this line at intermediate temperatures) and should be the symmetrybroken state for the J 1 -J 3 model in Ref.
12. An ab inito study of the breathing pyrochlore material LiGaCr 4 O 8 findsJ 3a = J 3b > 0 and SLP-X as the corresponding low-temperature state [26]. They also find SLP-X to be stabilized at intermediate temperatures for LiInCr 4 O 8 , whereJ 3a ≈ 1.8J 3b > 0, which could be explained by the finite-temperature extension of the SLP-X phase due to its collinearity.
The J 2 and AF J 3a bonds favor states which do not satisfy the tetrahedra conditions. Nevertheless, we find that the SLP state dominates a large portion of the exchange coupling space, particularly in the region J 1 > J 3b >J 3a for AF J 1 and J 3b . This region might be relevant for the Gd 2 B 2 O 7 class of materials (B is a nonmagnetic cation) [9]. These are believed to be described as classical pyrochlore Heisenberg AFs with further-neighbor and dipole-dipole inter-actions [27][28][29]. While we have not considered the dipoledipole interaction, we note that a recent experiment on Gd 2 Ti 2 O 7 suggests a partially ordered state with two ordering vectors at different L 1BZ points [30]. For SLP states such L-peaks are most likely to occur at intermediate temperatures near the line J 2 = J 3b > 0 where there are line minima in-between the L 1BZ points.
The pyrochlore spinel materials are expected to have J 3a at the same order of magnitude as J 3b [10,11]. As we find the SLP phase to be stable for J 3a ≤ J 3b for AF J 3b , it could be of relevance also for this class of materials. Especially in the ferrites AFe 2 O 4 , the two types of third nearest-neighbor couplings have been suggested to be of comparable strength [11]. Inclusion of an AF J 2 could also help to push the system towards the SLP phase.
We envision also the general SLP states to be realized on the breathing pyrochlore lattice when both the "nearest"-neighbor couplings J 1 and J 1 are sufficiently strong AF, such that the tetrahedra conditions are satisfied. For further work it would be interesting to study the stability of the SLP state on the breathing pyrochlore lattice as well as its stability when adding dipolar interactions and/or anisotropy. We also note that the extension of NBT to the pyrochlore lattice can be used to study other interesting problems such as possible symmetry breaking in the J 2 = J 3a model [13].