Signatures of a magnetic superstructure phase induced by ultrahigh magnetic fields in a breathing pyrochlore antiferromagnet

Significance The search for exotic magnetic states, such as a quantum spin liquid and a magnetic superstructure, in geometrically frustrated magnets has been a central research topic in the recent 30 y. Theoretically, the mutual coupling of spin and lattice degrees of freedom has been proposed to induce rich magnetic phases in various spin models. Here, we observe an unconventional multistep magnetostructural transition in a breathing pyrochlore antiferromagnet in ultrahigh magnetic fields, signaling the emergence of a magnetic superstructure phase with a periodic array of collinear 3-up-1-down and canted antiferromagnetic states. This finding could be attributed to the interplay between the spin–lattice coupling and breathing anisotropy.

The mutual coupling of spin and lattice degrees of freedom is ubiquitous in magnetic materials and potentially creates exotic magnetic states in response to the external magnetic field. Particularly, geometrically frustrated magnets serve as a fertile playground for realizing magnetic superstructure phases. Here, we observe an unconventional two-step magnetostructural transition prior to a half-magnetization plateau in a breathing pyrochlore chromium spinel by means of state-of-the-art magnetization and magnetostriction measurements in ultrahigh magnetic fields available up to 600 T. Considering a microscopic magnetoelastic theory, the intermediate-field phase can be assigned to a magnetic superstructure with a three-dimensional periodic array of 3-up-1-down and canted 2-up-2-down spin molecules. We attribute the emergence of the magnetic superstructure to a unique combination of the strong spin-lattice coupling and large breathing anisotropy. frustrated magnetism | breathing anisotropy | spin-lattice coupling | ultrahigh magnetic fields Superstructures in crystalline solids, where the unit cell of atomic arrangements or electronic states is an integer multiple of the original primitive cell of the lattice, have received considerable attention because of their complexity as well as their association with exotic physical properties. Prominent examples are surface superstructures such as the Si(111)-7 × 7 state (1), charge/orbital ordering in perovskite manganites (2), and charge density wave states in van der Waals transition-metal dichalcogenides (3) and topological kagome superconductors (4). For frustrated spin systems, a variety of quantum-entangled magnetic superstructures can appear in the external magnetic field. Of particular interest are a series of magnon crystals in a spin-1/2 kagome Heisenberg antiferromagnet (5) and successive transformations of singlet-triplet superstructures in a spin-1/2 orthogonal-dimer Heisenberg antiferromagnet (6). These phenomena are accompanied by a cascade of fractional magnetization plateaus (7)(8)(9)(10)(11). When a spin system interacts with the lattice degrees of freedom, the structural instability can facilitate the formation of a magnetic long-range order, potentially yielding a spin-lattice-coupled superstructure. The concept of a spin-lattice coupling was first theoretically proposed as the spin-Peierls transition in a spin-1/2 Heisenberg antiferromagnetic (AFM) chain (12,13), which was then demonstrated in a number of quasi-one-dimensional compounds (14,15). In this mechanism, the spin dimerization is induced by individual atomic displacements with a doubling of the unit cell, where the system gains exchange energy at the cost of elastic energy. A similar phase transition was subsequently found in three-dimensional semiclassical spin systems represented by chromium-based spinels ACr 2 O 4 (A = Zn, Cd, and Hg) (16)(17)(18)(19). At low temperatures, the strong geometrical frustration inherent in the corner-sharing tetrahedral network, i.e., the pyrochlore lattice, of spin-3/2 Cr 3+ ions is relieved by the local tetrahedral distortion (called the spin Jahn-Teller effect), resulting in a 2-up-2-down magnetic long-range order (16,19) (Fig. 1 A and B). Interestingly, the application of a magnetic field induces another spin-lattice-coupled long-range order, a 3-up-1-down state (19,20) (Fig. 1C). This phase transition is accompanied by a steep magnetization jump (21)(22)(23) and giant magnetostriction (24) followed by a robust half-magnetization plateau, which can be explained in the framework of microscopic magnetoelastic theories incorporating local phonon modes (25,26). These findings have triggered further studies on the spin-lattice-coupling physics for various spin models (27)(28)(29)(30)(31)(32)(33).
Recently, the introduction of a breathing anisotropy, i.e., a spatial modulation of magnetic interactions, has been recognized as a new approach to control the ground state of various spin models (34)(35)(36). In this context, the emergence of spin-lattice-coupled superstructures in addition to the conventional 2-up-2-down and 3-up-1-down states (not superstructures)

Significance
The search for exotic magnetic states, such as a quantum spin liquid and a magnetic superstructure, in geometrically frustrated magnets has been a central research topic in the recent 30 y. Theoretically, the mutual coupling of spin and lattice degrees of freedom has been proposed to induce rich magnetic phases in various spin models. Here, we observe an unconventional multistep magnetostructural transition in a breathing pyrochlore antiferromagnet in ultrahigh magnetic fields, signaling the emergence of a magnetic superstructure phase with a periodic array of collinear 3-up-1-down and canted antiferromagnetic states. This finding could be attributed to the interplay between the spinlattice coupling and breathing anisotropy.
is theoretically proposed in the breathing pyrochlore lattice where the neighboring tetrahedra differ in size in an alternating pattern ( Fig. 2A) (32). Here, we report the possible realization of this theoretical prediction in a model compound of the breathing pyrochlore antiferromagnet, LiGaCr 4 O 8 (36)(37)(38)(39)(40). By means of state-of-the-art magnetization and magnetostriction measurement techniques available in ultrahigh magnetic fields of up to 600 T, we demonstrate that LiGaCr 4 O 8 exhibits a two-step magnetostructural transition in an intermediate field range between 150 T and 200 T, followed by a half-magnetization plateau up to ∼420 T. The effective spin Hamiltonian for LiGaCr 4 O 8 is established based on densityfunctional-theory (DFT) calculations and classical Monte Carlo (MC) simulations. We show that the combination of the strong spin-lattice coupling and large breathing anisotropy can be responsible for stabilizing the intermediate-field phase, which we assign to a tetrahedron-based superstructure with a three-dimensional periodic array of 3-up-1-down and canted 2-up-2-down spin molecules.

Results
Breathing Pyrochlore Antiferromagnet LiGaCr 4 O 8 . We first introduce the basic properties of the present target LiGaCr 4 O 8 (36). The crystallographic ordering of nonmagnetic cations Li + and Ga 3+ identical to the zinc-blende structure (Fig. 2B) leads to the breathing pyrochlore Cr network, as shown in Fig. 2A. The resultant two inequivalent nearest-neighbor Cr-Cr bonds, whose lengths are r = 2.970 Å and r′ = 2.867 Å at room temperature (36), are characterized by two distinct AFM exchange couplings J and J′: [1] where S i and S j denote the classical Heisenberg spins, and 〉i, j〈 A and 〉i, j〈 B stand for nearest-neighbor sites within tetrahedra A and B, respectively ( Fig. 2A). Note that J is equal to J′ for the conventional Cr spinel ACr 2 O 4 . For LiGaCr 4 O 8 , J′ was initially believed to be stronger than J (36) judging from the empirical relationship between the strength of the AFM exchange coupling and the lattice constant (equivalently, the Cr-Cr bond length) in ACr 2 O 4 (42), though the actual situation does not seem so simple. From the structural point of view, the inversion symmetry breaking at the local Cr site would give rise to the anisotropic deformation of Cr 3d 3 orbitals. Indeed, previous DFT energy mapping revealed the dominance of J at room temperature: J/k B = 100 K and J′/k B = 66.2 K (43), k B being the Boltzmann constant. Furthermore, the breathing anisotropy J′/J can be strongly temperature dependent as exemplified in a related compound LiInCr 4 O 8 (43).
In order to evaluate J′/J for LiGaCr 4 O 8 at low temperatures, we measured powder X-ray diffraction patterns at 20 K and performed the Rietveld analysis (SI Appendix, Note 1) and DFT energy mapping (Materials and Methods). Since the coexistence of tetragonal and cubic phases are reported below the ordering temperature T N = 14 K (38) (SI Appendix, Note 2), we can expect reliable information on the exchange parameters from the structural refinement at a temperature slightly higher than T N . Notably, tetrahedra A and B are found to be reversed in size (r = 2.907 Å and r′ = 2.919 Å) while the cubic F-43m structure is preserved at 20 K. Fig. 2C shows the DFT energy mapping up to the third nearest-neighbor exchange couplings ( Fig. 2A). Here, the experimental Weiss temperature Θ W = -659 K (36) is used to fix the on-site interaction strength U. We find a significant reduction of J′/J on cooling: J′/J = 0.66 at room temperature (43) whereas J′/J = 0.04 at 20 K (J/k B = 166.6 K and J′/k B = 6.9 K), indicating that LiGaCr 4 O 8 is located close to the limit of decoupled AFM tetrahedra at low temperatures. The estimation of J′/J is nearly independent of U (J′/J = 0.04 ∼ 0.10 for U = 1.25 ∼ 2.75 eV), and the second and third nearest-neighbor exchange couplings are negligibly weak.

Physical Property Measurements in Ultrahigh Magnetic Fields.
As suggested from the large negative Weiss temperature Θ W = -659 K (36), ultrahigh magnetic fields of several hundreds of tesla are required to exceed the strong AFM exchange interactions in LiGaCr 4 O 8 . Accordingly, we employed two kinds of destructivetype pulsed megagauss generators: an electromagnetic flux compression (EMFC) system (44) (Fig. 3A) and a single-turn coil (STC) system (45) (Fig. 3B). The former allows physical property measurements up to ∼600 T only in the field-increasing process, whereas the latter offers the data up to ∼200 T both in the field-increasing and -decreasing processes. Typical magnetic-field waveforms generated by the EMFC and STC systems are shown by gray lines in Fig. 3 C and D, respectively. Physical property measurements on magnetic insulators in the EMFC system were so far limited to the Faraday rotation and optical absorption measurements (23). These magneto-optical detection schemes, however, could not be applied to the present polycrystalline LiGaCr 4 O 8 samples. In this work, we have extended two existing techniques implemented in the STC system to the EMFC system to observe successive phase transitions of LiGaCr 4 O 8 up to ∼600 T: i) a magnetization measurement using the induction method (31,46,47) and ii) a magnetostriction measurement using the optical fiber-Bragg-grating (FBG) method (48,49) (SI Appendix, Notes 3 and 6). Fig. 3E summarizes the magnetization data of LiGaCr 4 O 8 measured at ∼5 K. In the STC system, a linear increase in the magnetization M with respect to the external magnetic field B is observed up to a generated maximum field B max of 145 T, where M reaches 0.70 μ B /Cr, suggesting that spins are smoothly canting from the 2-up-2-down AFM ground state (38). Here, the absolute value of M is calibrated using the magnetization data obtained in a SQUID magnetometer MPMS up to 7 T (SI Appendix, Note 2) and in a nondestructive pulsed magnet up to 51 T (Fig. 3E). Upon the application of a higher field using the EMFC system, dM/dt anomalies with a small hump and a subsequent large hump are observed at B c1 = 159 T and B c2 = 171 T, respectively (Inset of Fig. 3C), indicating a two-step metamagnetic transition. We ensure the reproducibility in three independent experiments with different setups (SI Appendix, Note 4). By subtracting the background component from the observed dM/dt profile (SI Appendix, Note 5), M as a function of B and its field derivative dM/dB between 150 T and 200 T are obtained as shown by red and pink lines, respectively, in Fig. 3E. The second magnetization jump, ΔM 2 , is approximately six times larger than the first one, ΔM 1 , if the B-linear component is subtracted (Fig. 3E). The absolute value of M at 200 T is 1.35 ± 0.07 μ B /Cr, which is close to half the expected saturation magnetization of ∼3 μ B /Cr given that the g-value is estimated to be g = 1.98 ∼ 2.08 (36,37). The present observation is most likely ascribed to the appearance of a half-magnetization plateau as observed in (21-23), given that the magnetization jump is blunted by thermal fluctuations or crystallographic disorders. The existence of an intermediate-field phase between B c1 and B c2 is supported by the magnetostriction measurements. Fig. 3F shows the magnetostriction data of LiGaCr 4 O 8 measured at 5 K using the STC system with B max = 167 T (Data #1) and 185 T (Data #2). The temporal evolution of the relative sample-length change ∆L/L for Data #2 is shown in Fig. 3D. Data #1 and #2 were taken in the same setup and are in excellent agreement with each other. As seen in Data #1, a rapid lattice expansion starts at around 160 T, which corresponds to the phase transition at B c1 . Furthermore, as seen in Data #2, the increase in ∆L/L is accelerated above B c2 . These observations indicate that both transitions at B c1 and B c2 are likely first order accompanied by a structural phase transition.

Two-Step Magnetostructural Transition at around 200 T.
Field Width of Half-Magnetization Plateau. Theoretically, the magnetic-field width of a half-magnetization plateau correlates with the strength of the spin-lattice coupling in the (breathing) pyrochlore-lattice Heisenberg antiferromagnet (25,32). Therefore, it is important to determine the termination field of the halfmagnetization plateau B c3 for further discussions, although our magnetization measurement could not detect additional phase transitions above 200 T due to poor sensitivity of our probe in the high-field region (SI Appendix, Note 4).
As an alternative way to determine B c3 , we measured the magnetostriction of LiGaCr 4 O 8 at 5 K in the EMFC system. Considering the relation ∆L/L ∝ M 2 derived from the magnetoelastic theory (24), there should be no lattice-constant change in the half-magnetization plateau region, followed by a significant volume expansion on the high-field side. Fig. 4 shows three ∆L/L-B curves (Data #3 ∼ #5) obtained in independent sample settings. The raw data are shown in SI Appendix, Note 7, and the analytical procedure is described in SI Appendix, Note 8. For Data #3, a plateau-like behavior is observed from 200 T up to B max = 360 T. For Data #4 and #5 with B max ≈ 600 T, on the other hand, an upturn behavior is clearly observed at around 420 T, signaling the occurrence of a phase transition from the 3-up-1-down to a higher-field spin-canted phase. We hence determine B c3 ≈ 420 T. Furthermore, another ∆L/L kink is observed at around 550 T, potentially reflecting a phase transition to a paramagnetic phase.
Classical Monte Carlo Simulations. To get more insights into the field-induced phase transitions of LiGaCr 4 O 8 , we employ a microscopic magnetoelastic model on the breathing pyrochlore lattice (32) incorporating the Einstein site phonons assuming the independent displacement of each site (26). The effective spin Hamiltonian H eff is [2] where  0 contains minimal exchange terms expressed as Eq. 1, h is the external magnetic field applied along the z-axis, and [3] where b (b′) is a dimensionless parameter representing the strength of the spin-lattice coupling between the nearest-neighbor sites within tetrahedra A (B), N A (i) (N B (i)) denotes a set of the nearestneighbor sites of site i within tetrahedra A (B), and e ij is the unit vector oriented from site i to j. Possible effects from the single-ion anisotropy can be neglected here because of the quenched orbital angular momentum of the Cr 3+ ion, as indicated by the g-value [g = 1.98 ∼ 2.08 (36,37)]. In fact, the magnetization curves of CdCr 2 O 4 [g = 2.06 (17) (18). As for phonon contributions, the displacement energies of atoms are quadratic in the Hamiltonian so that they can be exactly traced out through the standard Gaussian integration (for the derivation process of Eq. 3, see SI Appendix, Note 9). Consequently, the phonon-mediated spin interactions  SLC consist of biquadratic terms favoring collinear spin configurations (first and second terms) and three-body quartic terms responsible for lifting the macroscopic degeneracy (third and last terms). Note that Eq. 2 with J′/J = 1 can reproduce a half-magnetization plateau with a 16-sublattice 3-up-1-down magnetic structure observed in HgCr 2 O 4 (19) and CdCr 2 O 4 (20), where up-up-up-down chains run along all the six equivalent [110] directions (26). The magnetic phase diagrams for the moderately breathing case of J′/J = 0.6 and 0.2 are reported in ref. 32.
As mentioned above, LiGaCr 4 O 8 is characterized by a large breathing anisotropy of J′/J = 0.04 just above T N according to the DFT calculations, although the situation would be more complex below T N (38). Here, we set J′/J = 0.1 as a typical value for the strongly breathing case in Eq. 2 and performed classical MC simulations to calculate the magnetization curves for various sets of the spin-lattice coupling parameters b and b′. Fig. 5 A-C shows the calculation results for a system size of N = 16L 3 spins with L = 12 in three typical cases: a) b = b′ = 0.05, b) b = b′ = 0.10, and c) b = b′ =0.15. Regardless of b and b′, phase I with an 8-sublattice canted 2-up-2-down state (Fig. 5D) and phase II with a 16-sublattice 3-up-1-down state (Fig. 5E) appear followed by a higher-field phase II' (the magnetic structure of phase II' is likely to be more complicated than a simple canted 3-up-1-down state, and its identification is beyond the scope of this work). An important difference manifests just below phase II. In the weak spinlattice coupling case of b = b′ = 0.05, the system undergoes a first-order transition from phase I to III with a 16-sublattice "1-up-1-down+V"-type spin correlation (Fig. 5F), which is continuously connected to phase II. With increasing the spin-lattice coupling, an intermediate-field phase once disappears for b = b′ = 0.10, and then another phase, phase IV, appears for b = b′ = 0.15. Remarkably, the magnetic structure of phase IV is a mixture of the canted 2-up-2-down and 3-up-1-down units of tetrahedra A in the ratio of 1:2, forming a long-range order with a 6 × 6 × 6 magnetic unit cell (Fig. 5 G and H). The detail of the ordering pattern is described in ref. 32. We confirm that, for system sizes in which L is not a multiple of six, this state cannot be obtained, and instead a metastable state with higher energy appears (SI Appendix, Note 10). Note that phase IV appears also for J′/J = 0.2 with larger values of b and b′ (b = b′ ≳ 0.20) but does not appear for J′/J = 0.6 with any values of b and b′ (32).

Discussion
From the above results, we find that the theoretical M-B curve for the strong spin-lattice coupling case of b = b′ = 0.15 (Fig. 5C) is qualitatively compatible with the experimental M-B curve of LiGaCr 4 O 8 (Fig. 3E) in that the two-step magnetization jump as well as the relatively wide half-magnetization plateau are reproduced. We hence believe that the experimentally observed intermediate-field phase prior to the half magnetization plateau would be a magnetic superstructure phase induced by the strong spin-lattice coupling and large breathing anisotropy. However, the magnetization in phase IV amounts to m/m s ≈ 0.4 (Fig. 5C), which is larger than the experimental value at B c2 , m/m s ≈ 0.3 ( Fig. 3E). Moreover, the present site-phonon model tends to underestimate the field width of the half-magnetization plateau (26,32), as demonstrated in ACr 2 O 4 with J′/J = 1 (21,22) and also in the present case (SI Appendix, Note 11). These quantitative discrepancies may arise from the missing incorporation of the macroscopic lattice deformation in the site-phonon model. Another microscopic magnetoelastic model assuming independent bond-length changes, i.e., the bond-phonon model, can resolve this issue (SI Appendix, Note 11), although it can describe neither magnetic long-range orders nor additional complicated magnetic transitions found in the site-phonon model (25,32). The development of an extended magnetoelastic Hamiltonian is necessary to fill the gap between experiment and theory. Finally, it is worth mentioning the relationship of the present findings to the magnetism of isolated tetrahedral clusters, namely J′ = 0. If we take into account quantum spins and neglect the spin-lattice coupling, the application of a magnetic field induces the quantization of a total spin number per cluster, which can be regarded as a spin crossover rather than a phase transition. As a consequence, the M-B curve exhibits fractional magnetization plateaus at m/m s = p/2n for spin-n/2 systems (p, n: integer values satisfying 0 < p < 2n). Such a magnetization behavior was indeed observed in a spin-1/2 breathing pyrochlore antiferromagnet Ba 3 Yb 2 Zn 5 O 11 (50) and a spin-3/2 tetrahedral-cluster compound Co 4 B 6 O 13 (51). In LiGaCr 4 O 8 , however, no plateau-like structure is observed at around m/m s = 1/6, i.e., M ≈ 0.5 μ B /Cr (Fig. 3E), suggesting that a spin crossover is masked by the intertetrahedral exchange coupling J′. Also, it is obvious that if we set J′ = 0 and b′ = 0 in the classical Heisenberg model with the spin-lattice coupling Eq. 2, no magnetic superstructure phase appears due to the absence of spin interactions within tetrahedra B. Three-body spin interactions across tetrahedra A and B ( Fig. 2A) originating from the site-dependent local phonons, corresponding to the last term in Eq. 3, play a crucial role in bringing about rich field-induced phases in the breathing pyrochlore system. In summary, we observe a two-step magnetostructural transition between 150 T and 200 T prior to a robust half-magnetization plateau in the breathing pyrochlore antiferromagnet LiGaCr 4 O 8 . Considering the magnetoelastic theory incorporating local phonon modes, the intermediate-field phase can be assigned to a spin-lattice-coupled superstructure with a three-dimensional periodic array of 3-up-1-down and canted 2-up-2-down spin molecules, which we attribute to the strong spin-lattice coupling and large breathing anisotropy. The present work, combining the exotic experimental observations with the microscopic magnetoelastic theory in the complicated three-dimensional frustrated magnet, paves the way for further verifications of intriguing physical phenomena originating from the spin-lattice coupling and/or breathing anisotropy, both of which can be relevant in magnetic materials regardless of the geometry of the underlying crystalline lattice. For instance, the formation of a 3 × 3 spin-lattice-coupled superstructure phase associated with a 1/9-magnetization plateau as well as a √ 3 × √ 3 one associated with a 1/3-magnetization plateau is theoretically predicted in a kagome-lattice Heisenberg antiferromagnet (33); its experimental verification as well as theoretical studies on the effect of introducing the breathing anisotropy remain open questions. In addition, we demonstrate that the electromagnetic induction and FBG methods are powerful tools to detect magnetic and structural phase transitions, respectively, in ultrahigh magnetic fields well above 100 T. These techniques would be applicable to a broad range of materials such as frustrated magnets, spin-crossover systems, heavy-fermion compounds, and superconductors, leading to further flourishing of high magnetic field science.

Materials and Methods
Sample Preparation. Polycrystalline samples of LiGaCr 4 O 8 were obtained from the same batch used in ref. 36, prepared by the conventional solid-state reaction method. The powder samples were formed into rod shape with ∼0.8 mm diameter and ∼1.5 mm length using epoxy (Stycast 1266) for the magnetostriction measurement.
Powder X-ray Diffraction Measurement and Structural Analysis. We performed the powder X-ray diffraction measurement on LiGaCr 4 O 8 at 20 K using a commercial X-ray diffractometer (SmartLab, Rigaku). The incident X-ray beam was monochromated by a Johansson-type monochromator with a Ge(111) crystal to select only Cu-Kα1 radiation. Only a tiny amount of impurity phases were found, ensuring that the sample is of high quality. The Rietveld analysis was performed using the JANA2006 program (52), confirming the cubic F-43m structure. Detailed results are shown in SI Appendix, Note 1.
DFT Calculations. Exchange parameters of LiGaCr 4 O 8 were estimated by the DFT-based energy mapping (53,54). We performed all electron DFT calculations using the full potential local orbital code (55). Note that this technique has proven to be very reliable for the breathing pyrochlore chromium spinels (31,43). We used the generalized gradient approximation (GGA) exchange and correlation functional (56). For the electronic structure calculations, we used the T = 20 K crystal structure with the F-43m space group determined in this work. We accounted for strong correlations on the Cr 3d orbitals by applying a GGA+U exchange correlation functional (57) for several different values of U and J H = 0.72 eV fixed following ref. 58.
Magnetization and Magnetostriction Measurements. The magnetization up to 7 T was measured using a SQUID magnetometer (MPMS; Quantum Design). The magnetization up to 51 T, 145 T, and 200 T was measured by the induction method in a nondestructive pulsed magnet, a horizontal STC system (45), and an EMFC system (44), respectively. The details of experimental setup, data, and analysis method for the magnetization measurements using the EMFC system are presented in SI Appendix, Notes 3-5, which include additional references (59)(60)(61)(62)(63). The longitudinal magnetostriction up to 185 T and ∼600 T was measured by the optical FBG method in the horizontal STC and EMFC systems, respectively. Here, a relative sample-length change ∆L/L was detected by the optical filter method with the resolution of 10 -5 to 10 -4 (48,49). The details of experimental setup, data, and analysis method for the magnetostriction measurements using the EMFC system are presented in SI Appendix, Notes 6-8. All of the experiments were performed at the Institute for Solid State Physics, University of Tokyo, Kashiwa, Japan.
Classical Monte Carlo Simulations. To identify various magnetic phases appearing in the spin Hamiltonian Eq. 2, we performed classical Monte Carlo (MC) simulations in which a spin vector at each lattice site is updated in conventional random and successive over-relaxation-like processes. In the former (latter) process, we tried to rotate a spin in a randomly proposed direction (by the angle π around the local mean field) by using the Metropolis algorithm. In our simulation, 2 × 10 6 MC site sweeps were carried out at each temperature and magnetic field under the periodic boundary condition, and the first half was discarded for thermalization. Observations were done in every 10 MC sweeps and the statistical average was taken over four independent runs. We often encountered various metastable states due to very complicated interactions in the spin-lattice coupling term Eq. 3, and low-temperature spin states obtained in the field-cooling, fieldincreasing, and field-decreasing processes were sometimes different. In such a situation, we compare the thermal average values of the energy of these states and regard the lowest-energy state as the equilibrium state. Since our cubic unit cell contains 16 sites, the total number of spins N is related to the linear system size L via N = 16L 3 . We have checked that the results for L = 6 and 12 are essentially the same and consistent with those for smaller L's, e.g., L = 1, where the thermalization is much easier (the spin configurations for phases I, II, and III shown in Fig. 5 can be described even in the small-size system of L = 1), so that only the result for L = 12 is shown in Fig. 5 A-C. Data, Materials, and Software Availability. All study data are included in the article and/or SI Appendix.