Theory of antiferromagnetic Heisenberg spins on breathing pyrochlore lattice

Spin-singlet orders are studied for the antiferromagnetic Heisenberg model with spin $S$>1/2 on a breathing pyrochlore lattice, where tetrahedron units are weakly coupled and exchange constants have two values $0<J' \ll J$. The ground state has a thermodynamic degeneracy at $J'$=0, and I have studied lattice symmetry breaking associated to lifting this degeneracy. Third-order perturbation in $J'$ for general spin $S$ shows that the effective Hamiltonian has a form of three-tetrahedron interactions of pseudospins $\tau$, which is identical to that previously derived for $S$=1/2, and I have calculated their matrix elements for general $S$. For this effective Hamiltonian, I have obtained its mean-field ground state and investigated the possibility of lattice symmetry breaking for the cases of $S$=3/2 and $1$. In contrast to the $S$=1/2 case,$\tau$'s response to conjugate field has a $Z_3$ anisotropy in its internal space, and this stabilizes the mean-field ground state. The mean-field ground state has a characteristic spatial pattern of spin correlations related to the lattice symmetry breaking. Spin structure factor $S(q)$ is calculated and found to have symmetry broken parts with amplitudes of the same order as the isotropic part.


Introduction
Frustrated magnets are a playground of the experimental and theoretical studies for the quest of new quantum phases [1,2]. Thermodynamic degeneracy of the classical ground-state manifold is the most important ingredient, and the question is how this degeneracy is lifted to select a unique quantum ground state if it exists. Several cases still select some types orders of magnetic dipole moments despite frustration effects, but there also exist three other cases from the viewpoint of symmetry breaking: (i) Spin rotation symmetry is broken, but the order parameter is not a conventional magnetic dipole but something more exotic like quadrupole or vector chirality. (ii) While spin rotation symmetry is not broken, another type of symmetry is broken like lattice symmetry. (iii) No symmetry is broken: this case corresponds to a spin liquid but the liquid behavior does not determine if spin gap is finite or zero. Valence bond crystal (VBC) state [3,4] is a representative example of the case (ii), and the lattice rotation and/or translation symmetry is broken. Affleck-Kennedy-Lieb-Tasaki (AKLT) state [5] is an example of the case (iii) and the spin gap is finite.
Generally speaking, exotic quantum phases are stabilized in frustrated magnets, if the frustration is strong enough and also if quantum fluctuations are large enough. The first condition is related to lattice geometry, and the Kagomé and pyrochlore lattices are the most typical examples in dimensions two and three, respectively. As for the second condition, factors enhancing quantum fluctuations include a high symmetry of the Hamiltonian, a low spatial dimensionality, and a small value of spin S.
c The Author(s) 2012. Published by Oxford University Press on behalf of the Physical Society of Japan. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The antiferromagnetic Heisenberg model on the pyrochlore lattice is a canonical example of frustrated quantum systems in three dimensions [6][7][8][9][10][11][12][13]. This lattice is a network of corner-sharing tetrahedron units, and thus each unit is highly frustrated. The pioneering mean-field analysis [6] showed a huge degeneracy of the semiclassical ground-state manifold, manifested by the presence of zero-energy excitations in all over the Brillouin zone [14]. Over a decade ago I studied the quantum limit S= 1 2 of this model and examined the possibility of exotic orders like scalar chirality [11]. This expectation came from the fact that the doubly degenerate ground states in each tetrahedron unit have opposite scalar chiralities. I derived an effective Hamiltonian in the subspace where every tetrahedron unit is within the two-dimensional local ground-state multiplet of spin singlet and analyzed that Hamiltonian. The result showed that the ground state does not have a scalar chirality order but it is a mixture of local singlet dimers and/or tetramers and thus the lattice symmetry is broken [11,12]. The spatial pattern of these dimers/tetramers is quite complicated and this is due to frustration in their configuration. Three-quarters of tetrahedron units have a specific favorable configuration of dimer pairs, but the remaining quarter of units have no favorable configuration. This is a frustration in the mean-field level, and quantum fluctuations select a uniform order of either dimer pairs or tetramers in the remaining part [12].
In my previous study for the S= 1 2 case, I introduced one parameter for controlling geometrical frustration. The original pyrochlore lattice was split into two parts: one is the set of pointing-up tetrahedron units and the other is the set of bonds connecting these units. The control parameter is the ratio of exchange constants for the two parts J /J, and I approached the original model (J /J=1) by perturbation starting from the decoupled limit J /J=0 [11,12]. A different split was also examined by another group [10] but my split had the advantage of keeping the tetrahedral lattice symmetry.
A few years ago, Okamoto et al. [15] noticed that Cr ions in the compounds LiGa 1−x In x Cr 4 O 8 (0 ≤ x ≤ 1) constitute a network corresponding to my previous perturbative expansion and named this sublattice breathing pyrochlore. Despite the same lattice structure, this system differs from my previous model in the point that Cr ions have spin S= 3 2 . The In end (x = 1) has largest breathing and shows an antiferromagnetic phase transition at the temperature T N =13-14 K [16,17]. This is much smaller compared with the modulus of the Weiss temperature |θ W |=332 K, showing the effects of frustration. This magnetic phase is very fragile against reducing In-ion concentration, disappearing at around x ∼ 0.9, and the system remains paramagnetic down to the lowest temperature 2 K. This may indicate that nonmagnetic ground state is stabilized in the breathing pyrochlore lattice. Therefore, it is interesting to reinvestigate the problem now for the case of S= 3 2 and also for general S, and I will examine how the results depend on S. It will turn out that for S > 1 2 there appears a generic anisotropy in dimer/tetramer configuration and this stabilizes a specific spatial modulation of spin correlations. This will be manifested in the wave-vector dependence of the energy-integrated spin structure factor S(q), and I will calculate its explicit form.

Model
In this paper, I study the ground state of a spin-S Heisenberg model on a breathing pyrochlore lattice, and apply the results to the S = 3 2 case, which may be realized in the Cr compound LiXCr 4 O 8 . The original pyrochlore lattice is a network of two types of corner-sharing tetrahedra, and they have the same size. In breathing pyrochlore lattice, one type of tetrahedra expand in size and the other type shrink, but neither of them change their shape of regular tetrahedron. See Fig. 1(a). The Hamiltonian to study is a spin-S Heisenberg model with antiferromagnetic interactions between nearest neighbor sites on the breathing pyrochlore lattice, and I will analyze the cases of S = 3 2 and 1 in detail. Exchange coupling has two values depending on bond length in the breathing lattice structure

2/28
where r denotes the position of small tetrahedron and i, j are the site labels as shown in Fig. 1(b). Six bonds in each small tetrahedron unit have strong interaction J > 0, while the units are coupled by long bonds with weak antiferromagnetic interaction 0 < J ( J). In the following, I will call small tetrahedra with strong bonds tetrahedron units, and then they are connected by weak bonds in large tetrahedra. Throughout this paper, I use N to denote the number of spins, and then the number of tetrahedron units is 1 4 N. Note that the Hamiltonian is invariant upon exchanging J and J due to the lattice symmetry.
I studied in Refs. [11,12] this model for the S = 1 2 case and discovered a complex pattern of spin dimers and tetramers in the ground state. I will later focus on the special cases of S = 3 2 and S = 1 afterwards, but first study this Hamiltonian for general S to compare the cases of different spin quantum numbers S.

Spin singlet states in one tetrahedron unit
Following Refs. [11,12], I employ an approach of degenerate perturbation in J /J. The first step is to understand eigenstates in the limit of decoupled tetrahedra (J =0). The Hamiltonian of a single tetrahedron unit at position r is Here, S unit (r) = ∑ 4 i=1 S i (r) is the total spin of this unit and its quantum number is an integer 0 ≤ S unit ≤ 4S. Therefore, the eigenenergies of H unit (r) are determined by the unit spin alone as E unit = J[ 1 2 S unit (S unit + 1) − 2S(S + 1)] and the ground-state manifold coincides with the entire space of S unit =0. These results have been well-known including the fact that the ground states in 3/28 each unit have degeneracy 2S + 1, and thus the ground-state degeneracy in the entire system is (2S + 1) N/4 where N is the number of original spins [11]. This corresponds to the residual entropy that is precisely 25% of the total entropy irrespective of the value of S.
The issue of this study is how the weak inter-unit interactions J release the macroscopic entropy and which type of ground state is selected. Since anything interesting happens in the spin-singlet space, the spin rotation symmetry has no chance to be broken, and it is the lattice symmetry that can be broken. To examine this issue, symmetry argument is useful and I will check how the (2S + 1)-fold ground states in each unit are transformed with operations of the point group symmetry. Each unit has the tetrahedral symmetry T d , and this point group has 5 types of irreducible representations (irreps) [18]: 2 one-dimensional ones (A 1 and A 2 ), 1 two-dimensional one (E), and 2 three-dimensional ones (T 1 and T 2 ).
To perform calculation, we need an explicit form of (2S + 1)-states in the space of S unit = 0. A convenient choice is the following one where l is an integer satisfying 0 ≤ l ≤ 2S. Here, φ (12) (l, l z ) is the wavefunction in which the two spins S 1 and S 2 couple to form a state with the composite spin l and its z-projection l z , while the remaining two spins S 3 and S 4 do the same in φ (34) (l, −l z ) but with the opposite z-projection. These two-spin wavefunctions are written with single-spin bases as where S, S, m, l z − m|l, l z is the Clebsch-Gordan (CG) coefficient 1 of combining two angular momenta [19], and these wavefunctions have the symmetry These two composite spins couple and finally form a total spin singlet in a tetrahedron unit. The prefactor C in Eq. (3) is also given by a CG coefficient and this is simple because the total spin is singlet These (2S + 1) wavefunctions {Φ l } constitute a complete orthonormal set in the S unit =0 space at each tetrahedron unit. I have classified these Φ l 's according to the T d point group symmetry. This was done by calculating the characters of the symmetry operations, and the details of calculation are explained in Appendix A. The result turns out interesting. For any value of individual spin S, all of the 2S + 1 ground states in the tetrahedron unit belong to A 1 , A 2 , and E irreps, while no ground states transform as the three-dimensional irrep T 1 or T 2 . I have calculated the multiplicity of these irreps and the result is Here, χ 2 =mod (2S − 1, 2) and χ 5 =mod (2S − 1, 3) − 1.

Effective Hamiltonian for general S
Now, I am going to derive an effective Hamiltonian that lifts the macroscopic degeneracy of the ground states in the limit of decoupled tetrahedra. To this end, one needs a degenerate perturbation in the weak interaction J , and I succeeded in this task for the S = 1 2 case after a lengthy calculation of many matrix elements [11,12]. For larger spins, the local Hilbert space increases its dimension, and this makes calculations more impracticable and difficult. Therefore, I took a different strategy and tried to simplify the formulation in perturbation as much as possible. I have achieved a huge simplification in the third-order perturbation, and this works for any value of spin S. With this simplified formulation, my previous result for the S = 1 2 case is also easily reproduced. The effective Hamiltonian is to be derived for describing dynamics in the low-energy subspace where all the local states at tetrahedron units are within the spin-singlet manifold, and let me comment on its validity. The use of such an effective Hamiltonian is justified under two conditions. The first condition is the presence of a finite spin gap ∆ s > 0. The size of the spin gap depends on the ratio of two exchange constants 2 , ∆ s = J∆( J J ; S) with∆(0; S) = 1, and the first condition is satisfied at least for small J J . The second condition is that the energy range of consideration should be smaller than the spin gap, ∆ E < ∆ s , where ∆ E is measured from the ground-state energy. For studying the high-energy region ∆ E > ∆ s , it is necessary to take account of the subspaces with S unit ≥ 1 on the same footing as the S unit = 0 subspace, which is beyond the approximation of this effective Hamiltonian.
Before demonstrating how perturbation calculation is simplified, I now introduce matrix elements necessary for that and discuss their symmetry. They are two-spin correlations defined for a pair of ground states in one tetrahedron unit 3 Here, µ, ν ∈ {x, y, z} are spin index, and the result is diagonal in spin space because spin-singlet states are rotationally invariant. For the same site correlation, it is trivial, simply f (ii) αβ =S(S + 1) Φ α |Φ β . Using the fact that the two states are both spin singlet, one can further prove important symmetries of f (i j) for i = j. For a site pair in a tetrahedron unit, let call the remaining two sites its conjugate site pair; e.g., for the site pair 1-2, its conjugate pair is 3-4. For the site pair i-j, let us define f and then F (i j) is identical to the value for its conjugate site pair: and the sum of these vanish These F's are a square matrix with dimension 2S + 1. These relations (10) and (11) will be referred to in the following as conjugate-pair equivalence and neutrality identity, respectively. 2 The spin gap is ∆ s = J in the decoupling limit J = 0. At small J /J, the ground-state energy has a correction starting from the order J 2 /J, while the S unit = 1 state can hop from one tetrahedron to neighboring ones with matrix element proportional to J . Therefore, the leading correction in the spin gap is the order J , and ∆ s = J − a(S)J + · · · . 3 The result that the matrix elements are diagonal in spin space does not depend on the choice of basis states. For example, one can use those in Eq. (3) or bases of the irreps of T d group. Before going to perturbation calculations, I quickly prove Eqs. (10) and (11). For the conjugatepair equivalence, it is sufficient to prove f (12)

5/28
αβ , and the following proof does not depend on the choice of site pair i-j. Let S unit be the total spin in the tetrahedron unit, S unit =∑ 4 i=1 S i , and here I drop the label of the unit position r, since all the calculations are limited in one unit. For any tetrahedron singlet state Φ α , the most important relation is Another relation to use is the identity S i · S j = 1 2 S 2 i j − S(S + 1), where S i j ≡ S i + S j is the composite spin of the pair. Then, the relation to prove is equivalent to Using S unit , the 3-4 site pair can be represented by quantities related to the 1-2 site pair This completes the proof of the conjugate-pair equivalence.
Using the conjugate-pair equivalence for three pairs, one can rewrite the neutrality identity as follows It is straightforward to prove this, since the left-hand side is nothing but Φ α | 1 4 S 2 unit |Φ β =0. Thus, the neutrality identity is also proved.
We are now ready to start a degenerate perturbation for constructing an effective Hamiltonian. In perturbation in J , first-order terms vanish and second-order terms only yield a constant energy shift for all states where 3 2 N is the total number of long bonds. Therefore, the leading terms that lift the degeneracy are third order ones, which was explicitly derived for the S = 1 2 case [11]. Beware that there are two types of third-order terms. One corresponds to perturbation paths on triangular loops made of weak bonds alone, and they contribute only a constant energy shift again, where N should be read as the number of the triangular loops, which is identical to the number of original spins. Fig. 3 One process in the third-order perturbation for the tetrahedron units A-C. This starts from the product state of three singlets Φ A α Φ B β Φ C γ , and each zigzag bond depicts a perturbation J S i (r) · S j (r ) to be operated. Shadowed units are tetrahedra excited to S unit =1, while all the others remain in S unit =0. Different orders of operating three perturbations generate 5 other processes, and all of them have the same contribution as the process shown here, since the three perturbations commute to each other.
The other type is what we need and corresponds to perturbation paths on hexagon loops each of which includes three weak bonds. Two examples of the latter type are shown in Fig. 2. The colored hexagon loop in the left panel includes the three tetrahedron units ABC and the corresponding thirdorder perturbation term is given by where v's are excited states of the system of the units ABC and ∆ ε ABC (v) is their excitation energy measured from the ground state value. Note that the conjugate-pair equivalence has been used for the unit A, f (34) = f (12) . One process is depicted in Fig. 3 and this corresponds to the matrix elements in Eq. (18a). The factor 3! comes from the fact that different orders of three S · S's have all the same contribution. The sums ∑ v 1 ,v 2 are originally taken over all the excited states of the units ABC, but matrix elements are nonvanishing only with those v's in which two units are excited to S unit =1 and one unit remains S unit =0. Therefore, the excitation energy is always ∆ ε ABC (v)=2J for any of those v, and one can move the energy denominators to outside the sum. I should emphasize that this is a special feature of the present model, and this simplifies calculations. After doing this, one can modify the sum ∑ v 1 ,v 2 now with including the ground states. This modification does not change the result, because the additionally included ground states have only zero matrix elements. 4 The modified sum 4 Recall Φ β |S µ j |Φ α =0 for any pair of states in the S unit =0 subspace. This is because S µ j |Φ α belongs to the S unit =1 subspace. for the virtual states is taken over the entire Hilbert space of the three units, and therefore we can safely drop this sum, (∑ excited states + ∑ ground states )|v v| = 1. Thus, Eq. (18b) is obtained, and it is straightforward to rewrite that into the final step (18d). There, the unit index A-C is added in the subscript to show explicitly which unit contributes to each f factor. The neutrality identity implies that only two of the operators F (1 j) 's are independent in each tetrahedron unit, and I choose the following two (14) ).
I have checked that these operators transform as basis of the E-irrep upon operations in the T d point group, and with e m = e(θ =2mπ/3) where e(θ ) = (cos θ , sin θ ). As will be shown in Eq. (28) and Appendix B, all of the operators {F 1 j } have eigenvalues The full effective Hamiltonian in the third order is obtained by repeating the same calculation for all the tetrahedron triads participating to the shortest hexagon loops in the breathing pyrochlore lattice. With these new operators τ τ τ, the effective Hamiltonian for the spin-S Heisenberg model (1b) is represented as Here the sum is taken over all the tetrahedron triads explained before. As explained in the caption of Fig. 4, the number of terms in H eff is 12 where N is the number of original spins and the factor 1 3 comes from the fact that each term of this three-unit interaction is counted three times. Corresponding to the choice of three units out of the four sublattices A-D, H eff has 4 types of terms and the prefactors e m 's in Eq. (22) should be chosen as follows depending on sublattices One should note that each e m appears once and only once in each triple product term in Eq. (22).
The prefactor e m is determined by how its tetrahedron unit is connected to the hexagon loop. It is e 0 if the unit is connected with the site pair 1-2 or 3-4, e 1 for 1-3 or 2-4, and e 2 for 1-4 or 2-3. This effective Hamiltonian derived for general S is identical to the one obtained in the previous studies for the S = 1 2 case [11,12]. The only but essential difference is that τ τ τ r now operates in the local singlet space which has the dimension 2S + 1. For S = 1 2 , τ τ τ are a half of Pauli matrices and c 0 = 1 4 , and the Hamiltonian (22) reduces to the effective model in Refs. [11,12] up to the numerical factors. 5 Expanding the triple products in H eff , it is again found that all the terms linear in τ τ τ vanish for any S.
Before proceeding to the next step, I briefly comment on the classical limit S = ∞, and explain that the present perturbative approach fails there. With increasing S to infinity, while the quantum Hamiltonian converges to the classical Heisenberg model, the quantum ground state does not continuously evolve to the ground state of the classical model. The reason is the following.
As explained at the beginning of this section, the use of the effective Hamiltonian is limited to the ground state and the low-energy sector of the original Heisenberg model where the excitation energy is smaller than the spin gap ∆ E < ∆ s = J∆ s ( J J ; S). In the S = ∞ limit, the classical Heisenberg model is defined with classical unit vectors s i (r)'s as H cl = J cl ∑ r,i< j s i (r) · s j (r) +J cl ∑ (i,r),( j,r ) s i (r) · s j (r ). To converge to this upon increasing S in the quantum Hamiltonian (1a), one needs to renormalize spin variables as S i (r) = Ss i (r), and this requires a proper scaling of the exchange constants where J cl and J cl are constants independent of S. This immediately implies that ∆ E < J cl S 2∆s ( J cl J cl ; S) → 0 with S → ∞, and the energy region of the effective model shrinks to zero in the classical limit. Thus, the low-energy region of the quantum Hamiltonian (1a) with finite S is not continuously connected to that in the classical Heisenberg model. In particular, the S = ∞ limit is a singular point for the ground state: the spin rotation symmetry is not broken in the ground state for any finite S, but it is broken in the classical ground state. Therefore, it is impossible to formulate an expansion of 1 S type starting from the classical limit. This contrasts with the case of magnetically ordered states, where the 1 S -expansion correctly describes the ordered ground state and magnon excitations in the corresponding quantum system. The classical Heisenberg model on the pyrochlore lattice is itself exotic, and the ground state is thermodynamically degenerate [6,20].

Spin-pair operators τ τ τ for general S
To analyze the effective Hamiltonian, one needs to know an explicit form of τ τ τ operators, and this is another challenge for S > 1 2 . In this section, I am going to calculate the matrix elements of τ τ τ in terms of the basis states {Φ l } defined in Eq. (3). It turns out useful to introduce uniform and staggered components of spin pair, and their ladder operators, The spin-pair wavefunction φ (12) (l, l z ) is an eigenvector of S 2 12 with the eigenvalue l(l + 1) for any l z . Therefore, this is also the case for our basis functions in the S unit = 0 subspace of tetrahedron unit, and This immediately leads to the matrix elements of τ 1 This matrix is diagonal and traceless, ∑ 2S l=0 (τ 1 ) ll = 0. The largest and smallest eigenvalues are (12) ) is more elaborate, since the operation of S 1 · S 3 or S 1 · S 4 hybridizes different Φ l 's and four-spin nature of the wavefunctions complicates its evaluation. As shown in Appendix B, F (13) is related to F (12) by a unitary transformation, but this needs an involved calculation of many 9 jor 6 j-symbols. I have found a practical way of directly calculating τ 2 for general S, and I explain this in the following.
Matrix element (τ 2 ) l l is given by the overlap integral between Φ l andΨ l ≡ 2(S 1 · S 3 − S 1 · S 4 )Φ l multiplied by factor 1 2 √ 3 . It is important to notice thatΨ l is in the subspace of S unit = 0. This is because the relation [S 2 unit , S 1 · S j ] = 0 leads to the eigenvalue equation Next, I rewriteΨ l to a symmetric form for simplifying further calculation. The definition gives Ψ l = (S 12 + N 12 · N 34 Φ l , and the conjugate-pair equivalence leads to another expressionΨ l = 2( Averaging these two, one obtains a more symmetric formΨ and I am going to calculate this. Now, let us examine more details ofΨ l . It reads in terms of pair wavefunctions as The goal is to express this in terms of our singlet basis functions {Φ l }. I have not been able to find the formula of operating N in the literature, and so I need to derive it. 10/28 I start with the part operated by N z operator. The definition of the pair wavefunction (4) leads to Among various recursion formulas of the CG coefficient, useful is the one that changes the composite angular momentum 6 [19], with and the projection l z does not change here. Since the coefficients on the right-hand side of Eq. (33) do not depend on m, this leads to the same formula for the pair wavefunction Operation of the ladder operators N ± 12 is more difficult to perform. Instead of their direct operation, it is useful to notice the following identity and operate these commutators instead. In this case, one knows all the necessary matrix elements. Operating S ± 12 changes l z by ±1, and N z 12 changes l by ±1. The result is N ± 12 φ (12) (l, l z ) = ± (l ∓ l z − 1)(l ∓ l z ) B S (l) φ (12) (l − 1, l z ± 1) ∓ (l ± l z + 1)(l ± l z + 2) B S (l + 1) φ (12) (l + 1, l z ± 1).
With these results, we go back to Eq. (31) and sum over l z on the right-hand side. This summation contains two types of products concerning the composite spin: φ (12) (l + ∆ l, ·) ⊗ φ (34) (l + ∆ l, ·), and φ (12) (l + ∆ l, ·) ⊗ φ (34) (l − ∆ l, ·), where ∆ l = ±1. Those of the latter type do not belong to the subspace of S unit = 0, since the two composite spins differ. As proved before,Ψ l is a wavefunction in the S unit = 0 subspace, and therefore these cross terms cancel to each other. The products of the former type contribute to the singlet components Φ l−1 and Φ l+1 asΨ This completes the calculation of τ 2 . The matrix elements are given by and this matrix is tridiagonal with zero diagonal elements. 6 Equation (33) used in the present work is a special case of Eq. (C.20) in Ref. [19], but the result in the reference is erroneous. f (x) there should be multiplied by factor 2. 11/28 6. Mean field theory of the effective model 6.

Mean field equation
The final step is the task of solving the effective Hamiltonian H eff . I do this by a mean field approximation at zero temperature. First, I examine in this section this problem for general S, and later in the following sections obtain explicit solutions for the cases of S = 3 2 and 1 and discuss the results in detail. This approach is equivalent to approximating the ground state by a product of local wavefunctions of all the tetrahedron units and those local wavefunctions are to be optimized. I further assume that the spatial pattern has the cubic unit cell with 16 original spins, corresponding to the four tetrahedron units A-D in Fig. 1, and the translation symmetry is not broken further. In this case, a trial product state reads as where R denotes the position of cubic unit cell. Each ψ X is a (2S+1)-dimensional trial wavefunction in the unit X, and {ψ A , · · · , ψ D } are to be determined by energy minimization. Equivalently, one may define the "order parameters" by a set of 4 two-dimensional real vectors { τ τ τ X } D X=A , where τ τ τ X = ψ X |τ τ τ|ψ X , and determine them by minimizing the mean field energy where c 0 = 5 4 in the S = 3 2 case and 1 3 S(S + 1) for general S. e m = (cos 2mπ 3 , sin 2mπ 3 ) as before and P XY denotes the operation that exchanges the units X and Y . Beware that τ τ τ X is related to two-spin correlation in the unit X This relation holds for general S.
Minimizing E MF with respect to ψ X reduces to an eigenvalue problem, H X MF ψ X = ε X ψ X , with the mean field Hamiltonian and the mean field at the unit A is given by Similar results are also obtained for the other units BCD. One needs to determine τ τ τ A , · · · , τ τ τ D self-consistently.

Single-unit problem
To solve the self-consistent equations in the mean field approach, one needs to calculate τ τ τ X for a given mean field field h X , and I now discuss this problem in more detail. Calculations up to this stage are common for general S. However, solutions of the single unit problem have different characters depending on the value of S, and the case of S = 1 2 is exceptional as I will show below. Eigenstates of H X MF are completely determined by the direction ζ of the mean field, h X /|h X | ≡ e(ζ ). Therefore, it is sufficient to define the dimensionless Hamiltonian and consider its ground state 12/28 ψ 0 (ζ ):H Here ε 0 (ζ ) is the dimensionless ground-state energy and this generally depends on the field direction ζ . The only exception is the case of S = 1 2 , where τ 1 and τ 2 are both half Pauli matrix and therefore ε 0 (ζ ) = − 1 2 for any ζ . Regarding the field direction dependence, it is important to notice that as will be proved later, the T d point group of the tetrahedron unit implies the following properties of the ground state energy This shows that the field anisotropy has the Z 3 symmetry, and also means that ε 0 (ζ ) is extreme for the 6 directions where denotes the derivative with respect to ζ . SinceH MF (0) = −τ 1 andH MF (π) = +τ 1 , the extreme values of ε 0 (ζ ) are given by τ 1 's largest and smallest eigenvalues obtained in Eq. (28) Note that the above arguments claim only local extremeness and do not guarantee these are global maximum or minimum in the whole ζ region. However, calculations for S = 3 2 and 1 show that they are the global maximum and minimum, and this suggests this also holds for general S.
The strength of the Z 3 anisotropy is characterized by the ratio of the minimum and maximum values of the ground state energy For S = 1 2 , this parameter reduces to R anis = 1, and the anisotropy completely vanishes. For larger S, the anisotropy monotonically increases with S and approaches 2 in the S = ∞ limit. One should note that this anisotropy parameter also gives the ratio of the maximum and minimum moduli of the order parameter vector R anis = max ζ | τ τ τ(ζ ) |/ min ζ | τ τ τ(ζ ) | = | τ τ τ(0) |/| τ τ τ(π) |. Now, I quickly prove the symmetries (46). The first and second equalities are related to the two symmetry operations P 5 and P 4 in the tetrahedron unit, respectively, introduced in Appendix A. P 5 is a 3-fold rotation about the site 1, while P 4 is a diagonal mirror operation that exchanges the sites 1 and 2. These operations transform τ τ τ operators as follows where t P n = P −1 n and Eq. (B4) is used for the first two relations. For the relations with P 4 , I have used (P 4 ) l l = (−1) 2S+l δ l l as well as Eqs. (28)  to rotation and mirror in the τ τ τ-space P 4 τ τ τ t P 4 = 1 0 0 −1 τ τ τ, P 5 τ τ τ t P 5 = t R 2 3 π τ τ τ, where R(θ ) ≡ cos θ − sin θ sin θ cos θ .
Using these relations, it is straightforward to show the following transformation for the Hamiltonian Since P 5 and P 4 are both orthogonal matrix, this guarantees that the 6 cases of ±ζ , ±ζ + 2π 3 , and ±ζ − 2π 3 have the identical energy spectrum, and the properties (46) are proved. The ground state wavefunctions are related as Let me also discuss the order parameter defined for the ground state, τ τ τ(ζ ) ≡ ψ 0 (ζ )|τ τ τ|ψ 0 (ζ ) . This expectation value is also transformed as shown in Eq. (51) with symmetry operations, and therefore manifests symmetry breaking of the T d point group. Once the ground-state energy is obtained, one can calculate this without using ψ 0 (ζ ). To this end, the Hellmann-Feynman theorem [21] is useful and this yields the result τ τ τ(ζ ) = −ε 0 (ζ )e(ζ ) − ε 0 (ζ )e (ζ ), where again denotes the derivative with respect to ζ . This means that the component parallel to the applied field has amplitude −ε 0 (ζ ). The transverse component has amplitude −ε 0 (ζ ), and thus does not vanish unless the mean field points to any of the symmetric directions ζ = 1 3 π×(integer).
14/28 7. Mean-field ground state of the S = 3 2 case I now begin investigating specific cases and start from the case of S = 3 2 , which is relevant for the compound LiXCr 4 O 8 [15][16][17], and will obtain the mean-field ground state of the effective Hamiltonian.
In this case, the S unit = 0 space has dimension 4 and its basis states belong to the A 1 -, A 2 -, and E-irreps of the T d point group Here, {Φ l } 3 l=0 were defined in Eq. (3).

τ τ τ operators in the S = 3 2 case and solution of a single unit problem
With this basis set it is straightforward to represent τ 1 and τ 2 . The results read One important point is that they have no matrix elements in the subspace spanned by Φ A 1 and Φ A 2 . This is a consequence of the fact that τ τ τ operators transform as bases of the E-irrep, since the product representations A 1 ⊗ E and A 2 ⊗ E contain neither A 1 or A 2 irrep. The result above shows an interesting difference between these two operators. To see this, let us divide the local S unit = 0 space to two subspaces V + and V − that are spanned by {Φ A 1 , Φ Eu } and {Φ A 2 , Φ Ev }, respectively. I should note that V − is the space of wavefunctions which change sign upon exchange of spins 1 and 2 (i.e., permutation P 4 in Appendix A) or equivalently 3 and 4, Φ ∈ V − → P 4 Φ = −Φ, while Φ ∈ V + → P 4 Φ = +Φ. Then, τ 1 has no finite matrix elements between V + and V − , while τ 2 's finite matrix elements are only between them. This is another manifestation of the transformation (50c).
The eigenvalues of τ 1 are − 5 2 , − 3 2 , 1 2 , 7 2 , while τ 2 's eigenvalues are ± 21/4 + The mean-field HamiltonianH MF (ζ ) is now explicitly represented by a 4 × 4 matrix. I have diagonalized it and found that its four eigenenergies are all non-degenerate for any direction ζ as shown in Fig. 5(a). The lowest eigenvalue is where g(ζ ) is a positive parameter given by  Fig. 5(a), ε 0 (ζ ) is minimum for the three field directions ζ = 2mπ/3 (m=0, 1, 2) while maximum for ζ = (2m + 1)π/3. The order parameter 15/28 τ τ τ(ζ ) is calculated from the formula (54) and its trajectory is plotted in Fig. 5(b). This implies that the size of the local order parameter is limited as | τ τ τ | ≤ 7 2 for any state in the S unit = 0 subspace, and − 5 2 ≤ e m · τ τ τ ≤ 7 2 , which is equivalent to − 15 4 = −S(S + 1) ≤ S 1 · S j ≤ 9 4 = S 2 . The neutrality identity (11) imposes a further constraint on the three S 1 · S j 's. Since their sum −S(S+1) agrees with the lower bound of S 1 ·S j , the partial sum of any two correlations is bounded from above where three j's are all different. This manifests no possibility of two ferromagnetic spin pairs in any tetrahedron unit. Almost always, just two pairs should be antiferromagnetic and the remaining one should be ferromagnetic in the ground state. The only exception is the case of a pair of spin-singlet dimers: one S 1 · S j is antiferromagnetic and the other two are zero. As far as | τ τ τ(ζ ) | is minimum at ζ = π, this result holds for general S, because the neutral identity and the τ 1 's lower bound are common for all S. However, beyond the mean-field approximation or at finite temperature, which I do not discuss in the present work, τ τ τ shrinks to a point in the interior region in Fig. 5(b), and it is also possible that all the three spin pairs are antiferromagnetic, if the shrinking is large.

Solution for a triad of tetrahedron units
With this result, let us to find the lowest-energy solution for one tetrahedron triad, e.g. ABC. The corresponding mean field energy is and the mean field is, for example, h A = −e 0 5 4 − e 1 · τ τ τ B 5 4 − e 2 · τ τ τ C . The lowest-energy solution is the one in which e(θ X ) · τ τ τ X =− 5 2 for all the units X's and its energy is E MF (ABC)/(4J eff ) = −(15/4) 3 ≈−52.73. This means that the order parameters are τ τ τ(ζ ) with ζ =π, 5 3 π, and 1 3 π for θ X =0, 2 3 π, and 4 3 π. Thus, the mean field points towards the direction opposite to e(θ X ) at all the units, and the three τ τ τ 's form an equilateral triangle. See Fig. 6(a). This solution has interesting spin correlations. In each tetrahedron unit, only one pair of bonds have a strong antiferromagnetic  Fig. 6 (a) The unique mean-field ground state for the triad of tetrahedron units ABC. (b) One of the mean-field ground states for the four units ABCD. This is also a solution in the bulk. 16/28 correlation and S i · S j X = 0 for the other bonds. It is important to notice that these strong bonds are parts of the hexagon loop that connects the three tetrahedron units in J perturbation. The strong bonds form a Kekulé pattern of antiferromagnetic correlation on the hexagon loop. The same behavior was already discovered in the S = 1 2 case, and this is the origin of stability of the state obtained [11,12,22].

Solution in the bulk
Let us examine the stability of this configuration when the fourth unit D is attached. This also gives a solution in the bulk. Since the effective interaction takes effect for each triad of tetrahedron units, this attachment increases the number of interacting triads from 1 to 4. For the solution obtained in Sec. 7.2, I have found that all the new interactions relating to D vanish as will be shown below. The mean field at the new unit D is and this vanish if τ τ τ X (X = A, B,C) are fixed as before, because each ( 5 4 − e m · τ τ τ X ) is zero in the equation above. Therefore, τ τ τ D is undetermined and may point to any direction with no energy cost. This does not affect the mean fields h A−C , as far as τ τ τ A−C are fixed to the values of the three-unit solution. For example, in Eq. (44) for h A , the parts e 0 · τ τ τ B = e 0 · τ τ τ C =− 1 2 τ 1,A = 5 4 remove the contribution of τ τ τ D . Therefore, this is a self-consistent solution of the four-unit problem, and its total energy of the four units is identical to that of the three-unit solution. Thus the energy per triad increases. This situation happens in the S = 1 2 case and its mean-field ground state in the four units ABCD is continuously degenerate such that the direction of one τ τ τ is arbitrary [11,12].
The situation completely changes in the S= 3 2 case. The solution above is one self-consistent solution, but there exists another solution with a lower energy. This difference comes from the 3-fold anisotropy in the τ τ τ space in the S = 3 2 case, i.e., dependence of the ground state energy on the field direction ε 0 (ζ ). I have numerically solved the mean-field equations for the 4 two-dimensional vectors { τ τ τ X } (X = A, B,C, D), and found that the ground state is unique except 12-fold degeneracy due to the T d symmetry. All the 12 solutions have paired order parameters. One solution has the following ground state at each unit Its order parameters pair up as and this corresponds to spin correlations units S 1 · S 2 S 1 · S 3 S 1 · S 4 A and D −3.1469 −2.8293 2.2261 B and C 9 4 Table 1 Spin correlations between neighboring tetrahedron units in one mean-field ground state for the S = 3 2 case. Each triple product corresponds to one term in E MF , and the four products are all negative. X 0 in the triple product π(X 0 )= S 1 · S 2 X S 1 · S 3 X S 1 · S 4 X denotes the label of missing sublattice X 0 = {A, B,C, D} − {X, X , X }. Let me explain the degeneracy of this mean-field ground state. In the solution above, two τ τ τ 's point along one of the trigonal axes e 0 . The other two τ τ τ 's are slightly tilted from another trigonal axis e 2 , but that tilt is small (about 1 60 π). This solution is degenerate in two ways. First, it is degenerate with respect to change τ 2 → −τ 2 . Secondly, the choice of two tetrahedron units are also arbitrary. One should note however that the trigonal axis to which their τ τ τ 's point depends on which units are chosen. For example, in the case shown above, the units B and C have τ τ τ pointing to e 0 , which is related to 1-2 and 3-4 site pairs as shown in Eq. (42). This corresponds to the lattice structure, in which the B unit is positioned from the nearest C unit along the direction of either 1-2 or 3-4 bond. Since there are 6 ways of choosing two from ABCD, the total degeneracy of the mean-field ground state is 2 × 6=12, and these solutions are related to each other by symmetry operations of the T d point group.
8. Spin correlation in the mean-field ground state in the S = 3 2 case I now use two-spin correlations S i · S j and reexamine symmetry breaking in the mean-field ground state obtained in the previous section. There are two types of these correlations: one type is the correlations on short bonds inside tetrahedron unit and the other is those on long bonds between neighboring units. As will be shown below, spin correlations are finite only between nearestneighbor sites within the perturbative approach used in the present work. However, they manifest spontaneous breaking of the point group symmetry in spin-singlet order. Spin correlations inside tetrahedron unit are already obtained during the calculation of the meanfield ground state, and their values are listed in Eq. (65). Recall that it is sufficient to see f (1 j) = S 1 · S j ( j=2, 3, and 4) because of the conjugate-pair equivalence. In the mean-field ground state for the S = 3 2 case, two pairs have antiferromagnetic correlations and the other one is ferromagnetic in all the units. In two units among the four (B and C in Eq. (65)), the two antiferromagnetic correlations have the identical value −3, while they differ in the other two units. 18/28 8.1. Spin correlations between neighboring tetrahedron units for general S Correlations between different units require a more elaborate calculation. This is because original spin degrees of freedom {S j (r)} are traced out and the effective Hamiltonian has only spin-pair operators {τ τ τ(r)} in each unit. In the effective Hamiltonian approach, correlations of traced-out degrees of freedom are to be calculated from the hybridization of ground-state wave function with excited states, and this generally needs additional careful perturbative calculation. For example, for the large-U limit of the half-filled Hubbard model, Bulaevskii et al. derived an expression of charge density and current in terms of spin operators [23]. In the present case, we can circumvent complicated calculation and obtain result quickly. The technique to use is Hellmann-Feynman theorem [21]. Let us temporarily generalize the original Hamiltonian such that weak couplings on long bonds are all different depending their positions, J → J i j (r, r ), and then its ground state energy is a function of these parameters, E gs [{J i j (r, r )}]. Its derivative with respect to one parameter is the corresponding spin correlation, and we finally set all the parameters to a uniform value J to come back to the original homogeneous Hamiltonian The approximation in the present work replaces the ground state energy E gs [· · · ] by ∆ E Here, the last term is the mean-field ground state energy of H eff and the first two terms are the energy shift in the second-and third-order perturbation Recall that the combination of neighboring units, r, r or r 1 , r 2 , r 3 , automatically fixes the positions of connected sites, i j or i jk. These energy shifts contribute a homogeneous part of spin correlations on long bonds where ρ ≡ J /J and the fact that each long bond is a part of two triangular loops is used for the part of ∆ E 0 . Non-uniform correlations come from the mean-field energy E MF [· · · ]. Since we need a result only for its first order derivative, one can use for E MF [· · · ] the value given from Eq. (41) by replacing J eff with those local values calculated from J i j (r, r )'s × −c 0 + e 0 · τ τ τ r 1 −c 0 + e 1 · τ τ τ r 2 −c 0 + e 2 · τ τ τ r 3 .
This replacement is exact up to the first order in each coupling in {J i j (r, r )}. The positions of connected sites, i 1 , · · · , i 2 , are also uniquely determined by the choice of three units r 1 , r 2 , r 3 . Note that 19/28 each weak coupling J i j (r, r ) appears just in two terms in the sum above. It is helpful to notice that one does not need to consider the contribution of order parameter deformation ∂ τ τ τ r /∂ J i j (r 1 , r 2 ) in Eq. (66). This is because the deformation couples to δ E MF /δ τ τ τ r , and this vanishes since the mean-field solution minimizes E MF . In the present case, all the bonds connecting the same pair of units XX have an identical value of spin correlations. For example, in Fig. 4(a), A1 unit is connected to surrounding four D units, and spin correlations are all the same: . This value is given by summing the contributions of constant energy shift and the mean field energy where π(X)'s are triple products defined in Table 1. The part of E MF contributes two terms π(B) and π(C) corresponding to 2 hexagon loops shown in Fig. 2.
For correlations of other pairs of sublattices, similar results are obtained and the combinations of the three units on the right-hand side are easily read from Eq. (41) One should note that this identity holds generally for any mean-field solution. In the solution obtained in the previous section for the S = 3 2 case, the order parameters pair up as τ τ τ A = τ τ τ D and τ τ τ B = τ τ τ C , which leads to π(A) = π(D) and π(B) = π(C) , but this degeneracy is not necessary for the above identity.
Like the case of the effective Hamiltonian, the results for the spin correlations obtained up to this stage hold for general S.

Results of the S = 3 2 case
I now apply the derived formula to the S = 3 2 case and calculate spin correlations on long bonds connecting neighboring tetrahedron units. For this case, the constant value is C inter = − 75 16 ρ + 1125 64 ρ 2 ≈ −4.6875ρ + 17.578ρ 2 . Triple products π(X)'s are already calculated and listed in Table 1. Using these results up to ρ 2 , the spin correlations on the long bonds are given as  Finally, I analyze spin structure factor. For that, I use the following definition where δ δ δ i and δ δ δ j are the position of the i-th and j-th spin in tetrahedron unit. Instead of considering structure factor for each i j-pair, I alternatively define this for q in the extended Brillouin zone, not limited to the reduced zone. S(q) has no Bragg peaks and is isotropic in the spin space, because of the absence of magnetic dipole order, and this is a smooth function of q. Naturally, one divides S(q) into two parts which correspond to correlations inside tetrahedron units and those between units: The first part is calculated from spin correlations on short bonds, and the result is given by S intra (q) = 15 4 +f (12) γ xy (q a S ) +f (13) γ zx (q a S ) +f (14) γ yz (q a S ), where the constant term is S(S + 1). Here, the form factor is γ µν (q) = cos(2q µ ) cos(2q ν ) with √ 2a s being the length of short bonds in tetrahedron units.f (i j) = 1 4 ∑ D X=A S i · S j X is the average spin 21/28 correlation between the spin pair i and j inside tetrahedron unit, and this is anisotropic in space depending on the bond direction; f (12) = −0.4485,f (13) = −2.9147,f (14) = −0.3870, and I have used the conjugate-pair equivalence,f 12 =f 34 etc. To see the spatial anisotropy in more detail, I rewrite the q-dependence and separate the part with cubic symmetry where and the amplitudes are given bȳ (12) +f (14) +f (13) The cubic part of the q-dependence is γ cub , and γ u and γ v represent spatially anisotropic components of spin correlations that transforms following E-irrep of the T d point group.
The contribution of correlations between neighboring tetrahedron units is generally given as S inter (q) = 1 4 (I AB + I CD )γ yz (q a L ) + (I AC + I BD )γ zx (q a L ) + (I AD + I BC )γ xy (q a L ) , where √ 2a L is the length of long bonds connecting tetrahedron units. I XY 's values are listed in Eq. (73). The identity (72) guarantees that the three γ's have a common amplitude, and its value is 1 2Ī inter = − 75 32 ρ + 5.2480ρ 2 . Adding up the two parts, one obtains the final result of the spin structure factor where the constant part S(S + 1) = 15 4 is factored out and again ρ = J J .
9. Mean field ground state of the effective model in the S = 1 case Now that we have found that the mean-field ground state in the breathing pyrochlore model differ between the two cases of S = 1 2 and 3 2 , a natural question arises about the intermediate case of S = 1. It turns out that the result is similar to the S = 3 2 case, and I will sketch calculations. In the case of S = 1, the S unit = 0 space at each tetrahedron unit has now dimension 3 and consists of basis states belonging to A 1 -and E-irreps The effective Hamiltonian is the one in Eq. (22) as before, but the constant is c 0 = 2 3 and the two operators are now represented as follows with the local basis set {Φ A 1 , Φ Eu , Φ Ev }, I have used the mean field approximation again for the S = 1 case, and it goes exactly the same as before. A necessary calculation is a solution of the eigenvalue problem for a 3 × 3 matrix of the  Fig. 8 One mean-field ground state of the effective Hamiltonian for S = 1. Trajectory is the ground state value of τ τ τ(ζ ) , and the expectation value ψ|τ τ τ|ψ calculated for any state ψ should be located on or inside this trajectory. Two bounds are τ 1 (ζ = 0) = 5 3 and τ 1 (π) = − 4 3 .
local mean field Hamiltonian (43). I have diagonalized the reduced HamiltonianH MF (ζ ), and found that its lowest eigenvalue is given as At the tetrahedron unit where the mean field h points to the direction e(ζ ), the order parameter is given by Eq. (54) with this new result. Its trajectory upon varying ζ from 0 to 2π is plotted in Fig. 8, and this has a shape of rounded triangle as in the S = 3 2 case. This manifests anisotropy in the order parameter space, but the anisotropy is smaller compared to the S = 3 2 case. As shown in Eq. (49), the anisotropy defined in the order parameter space is R anis = 5 4 , which is smaller than R anis = 7 5 for S = 3 2 . The mean-field ground state for a tetrahedron triad does not depend on the S value and the order parameters at the units ABC have the same configuration as the one shown in Fig. 6(a). A mean-field ground state in the tetrahedron quartet ABCD is also a solution in the bulk. Minimizing the mean field energy (40) with respect to the four order parameters τ τ τ 's, I have searched ground states and found 12 solutions. One of them is These values are plotted in Fig. 8, and this solution has the same nature of the solution in the S = 3 2 case in Fig. 6(b). The spin correlations are as follows The relative difference between these two values is larger compared to the result for the S = 3 2 case shown in Table 1. This is because | τ τ τ A,D | is shorter than | τ τ τ B,C | by about 1.8% here, while the reduction is only 0.5% in the S = 3 2 case. 23/28 I have repeated the same calculation as in Sec. 8.3 and calculated spin structure factor for the S = 1 case S(q) = 2 1 − 1 3 γ cub (q a S ) − 1 3 ρ − 0.4588ρ 2 γ cub (q a L ) + 0.2022γ u (q a S ) + 0.3955γ v (q a S ) , (88) where again ρ = J J . Compared to the S = 3 2 case, the structure factor has smaller amplitudes in its anisotropic parts γ u and γ v . This is related to the fact that the anisotropy R anis in the order parameter is smaller for S = 1.

Summary
In this paper, I have studied the ground state of the antiferromagnetic spin-S Heisenberg spin model on a breathing pyrochlore lattice, and examined a spontaneous breaking of lattice symmetry in the spin-singlet subspace. This lattice has two types of bonds, short and long, and the ratio of the corresponding exchange couplings J /J controls frustration. In the limit of J =0, the ground state is thermodynamically degenerate and the main issue is what is the ground state when 0 < J J and what type of spatial pattern do spin correlations show in the symmetry broken ground state.
Based on the third-order perturbation in J , I have derived an effective Hamiltonian for general S, and examined spin-singlet orders with broken lattice symmetry. It is noticeable that the effective Hamiltonian has a form of three-tetrahedron interactions that is identical to the one of the S= 1 2 case previously studied, and I have shown this with the help of newly found two identities of spin-pair operators. This Hamiltonian is represented in terms of two types of pseudospin operators τ τ τ, and they describe nonuniform correlations of four spins inside the unit. Despite of the identical form of the Hamiltonian, the dimension of τ τ τ's local Hilbert space increases as 2S + 1 with spin. Since their matrix elements were not known, I have used algebras of four-spin composition and calculated these matrix elements for general S.
Using these results, I have analyzed the effective model by a mean-field approximation and investigated its ground state for the special cases of S= 3 2 and 1. I have found that the response of pseudospins has a Z 3 anisotropy in (τ 1 , τ 2 ) space when S > 1 2 , and this has a critical effect on pseudospin orders. In contrast to the previous case of S= 1 2 , the ground state of four tetrahedron units has only a few stable configurations, and they have no further frustration when the configuration is repeated in space. Thus, this is the ground state of the entire system within the mean-field approximation. Actually, the ground state is uniquely determined in the sense that multiple solutions are related to each other by the lattice symmetry.
I have studied in detail spatial pattern of spin correlations in the ground state. Each tetrahedron unit exhibits one of the two types of internal spin correlations as shown in Fig. 7. Among six bonds in each unit, four bonds have either strong or weak antiferromagnetic spin correlations, while the other two bonds have ferromagnetic correlations. As for correlations between different units, they are always antiferromagnetic. I have also calculated the spin structure factor S(q), which may be compared to the energy-integrated value of neutron inelastic scattering. It is found that the amplitude of symmetry broken partsf u ,f v is comparable to the symmetric partf cub .
Quantum fluctuations between different units are neglected in the present work, but two features may justify this approximation. One is the three-dimensionality of the lattice, and quantum fluctuations are much smaller than in lower-dimensional cases. The second reason is the presence of the Z 3 anisotropy. A half of the local order parameters point to one of the directions that minimize the anisotropy energy, and the other half also only slightly tilt from another lowest-energy direction. Therefore, their fluctuations are expected not so large.

24/28
Finally, I make a comment on the implication of the present result for interpreting experiment results in the LiXCr 4 O 8 compounds. It is not fair to compare the theoretical results of this work to experimental work, since the present theory does not take account of two important features in these materials. One is the magnetoelastic effects, i.e., coupling of spins and phonons. This is particularly important because the ground state of the spin system breaks the lattice translation and rotation symmetries. Coupling to the corresponding phonon modes should have a large contribution to the ground state energy, and this may change relative stability among different quantum states. The other is the spin anisotropy due to large S-value. It is legitimate to assume that orbital angular moments of Cr ions are quenched as a starting point, but fluctuations in Cr valency are not zero and this leads to anisotropic corrections to the Heisenberg couplings. I believe that despite these points the results and predictions in the present work are useful for discussing magnetic parts of possible symmetry breaking in the breathing pyrochlore materials.  Φ for S ≤ 4. χ 2 =χ(IC 4 ) and χ 5 =χ(C 3 ). The multiplicity is 0 for the T 1 -and T 2 -irreps.