Breathing chromium spinels: a showcase for a variety of pyrochlore Heisenberg Hamiltonians

We address the long-standing problem of the microscopic origin of the richly diverse phenomena in the chromium breathing pyrochlore material family. Combining electronic structure and renormalization group techniques we resolve the magnetic interactions and analyze their reciprocal-space susceptibility. We show that the physics of these materials is principally governed by long-range Heisenberg Hamiltonian interactions, a hitherto unappreciated fact. Our calculations uncover that in these isostructural compounds, the choice of chalcogen triggers a proximity of the materials to classical spin liquids featuring degenerate manifolds of wave-vectors of different dimensions: A Coulomb phase with three-dimensional degeneracy for LiInCr4O8 and LiGaCr4O8, a spiral spin liquid with two-dimensional degeneracy for CuInCr4Se8 and one-dimensional line degeneracies characteristic of the facecentered cubic antiferromagnet for LiInCr4S8, LiGaCr4S8, and CuInCr4S8. The surprisingly complex array of prototypical pyrochlore behaviors we discovered in chromium spinels may inspire studies of transition paths between different semi-classical spin liquids by doping or pressure.


INTRODUCTION
Over the past 30 years, materials with magnetic moments on the vertices of networks of corner-shared triangular or tetrahedral units have played center stage in the experimental search for exotic magnetic states driven by competing, i.e. frustrated, interactions. 1 In this context, the pyrochlore lattice of cornershared tetrahedra has emerged as a quintessential example of high frustration in three dimensions with a multitude of new phenomena having been uncovered. 2 In particular, for rare-earth pyrochlore oxides with a trivalent 4f magnetic ion, there is a wealth of stoichiometric compounds available in large single-crystal form, necessary for detailed neutron scattering studies. This, along with the synergetic dialog between theory and experiments, has led to the experimental discovery and the theoretical rationalization of numerous phenomena like spin ice physics, magnetic moment fragmentation, order-by-disorder, fragile splayed-ferromagnetism and candidate spin liquid phenomenology. 3,4 In contrast with the above rare-earth systems characterized by strongly anisotropic interactions at the Oð1Þ K scale, there is a paucity of magnetic compounds involving transition metal ions on a pyrochlore network and coupled via a primarily isotropic Heisenberg exchange of a high (Oð10 2 Þ K) energy scale. The availability of such compounds in single-crystal form would empower researchers to explore and possibly discover novel collective behaviors existing over a wider and more experimentally accessible temperature window. Specifically, this would allow exposing new ways in which perturbations, which might have appeared inconsequential from a cursory assessment, ultimately emerge as the dictating forces behind the low-temperature and low-energy scale physics of the material. In that regard, the recent successful synthesis of magnetic pyrochlore fluorides ABM 2 F 7 (A ¼ Na; B ¼ Ca, Sr; M ¼ Ni, Co, Fe, Mn) in large single crystals is indeed an exciting and most welcome development in the field of highly frustrated magnetism. 5 There is, however, a short-coming with the latter compounds whose implications remain to be fully ascertained: The cation disorder on the non-magnetic A and B sites is expected to randomize the M-M superexchange and, most likely, ultimately drive the low-temperature state of these systems to a semi-classical spin glass phase, albeit with some nontrivial spin dynamics. 5 In this context, the disorder-free so-called breathing chromium spinels, which define a fairly broad range of materials (e.g. LiInCr 4 O 8 , LiGaCr 4 O 8 , LiInCr 4 S 8 , LiGaCr 4 S 8 , CuInCr 4 S 8 , and CuInCr 4 Se 8 ), [6][7][8][9][10][11][12][13][14][15][16] constitute a significant opportunity for carrying out the above research program in and below the Oð10 2 Þ K temperature scale.
The introduction of a breathing degree of freedom to the regular pyrochlore lattice, characterized by the alternation of small and large tetrahedra, results in two inequivalent nearest-neighbor exchange couplings, with J for the small and J 0 for the large tetrahedron, and adds a new level of complexity and thus richness to the regular pyrochlore Heisenberg Hamiltonian. 1 In the breathing chromium spinels, of chemical formula AA 0 Cr 4 X 8 (A = Li, Cu; A 0 = Ga, In; X = O, S, Se), the magnetic Cr 3+ spin S ¼ 3=2 resides on a breathing pyrochlore lattice (see Fig. 1). Of particular interest, thanks to their significantly different size, the ions at the A and A 0 sites do not chemically admix (i.e. they order), and the breathing chromium spinels are not subject to the A=B site disorder and consequential exchange randomness of the aforementioned ABM 2 F 7 pyrochlores. While an ordered arrangement of the A and A 0 cations was first reported over half-a-century ago, 17 it is only recently that breathing chromium spinels have become a 1 subject of focused attention. 6 A variety of experimental results have been reported, including a complex sequence of structural and magnetic transitions in LiInCr 4 O 8 and LiGaCr 4 O 8 , 8,9,11,12 as well as a half-magnetization plateau in high magnetic field measurements on LiInCr 4 O 8 . 13 On the theoretical front, the properties of the classical nearestneighbor breathing pyrochlore Hamiltonian have been investigated, 18 and spin-spin correlations in the ground state 19 and the effects of spin-lattice coupling 20 have been explored. Possible implications of the topological aspects of magnetic excitations in the breathing pyrochlore lattice have been pointed out. 21,22 However, to the best of our knowledge, there appears to have been no attempt to (i) provide a microscopic determination of the Cr-Cr exchange interactions and to (ii) investigate using an unbiased theoretical framework the implications of the concurrent thermal and quantum fluctuations operating in these compounds. Such an endeavor is necessary to gain a deeper microscopic understanding of the already available experimental results beyond the simplest interpretation currently provided, as well as identifying new and exciting avenues of investigation in these materials-this is the purpose of the present work. As we shall discuss later on, different signs (ferromagnetic or antiferromagnetic) exchange couplings J and J 0 are realized in the oxides (X ¼ O), sulfide (X ¼ S) and selenides (X ¼ Se) breathing chromium spinels. This provides an uncharted territory for the investigation of the nontrivial role played by long-range interactions beyond J and J 0 in these materials. This leads us to propose that breathing chromium spinels constitute a promising platform to explore the role of strong competing long-range interactions in a semi-classical (S ¼ 3=2) regime. In particular, we predict these to engender new forms of nontrivial correlations in the cooperative paramagnetic regime characterized by different classically degenerate ground state manifolds of different dimension and intricate magnetic orders at low temperatures. Examples include emergent effective Heisenberg antiferromagnetism on the face-centered cubic lattice, 18 a Coulomb phase, 23 and a spiral spin liquid. 24

RESULTS
Breathing pyrochlore Hamiltonians Room temperature structures. We first perform electronic structure calculations for all presently known room temperature structures of chromium breathing spinels. We use the structures of Okamoto et al. 6 for LiInCr 4 O 8 and LiGaCr 4 O 8 , the structures of Okamoto et al. 16 for LiInCr 4 S 8 , LiGaCr 4 S 8 , and CuInCr 4 S 8 and that of Duda et al. 7 for CuInCr 4 Se 8 . The main exchange connectivity of all structures is illustrated in Fig. 1. We follow the experimental literature 6,8,9,11 in referring to the exchange couplings within the small and large tetrahedra as J and J 0 , respectively. As in the isotropic pyrochlore lattice, 25 there are 12 second nearest neighbors J 2 and 12 third nearest neighbors which split into two symmetry inequivalent classes; six of these, named J 3a , have an in-between Cr 3+ ion, and the six others, named J 3b , do not. Figure 2 presents the first main result of our work. The exchange couplings as defined in the Heisenberg Hamiltonian of the form  (Fig. 2d-f), the large tetrahedron is characterized by ferromagnetic exchange J 0 . However, the biggest surprise here is the presence of substantial second-neighbor and third-neighbor couplings. As each Cr 3+ ion has a total of 24 such bonds, compared to six J or J 0 bonds, they can have a substantial influence on the behavior of the materials. Interestingly, the fact that the experimentally determined Curie-Weiss temperature of LiGaCr 4 S 8 is small is thus readily explained by a cancellation of ferromagnetic J and J 0 by antiferromagnetic J 2 and J 3a , J 3b couplings rather than by opposite signs of J and J 0 as hypothesized in ref. 16 . Opposite signs of J and J 0 are only found for CuInCr 4 S 8 . Finally, the selenide CuInCr 4 Se 8 (Fig. 2c) is dominated by ferromagnetic J and J 0 couplings; however, even for this compound, the third neighbor couplings J 3a and J 3b are antiferromagnetic and, most importantly, non-negligible.
Temperature dependence. For most materials, only a room temperature structure is on record. However, ref. 9   We use J, J 0 to label the two nearest-neighbor exchange paths within small (blue) and large (red) tetrahedra, respectively. Note that J 2 connects Cr 3+ ions of small and large tetrahedra so that there is only one kind of such a coupling. J 3a and J 3b are the two symmetry inequivalent third nearest neighbors; the J 3a features an in-between Cr 3+ ion while the J 3b does not.
P. Ghosh et al.
Even though the low-temperature structure of LiInCr 4 O 8 has a reduced breathing anisotropy of r 0 =r ¼ 1:029, compared to r 0 =r ¼ 1:051 at room temperature, the exchange couplings J and J 0 become nearly identical (see Supplementary Fig. 1), with the room temperature exchange anisotropy of J 0 =J ¼ 0:37 having been eliminated. For LiGaCr 4 S 8 , while the low temperature structural anisotropy is still significant at r 0 =r ¼ 1:070, the two ferromagnetic couplings J and J 0 become very similar near the U value determined by the experimental θ CW (see Supplementary Fig. 1). Thus, we find that the nearest-neighbor exchanges J and J 0 remain antiferromagnetic for LiInCr 4 O 8 and ferromagnetic for LiGaCr 4 S 8 , but in both cases they become very similar at low temperatures. At the same time, the further neighbor couplings remain very small for the oxide and substantial for the sulfide. Energy mapping clearly identifies three different types of breathing pyrochlore Hamiltonians: Mostly antiferromagnetic oxides, sulfides with ferromagnetic large tetrahedron and antiferromagnetic longer range ðJ 2 ; J 3a ; J 3b Þ exchange, and the mostly ferromagnetic selenide perturbed by crucial antiferromagnetic longer range interactions. Thus, we organize our pseudofermion functional renormalization group (FRG) analysis of the Hamiltonians into oxide, sulfide, and selenide sections.

Pseudofermion FRG calculations
Oxides. For LiInCr 4 O 8 and LiGaCr 4 O 8 , the J and J 0 couplings forming the tetrahedra are both antiferromagnetic. With only these two couplings, one would expect the typical bow-tie features with the associated pinch points in the reciprocal space susceptibility at sufficiently low temperatures (see ref. 27 ). In this previous study, we showed that the inclusion of even very weak second nearest-neighbor couplings can dramatically alter the susceptibility profile of the nearest-neighbor ground state phase and destroy the pinch points and bow-ties. 28 As stated above, for both breathing pyrochlore oxides, our DFT analysis reveals the presence of additional antiferromagnetic J 2 , J 3a , and J 3b couplings (see Table 1). At the classical level, it is known that when the dominant interactions J and J 0 are large and antiferromagnetic such that the four spins at the vertices of each tetrahedron can be assumed to sum to zero approximately, then a third neighbor coupling J 3a in an ideal pyrochlore lattice is equivalent at temperature T ( minðJ=k B ; J 0 =k B Þ to a J 2 coupling of opposite sign. 28,29 Thus, effectively, both chromium oxides are approximately described by a pyrochlore Hamiltonian with antiferromagnetic J, J 0 and ferromagnetic J 2 : ðJ 2 À J 3a Þ=J ¼ À0:025 for LiInCr 4 O 8 and À0:011 for LiGaCr 4 is the overall nearest-neighbor energy scale, listed in Supplementary Table 1). According to ref. 27 , we expect that the spectral weight moves away from the pinch point and forms two symmetrical maxima, resulting in a hexagonal cluster pattern of intensities in the ½hhl plane. Indeed, using the exchange couplings obtained for the room temperature structure (see Table 1), we employ PFFRG to calculate the magnetic susceptibility profiles which are shown in Fig. 3a, b, d, e. They are found to closely resemble the susceptibility of the isotropic J 1 −J 2 pyrochlore magnet for antiferromagnetic J 1 and ferromagnetic J 2 . 27 It is worth noting that, due to the relatively weaker effective second nearestneighbor coupling in LiGaCr 4 O 8 , the pinch points and bow-ties are better preserved at low temperature for LiGaCr 4 O 8 (see Fig. 3d) in comparison to LiInCr 4 O 8 (Fig. 3a). However, from our PFFRG calculations we find that for both oxides, and at all temperatures down till the RG flow breakdown, the maxima of the susceptibility are always located at the high-symmetry W point in the extended Brillouin zone, i.e., at a q ¼ 2π a ð2; 1; 0Þ-type ordering vector (see Supplementary Note 2), as found in ref. 27   isotropic S ¼ 3=2 Heisenberg antiferromagnet. A direct space schematic illustration for the spin configuration corresponding to this order is shown in Fig. 4a. In direct space, for this spin configuration, the pyrochlore lattice decomposes into two pairs of fcc sublattices. Within each pair, the two fcc sublattices are separated from each other by a nearest-neighbor vector which is perpendicular to the 2π a component of the ordering vector, e.g., for a 2π a ð2; 1; 0Þ-ordering vector, one pair of fcc sublattices is separated by the vector a 4 ð1; 0; 1Þ and the other pair is separated by the vector a 4 ð1; 0; À1Þ. All fcc sublattices show a stripe antiferromagnetic order, i.e., spins are antiferromagnetically aligned in the 2π a -direction and ferromagnetically in the perpendicular directions. Within each pair, the neighboring spins of different sublattices are aligned antiparallel in spin space. The two pairs can be rotated freely with respect to each other. In reciprocal space, this configuration leads to main Bragg peaks at all 2π a ð2; 1; 0Þ-type vectors which share the same 2π a component, e.g., 2π a ð2; 1; 0Þ, 2π a ðÀ2; 1; 0Þ, 2π a ð0; 1; 2Þ, 2π a ð0; 1; À2Þ, 2π a ð2; À1; 0Þ, 2π a ðÀ2; À1; 0Þ, 2π a ð0; À1; 2Þ, and 2π a ð0; À1; À2Þ. In the reciprocal plane perpendicular to the 2π a -direction, the relative orientation of the sublattices results in subdominant Bragg peaks with half the intensity of the dominant ones at the 2π a ð1; 0; 1Þ-type vectors in this plane, i.e., at 2π a ð1; 0; 1Þ, 2π a ðÀ1; 0; 1Þ, 2π a ð1; 0; À1Þ, 2π a ðÀ1; 0; À1Þ. For the classical model of LiGaCr 4 O 8 , this state is the exact ground state given by the Luttinger-Tisza method and confirmed with the iterative minimization procedure. However, we find that this state is strongly dependent on the interplay between the two kinds of third nearest-neighbor couplings. This leads to the finding that for LiInCr 4 O 8 , due to the larger difference between J 3a and J 3b , an incommensurate order on the sublattices is stabilized as the ground state found by iterative minimization. We would like to emphasize that the stabilization of the orders discussed above rests crucially on the presence and the relative magnitude of the third nearest-neighbor couplings. Indeed, classically, the pure J 1 -J 2 -model on a regular J 1 J ¼ J 0 pyrochlore lattice would feature a q ¼ 0 state, completely different from the 2π a ð2; 1; 0Þ-order found here. To make contact with experimental results, it is useful to assess the temperature dependence of the structure factor as calculated within PFFRG. The relation T ¼ ð2π=3ÞSðS þ 1ÞΛ can be used to relate the infrared cut-off Λ employed in the PFFRG framework to temperature T. 27,30 The relation is obtained by comparing the classical limit (S ! 1) of PFFRG, where only the RPA diagrams contribute, i.e. a mean-field description, with the conventional spin mean-field theory formulated in terms of temperature T instead of Λ. The resulting estimated ordering temperatures are given in Supplementary Note 2. As explained in detail in the Appendix A of ref. 27 , in the S ! 1 limit, the absence of higher diagrammatic orders in 1=S in the one-loop PFFRG used here implies that, in particular for the nearest-neighbor pyrochlore antiferromagnetic Heisenberg model, a spurious divergence of the susceptibility at finite temperature occurs which is not expected from what is well established for the classical problem. 23,27 Thus, our estimates of the ordering temperatures in all of the pyrochlore magnets considered here are systematically too high. 31 Our comparisons to an experimental temperature T are shown for a Λ obeying Λ=Λ c ¼ T exp =T exp c . It is expected that multi-loop implementations of PFFRG will lead to significant improvements for the calculated temperature scales; 32,33 however, these advanced schemes are still under development.
We have also performed the PFFRG calculation for the LiInCr 4 O 8 Hamiltonian corresponding to its low-temperature structure (see Table 1, last but one line). The calculated susceptibility profile (see Supplementary Fig. 4) is seen to be virtually indistinguishable from that obtained from the room temperature structure shown in Fig.  3a, b even though the breathing anisotropy of the Hamiltonian is significantly different, having evolved from J 0 =J ¼ 0:37 at room temperature to J 0 =J ffi 1 at T ¼ 20 K. The general observation of the magnetic response being largely independent of the breathing anisotropy can be understood from a classical picture. As discussed in the Supplementary Note 3, the eigenvalues of the interaction matrix J ðkÞ for a J−J 0 -only model exhibit flat bands at the lowest energies which determine the typical bow-tie features of pyrochlore antiferromagnets. Interestingly, these flat bands persist for any J > 0 and J 0 > 0 18 (see Supplementary Fig. 7), rendering the momentum space profile of the magnetic susceptibility independent of J 0 =J. As a consequence, longer range J 2 , J 3a , J 3b have the opportunity to play a crucial role in determining the magnetic behavior of the system. In the present case, however, we find that the effective second-neighbor coupling of LiInCr 4 O 8 at 20 K given by ðJ 2 À J 3a Þ=J ¼ À0:030 is nearly unchanged compared to room temperature (see Supplementary Table 1) thus leading to almost identical magnetic responses.
It is now interesting to compare and verify the accuracy of the Hamiltonian determined for LiInCr 4 O 8 against polarized neutronscattering data from powder samples of Okamoto et al. 10 . It was shown by Benton et al. 18 that within the self-consistent Gaussian approximation (SCGA), the experimental estimate J 0 =J ¼ 0:1, or even a range 0:05 t J 0 =J t 0:15, is consistent with the neutron data. We also calculate the product jFðQÞj 2 SðQÞ of the magnetic form factor jFðQÞj 2 and the powder averaged structure factor SðQÞ using FðQÞ ¼ À 0:3094e À0:0274ð Q 4π Þ 2 þ 0:36804e À17:0355ð Q 4π Þ 2 þ 0:6559e À6:5236ð Q 4π Þ 2 þ 0:2856; (2) Material for all the materials considered here. 34 The SðQÞ is calculated using both the quantum S ¼ 3=2 PFFRG and classical Monte Carlo (CMC) methods. We show in Fig. 5a 18 is also likely rooted in the fact that in this approach the "hard constraint", namely, that the spin vector on each site has magnitude S, is relaxed and only implemented on average. This may tacitly incorporate some effects of quantum fluctuations by allowing for fluctuations of the classical moments. 35 At T ¼ 150 K, the neutron data essentially follows the Q dependence of the magnetic form factor; 10 while the computed scattering intensities for the room temperature and T ¼ 20 K Hamiltonians are consistent, the agreement with the experimental T ¼ 150 K neutron data is neither satisfactory for the PFFRG nor for the CMC results. A similar observation has already Sulfides. We have performed PFFRG calculations for all three chromium sulfides, LiInCr 4 S 8 , LiGaCr 4 S 8 and CuInCr 4 S 8 . As we find very similar susceptibility profiles for the three compounds, we show the one for CuInCr 4 S 8 in Fig. 3g, h and the other two in Supplementary Fig. 4. CuInCr 4 S 8 has a ferromagnetic J 0 which is the strongest among all the interactions present in this material. At the low temperatures that we are considering here, we can therefore assume that the four spins on the large tetrahedra point in approximately the same direction, thereby behaving as effective spin-6 entities. Since the large tetrahedra are arranged in a face centered cubic (fcc) magnetic lattice (see Supplementary  Fig. 5), the system can effectively be mapped onto a spin-6 Heisenberg Hamiltonian on a fcc lattice with renormalized nearest-neighbor couplings J fcc These couplings are all antiferromagnetic with values J fcc 1 ¼ 1:11 K for LiInCr 4 S 8 , J fcc 1 ¼ 0:95 K for LiGaCr 4 S 8 and J fcc 1 ¼ 2:56 K for CuInCr 4 S 8 . Interestingly, the nearest-neighbor antiferromagnetic classical Heisenberg model on the fcc lattice features a subextensive (O½L) ground state degeneracy in the classical limit (S ! 1), with associated wave vectors of the form q ¼ 2π a ð1; δ; 0Þ and symmetry-related q. 36 This degeneracy appears as lines of strong intensity in the magnetic response. However, as further explained in Supplementary Note 5, the susceptibility profile of our pyrochlore model differs from that of the regular fcc antiferromagnet, since it is modulated by a form factor arising from the presence of a tetrahedral basis of our effective fcc lattice. For the structure factors in the ½hk0 plane, this leads to a square ring-type structure formed by lines of strong and almost constant intensity 18 (see Fig. 3h and Supplementary Fig. 4). The point q ¼ 2π a ð1; 0; 0Þ (and symmetry-related q-space positions) is special because it resides at the junctions of the ground state manifold  q ¼ 2π a ð1; δ; 0Þ and q ¼ 2π a ð1; 0; δÞ. Because of this, it is found that collinear ordered states with the ordering wave vector of the type q ¼ 2π a ð1; 0; 0Þ are selected by thermal, 36 as well as quantum 37 order-by-disorder mechanism. As the PFFRG method incorporates the combined effects of thermal and quantum fluctuations, and thus of their selection effect, it is interesting to note that our analysis reveals maxima at q ¼ 2π a ð1; 0; 0Þ type positions for all three sulfides (see Supplementary Note 2). The corresponding direct space spin configuration is illustrated in Fig. 4b. The astonishing property that all sulfides show almost identical susceptibility profiles even though the ratio J=J 0 varies significantly between the three compounds can be understood by considering the eigenvalues of the interaction matrix J ðkÞ for a pyrochlore system with J > 0 and J 0 < 0. As shown in Supplementary Fig. 7, the lowest band exhibits line-like degeneracies, regardless of the size of J=J 0 , leading to a magnetic response which is largely independent of this ratio. 18 Plumier et al. 38,39 studied CuInCr 4 S 8 using unpolarized neutrons. In order to compare our calculations to this experiment, we subtract the room temperature neutron spectra from the low temperature spectra in order to remove the background. In Fig. 5c, d, this data is compared to the calculated structure factor SðQÞ weighted with the magnetic form factor jFðQÞj 2 of Cr 3+ ions, i.e., jFðQÞj 2 SðQÞ. For the 40 K neutron data, the experiment matches well with our numerical data from PFFRG up to an arbitrary scale factor. At 85 K, while the peak positions match very well between theory and the experimental intensity profile, the intensity has some deviation at low Q. This could be partly due to the uncertainty in the mapping of the flow parameter Λ to temperature. Since 85 K is about 2.4T c for CuInCr 4 S 8 , the calculated structure factor is compared at 2.4Λ c . Nevertheless, the comparison of the spin-spin correlation profile to the available experimental information lends strong support to the validity of the Hamiltonian we determine for CuInCr 4 S 8 . We have verified that all subleading couplings J 2 , J 3a and J 3b as given in Table 1  Selenide. The material CuInCr 4 Se 8 is the only selenide breathing chromium spinel with detailed structure reported for which we can perform calculations. At a first glance, the Hamiltonian parameters (Table 1) appear similar to those for LiGaCr 4 S 8 , with a ferromagnetic large tetrahedron. However, the PFFRG result shown in Fig. 3j-l displays a completely different spin susceptibility. In contrast to a repetition of the broken high intensity lines representative of the one-dimensional sub-extensive degeneracy that we discussed in the previous subsection, we find the twodimensional degeneracy of a sphere in reciprocal space. The ferromagnetic large tetrahedra characterized by the J 0 coupling form a large (S ¼ 6) moment fcc lattice with an effective nearestneighbor interaction given by J fcc 1 ¼ ðJ þ 4J 2 þ 2J 3a þ 2J 3b Þ=16 which evaluates to J fcc 1 ¼ À0:55 K. The ferromagnetic nature of J fcc 1 in the selenide is indicative of the crucial difference between the selenide and the sulfides (for which the J fcc 1 are antiferromagnetic): in contrast to the sulfides, nearly negligible J 2 and smaller J 3a and J 3b couplings do not fully compensate the substantial ferromagnetic J anymore. Nonetheless, clearly no simple ferromagnetic order is realized in CuInCr 4 Se 8 . Rather, the two third-neighbor couplings (taken together) substantially perturb the FM large tetrahedra, which leads to the appearance of a two-dimensional degeneracy of classical ground states. At the ordering temperature (Supplementary Note 2), the system develops an incommensurate spiral magnetic order with a helix pitch vector q % 2π a ð0:521; 0; 0Þ (and symmetry-related points ð0; q; 0Þ and ð0; 0; qÞ). The corresponding magnetic susceptibility profile obtained from PFFRG at the ordering temperature is shown in Fig. 3j-k, wherein one observes peaks in SðqÞ at q ≠ 0. The corresponding direct space spin configuration is illustrated in Fig.  4c. It is worth emphasizing that the shift of the spectral weight away from q ¼ 0 is the combined effect of both J 3a and J 3b couplings of a comparable magnitude. 40 Above the ordering temperature, and within the cooperative paramagnetic regime, we observe that the spectral weight is distributed rather uniformly over the spiral surface suggesting the presence of an approximate spiral spin liquid on the pyrochlore lattice stabilized concomitantly by thermal and quantum fluctuations. While an incommensurate long-range ordered spiral phase of ðq; 0; 0Þ type has been observed for the material ZnCr 2 Se 4 with pyrochlore magnetic lattice, 40 it is only very recently that a spiral spin liquid has been observed in the chromium spinel MgCr 2 O 4 . 41 The spin spiral surface found in the cooperative paramagnetic regime in a PFFRG calculation, is also found to be present in the corresponding classical model at T ¼ 0, as revealed by a Luttinger-Tisza analysis. However, the Luttinger-Tisza eigenstates on this spiral surface generically do not represent real normalizable spin configurations, i.e., the spin length on different fcc sublattices are not equal in magnitude. It turns out that, it is only for the ðq; 0; 0Þ-type ordering vectors (with q % 2π a ð0:45; 0; 0Þ) that the minimal energy can be achieved by a configuration of normalized spins. This ground state configuration then breaks the symmetry of the lattice as it is only governed by one of the three ðq; 0; 0Þ-type vectors. Indeed, our CMC calculations give an ordering vector of q ¼ 2π a ð0:40 ± 0:04; 0; 0Þ at the transition temperature. Although the spiral surface is inaccessible at T ¼ 0 due to the violation of the spin length constraint, thanks to thermal and/or quantum fluctuations it can in principle be made accessible. Indeed, for CuInCr 4 Se 8 , our PFFRG analysis shows that thermal and quantum fluctuations taken together are able to restore a well-defined spin spiral surface (see Fig. 3j, k) right above the ordering temperature. In contrast, our CMC simulations, right above the ordering temperature, find a highly nonuniform distribution of spectral weight, and thus an absence of spiral spin liquid, however, as the temperature is increased to $ 1:5 times the ordering temperature, a spiral surface appears (compare Supplementary Fig. 8 with Fig. 3j, k above). These findings lend support to the role played by quantum fluctuations in aiding the stabilization of a spiral spin liquid. It is worth noting that within the ordered phase, this single-q state does not however represent a single planar spiral throughout the entire lattice. Instead, similar to what is found in CuInCr 4 S 8 , the lattice can be divided into two pairs of fcc sublattices. Within each pair, the sublattices are connected by a nearest-neighbor vector perpendicular to the ordering vector, and the planar spiral orders within each pair are of equal pitch and phase. However, the two pairs of fcc sublattices are found to be out of phase with respect to each other by~5 . This offset in phase is a direct and unique consequence of the breathing anisotropy, and is found to vanish in the isotropic pyrochlore lattice.

DISCUSSION
One of the important findings of this study is the significant and heretofore unappreciated role played by long-range exchange interactions on the breathing pyrochlore lattice. While we have previously pointed out the sensitivity of the spin structure factor to the sign and size of the subleading couplings in the J 1 −J 2 quantum pyrochlore Heisenberg antiferromagnet, 27 we have found in the present work for the chromium breathing pyrochlores an impressive demonstration of the decisive role played by longer-range couplings, i.e., the two kinds of symmetry inequivalent third-nearest-neighbor couplings J 3a and J 3b . On the other hand, the effects of breathing anisotropy are shown to be surprisingly minor. Indeed, we observe that in LiInCr 4 O 8 , the effective ferromagnetic second-neighbor coupling of 2.5% of the average nearest-neighbor coupling far outweighs the substantial room temperature breathing anisotropy of J 0 =J ¼ 0:37 in its effect on the structure factor profile. In the case of the sulfides, although the J 0 =J ratios vary between 93 and À1:8 among the three compounds, the susceptibility profiles appear all very similar. In this case, this is rooted in the effective mapping of the system to an emerging effective nearest-neighbor fcc Heisenberg antiferromagnet with the fcc lattice sites being occupied by ferromagnetic tetrahedra featuring a large S ¼ 6 magnetic moment. In this case, we show that the spin structure factor displays line-like degeneracies which are, once more, largely independent of the breathing anisotropy ratio J 0 =J. Finally, for the dominantly ferromagnetic selenide CuInCr 4 Se 8 , we demonstrate that the combined effect of the two third-neighbor couplings J 3a and J 3b drastically perturbs the ground state away from the simple ferromagnetic order. Interestingly, we find that the perturbed state corresponds to an approximate spiral spin liquid above the ordering temperature and within the cooperative paramagnetic regime, where the individual wave vectors form a sphere-like manifold in reciprocal space. Our PFFRG analysis indicates that the combined effect of quantum and thermal fluctuations only leads to a weak order-by-disorder selection into an incommensurate spiral state. Approximate spiral spin liquids on three-dimensional lattices are so far known only on the diamond lattice in MnSc 2 S 4 , 24,31,42 with an occurrence having also been very recently reported for the pyrochlore material MgCr 2 O 4 . 41 We have theoretically investigated six chromium spinels featuring crystalline, and thus magnetic exchange breathing anisotropy. The Hamiltonians which were determined from density functional theory calculations and investigated using the pseudofermion FRG method showcase a colorful selection of magnetic properties arising from frustrated interactions. We find that the oxide compounds are in a perturbed Coulomb phase, 23 reminiscent of the famous spin ice materials. The sulfides display an effective fcc lattice Heisenberg antiferromagnetic-type behavior, and the selenide features incommensurate magnetic correlations. As a unifying feature, all materials are found to be close to a classical degeneracy which shows up as different variants: The oxides are near a phase featuring an extensive number of classical ground states scaling exponentially in the volume of the system. The selenide is approximately degenerate on a two-dimensional surface with a manifold of states scaling exponentially in L 2 (where L is the linear dimension of the system). Finally, the sulfides are characterized by approximate line-like degeneracies with the number of ground states scaling exponentially in the linear dimension L. This variety of different behaviors in a family of related materials represents an attractive feature which promises a wealth of opportunity for future investigations; on the one hand, doping series interpolating between different types of ground states could be very interesting. On the other, it is an intriguing question how other known but only sketchily characterized breathing chromium spinels like CuGaCr 4 S 8 , AgInCr 4 S 8 , CuGaCr 4 Se 8 , AgInCr 4 Se 8 43 fit the picture outlined in this study. Of particular importance are experimental investigations with single crystals which would allow to assess our theoretical predictions.

Energy mapping method
We determine the electronic structure and total energies of the breathing pyrochlores using the full potential local orbital (FPLO) basis set 44 and the generalized gradient approximation (GGA) functional. 45 The GGA+U 46 correction to the exchange and correlation functional is used to deal with the strong electronic correlations on the Cr 3þ 3d orbitals. Two parameters enter these calculations, the on-site interaction strength U and the Hund's rule coupling J H . As the Hund's rule coupling is an intra-atomic interaction, we do not expect it to vary much from one compound to the other, and thus henceforth fix its value to J H ¼ 0:72 eV as commonly employed for Cr 3+ . 47 Meanwhile, the on-site interaction U was fitted from the experimentally determined Curie-Weiss temperatures as explained below. Heisenberg Hamiltonian parameters are extracted using the energy mapping approach: 31,[48][49][50] For this purpose, a 3 1 1 supercell of the primitive unit of the cubic F43m structure with Cm space group is constructed. In this structure, there are 9 independent spins (out of 12 in total), allowing 75 spin configurations with different energies. Total energies are converged using 6 6 6 k points. An arbitrary subset of 19 spin configurations allows us to fit six exchange interactions, using a Heisenberg Hamiltonian of the form of Eq. (1). Note that we count every bond only once. The quality of the fit is seen to be excellent (see Supplementary Fig. 2 for an example), indicating that the method works extremely well for the chromium breathing pyrochlores, may they be dominantly antiferromagnetic or ferromagnetic. Now, the U parameter of the GGA+U functional can be fixed: We calculate the Hamiltonian parameters for a small set of U values, which we choose such that the Curie-Weiss temperatures Fig. 1) calculated from these Hamiltonians cover a range of temperatures which includes the experimental θ CW . We find the U value from the condition θ theor CW ¼ θ exp CW by interpolation (see Fig. 2). For the breathing pyrochlores, we find the θ theor CW to be monotonous functions of U and thus to allow for a unique solution. Some more discussion on θ CW and U can be found in Supplementary Note 1.

Pseudofermion FRG method
In the PFFRG scheme, 51 a general spin-S Heisenberg Hamiltonian, e.g. Eq. (1), is treated by first re-expressing the spin operator at each lattice site by 2S spin-1/2 degrees of freedom using Abrikosov pseudofermions. Aside from the desired Hilbert space of spin-S states this approach also introduces unphysical states with lower spin magnitudes. Since these states contribute less to the energy of the system compared to the physical ones, they act as excitations and, thus, naturally get excluded at low temperatures in the RG process. 27 Next, the pseudo-fermionic Hamiltonian is treated using the FRG formalism. [51][52][53][54] Here, the bare fermionic Green's function in Matsubara space is equipped with an infrared frequency cut-off Λ suppressing the fermionic propagation at low frequencies. The PFFRG flow equations then follow from the Λ dependence of the m-particle irreducible vertex functions which are given by an exact and infinite hierarchy of coupled integro-differential flow equations. 52 For a numerical implementation, this hierarchy however needs to be truncated, which is done by neglecting three-particle and higher vertex functions. Still, and most importantly, the truncated version of the method exactly sums up the full set of fermionic Feynman diagrams in two limits taken separately, namely, the large-S limit 53 and the large-N limit. 54 In doing so, it provides an unbiased analysis of the competition between magnetic ordering and disordering tendencies. 27,[51][52][53] The flow equations are numerically solved in direct space and, after the Fourier transform of the spin-spin correlations, one obtains the static susceptibility in reciprocal space SðqÞ as a function of Λ. Note that the RG parameter Λ can be interpreted as the temperature T 27 and is therefore used as a proxy to investigate the thermal evolution of the magnetic response. In the present study, we evaluate spin correlators within a cube of edge length of 5a, where a is the cubic lattice constant, which incorporates a total of 2315 correlated sites, producing well-converged results with a high q-space resolution. Furthermore, we approximate the frequency dependence of the vertex functions by discrete grids containing 64 points for each frequency variable. An ordered system is identified by cusps or kinks in the susceptibility flow signaling a magnetic instability towards either magnetic or valence bond order. 51 In contrast, a smooth flow of the susceptibility down to Λ ! 0 indicates a magnetically disordered state, that is, a putative (quantum) spin liquid. 51

Luttinger-Tisza method
In the Luttinger-Tisza method 55-57 the ground state and energy spectrum of a classical spin system are approximately calculated. The classical limit of Eq. (1) consists of replacing the spin operators by classical vectors, which are normalized (the so called strong constraint). To find an approximate ground state, this constraint is relaxed to be fulfilled only on average over the whole system. This allows to Fourier transform the interaction matrix J ij P. Ghosh et al.
with respect to the underlying Bravais lattice of the actual spin lattice, leading to a matrix J αβ ðkÞ, where α; β label the basis points of the lattice.
The eigenvalues λðkÞ of this matrix subsequently give the energies of spin spiral states with wave vector k, the aforementioned energy spectrum, and the absolute value of the corresponding eigenvector components give the relative length of the vector spins on the different basis points. This means, that for such a state to fulfill the strong constraint, all components have to have equal absolute value.
If the lattice is of Bravais type, there is only one basis point and therefore the Fourier transformed interaction matrix is a scalar, which means this condition is trivially fulfilled, rendering the Luttinger-Tisza method exact on such lattices. In this case, we can also, through the equivalence of the Luttinger-Tisza method and PFFRG for S ! 1, 53 calculate the classical spin susceptibility χðkÞ $ 1 2 3T À 1 JðkÞ : Iterative minimization method By employing the iterative minimization method, 27,57 we find the ground state of the classical version of Eq. (1) obtained by replacing the spin operators by normalized vectors. In this scheme, we typically start from a random spin configuration, randomly select one vector-spin S i at a time and rotate it to make it antiparallel to its local field thus minimizing the energy. In a given sweep, we update all the spins contained in the configuration so that on average each spin is updated once. For our calculations, we use 32 cubic unit cells of the pyrochlore lattice in each direction with periodic boundary conditions, totaling to 16 32 3 ¼ 524; 288 spins in the system. Convergence is reached once the total energy difference per sweep is <10 À10 K. To reduce the likelihood of the system getting stuck at a local energy minimum in configuration space, we carry out the minimization for at least 20 different random initial configurations. We also use specific non-random configurations, e.g. spin spirals, as a starting point, to test exact spin parameterizations and, possibly, find lower minimal energies than those achieved by starting with arbitrary random spin configurations.

CMC method
We have performed CMC simulations for the model Hamiltonian parameters of the compounds as obtained from DFT given in Table 1, employing the standard single spin-flip technique. For each Monte Carlo update, we perform the Metropolis moves as many times as the number of spins in the configuration such that, on average, each spin is updated once. Starting from a random spin configuration and after allowing for the system to reach thermal equilibrium, we evaluate the Fourier transform of the equal-time spin-spin correlator, i.e., the static structure factor SðqÞ for 50 different configurations where each configuration is separated by 100 Monte Carlo updates. At lower temperatures, in order to increase the number of accepted moves, we restrict the newly generated spins to a small cone around its predecessor. The simulations are done for a system of 32 32 32 cubic unit cells of the pyrochlore lattice with periodic boundary conditions in each direction, amounting to 16 32 3 spins. The entire simulations are independently carried out 15 times and the average of the results is taken for suppression of numerical noise.

DATA AVAILABILITY
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.