{'fraction_non_alphanumeric': 0.038039896962798775, 'fraction_numerical': 0.02054753489486611, 'mean_word_length': 4.339357108587878, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A cascade of magnetic field induced spin transitions in LaCoO3

We present magnetization and magnetostriction studies of the insulating perovskite LaCoO3 in magnetic fields approaching 100 T. In marked contrast with expectations from single-ion models, the data reveal two distinct first-order spin transitions and well-defined magnetization plateaux. The magnetization at the higher plateau is only about half the saturation value expected for spin-1 Co3+ ions. These findings strongly suggest collective behavior induced by strong interactions between different electronic -- and therefore spin -- configurations of Co3+ ions. We propose a model of these interactions that predicts crystalline spin textures and a cascade of four magnetic phase transitions at high fields, of which the first two account for the experimental data.

We present magnetization and magnetostriction studies of the insulating perovskite LaCoO3 in magnetic fields approaching 100 T. In marked contrast with expectations from single-ion models, the data reveal two distinct first-order spin transitions and well-defined magnetization plateaux. The magnetization at the higher plateau is only about half the saturation value expected for spin-1 Co 3+ ions. These findings strongly suggest collective behavior induced by strong interactions between different electronic -and therefore spin -configurations of Co 3+ ions. We propose a model of these interactions that predicts crystalline spin textures and a cascade of four magnetic phase transitions at high fields, of which the first two account for the experimental data.

PACS numbers:
Spin state crossovers induced by changes in the electronic configuration of transition metal ions can dramatically alter several materials properties [1,2]. One example is the pressure induced spin transition of ferric ions in magnesium silicate perovskite that occurs in the Earth's lower mantle [3,4]. Apart from being a major constituent of the Earth's mantle, perovskites exhibit a variety of behaviors which include colossal magnetoresistance [5,6] and high temperature superconductivity [7]. For certain insulating perovskites containing Fe or Co ions, the two lowest energy multiplets of the 3d electrons can be separated by a small energy gap, ∆, owing to competition between Hund's coupling and crystal field splitting. A external magnetic field, H, can change the electronic configuration if the lowest energy multiplet has lower spin than the first excited multiplet. In this situation, structural transitions with unusually large magnetoelastic responses can be induced for H ∆/gµ B (where µ B is Bohr's magneton and g is the g-factor).
The energy gap ∆ separating low-and high-spin electronic configurations of 3d ions in perovskites is typically of order 10-1000 meV. Controlled switching between these electronic configurations with magnetic fields therefore requires extremely large fields of order 100-10000 T (assuming g ≈ 2). Therefore, such studies have proven difficult to date. The recent development of nondestructive pulsed magnets with peak fields of 100 T and with long (ms-timescale) pulse durations now provides an opportunity to realize field-tuned spin state transitions in compounds with ∆ ∼ 12 meV. Such studies are highly desirable, as they would provide an essential complement to pressure-induced spin cross-over phenomena that have been studied in a number of materials [8][9][10]. Moreover, high-field studies allow one to test whether single-ion models suffice to describe spin-crossover phenomena in perovskites, or whether interactions between 3d ions are important.
To this end we study the insulating perovskite cobaltite (LaCoO 3 , or LCO), whose octahedrally-coordinated Co 3+ (3d 6 ) ions are natural candidates to explore fieldinduced transitions of the electronic configuration and spin [11,12]. A small energy gap ∆ ≈ 12 meV separating the spin singlet ground state (6 electrons in the t 2g orbitals) from the lowest energy magnetic configuration results from the competition between Hund's coupling and crystal field splitting [13]. While the Co 3+ ions are in their S = 0 state at low temperatures, thermal activation to a magnetic S = 0 state occurs above ∼30 K, giving rise to a paramagnetic response. Although considerable work has been carried out on this thermally-induced spin crossover [14][15][16], the spin value S of the lowest energy excited multiplet is still controversial. A single-ion model allows to explain the temperature dependence of the magnetic susceptibility and the measured g-factor (g = 3.4) with an S = 1 triplet [15]. However, field-dependent magnetization studies have recently suggested that interactions between Co ions may also play a crucial role in LCO: a field-induced gap-closing study shows a jump in the magnetization M to a plateau value of 0.5µ B just above 60 T [17]. This value is much lower than the 2µ B (4µ B ) expected from saturated S = 1 (S = 2) Co 3+ ions within a single-ion model. From now on we will assume that S = 1, because this assumption leads to a simpler explanation of the observed plateaux.
Very recently Platonov et al. [18] have reported the detection of several magnetic field-induced phase transitions in LCO single crystals at 4.2K. The samples are subjected to very high magnetic fields (≈ 500 T) produced by rapid (≈ 15 µs) explosive compression of magnetic flux. The results show that the magnetization per Co 3 + ion starts to rise significantly for H > 50T, where M/n Co 0.4µ B , reaching a plateau value of M/n Co 1.4µ B at 140 T, and a maximum value of 3.5µ B at 500 T. It is suggested that the observed plateau and other features may be linked to antiferromagnetic interactions. [18] Smoothing of the first order transition near 60 T [17] is attributed to the rapid variation of the magnetic field. No firm conclusions on the Co 3+ spin state are presented, although it is suggested that the S = 2 state may be important at the highest fields reached.
Here we measure the magnetization (M ) of LCO in longer-duration non-destructive pulsed magnetic fields approaching 100 T [19]. In contrast to expectations from single-ion models, we observe multiple magnetization steps and plateaux, which indicate the relevance of inter-ion coupling. Moreover, magnetostriction studies complement the magnetization results, and reveal large lattice changes induced by a combination of two factors: a) the S=0 and S=1 electronic configurations of the Co 3+ ions have different volume and b) the S=1 configuration is Jahn-Teller active. Susceptibility measurements made at higher fields in a single turn pulsed field experiment reveal that the second plateau persists at least up to 140 T. We propose a model that allows for collective behavior of the Co 3+ ions and reveals a rich phase diagram in which the plateaux originate from different stable crystalline textures of S = 1 Co 3+ ions in a background of S = 0 configurations. The orbital degree of freedom of the S = 1 Co ions also leads to orbital ordering that should coexist with the magnetic ordering.
The single crystals of LCO are grown in a floating zone furnace at Argonne National Laboratory, with pieces being cut and polished into either small tapered cylinders for the magnetization M measurements or 1 × 1 × 3mm bar shaped samples for the magnetostriction experiments. Magnetization curves are obtained as a function of µ 0 H, using an extraction magnetometer [20], with pulsed fields up to 97 T at temperatures in the range 1.6-55 K. High precision magnetostriction measurements are made using a fiber-optic strain gauge [21,22]. In all experiments, singe crystals are oriented with their [010] pseudo-cubic axes parallel to H. Figure 1(a) shows M versus µ 0 H at T = 1.7 K revealing two transitions, each accompanied by increases ∆M ∼ 0.5 µ B /Co, and two plateau regions consistent with two distinct magnetic phases. The second plateau occurs at M ∼ 0.9 µ B /Co which is well below the magnetization of 2 µ B /Co expected from saturated Co ions. Both transitions exhibit hysteresis in this pulsed field experiment. The location of the transition between rising and falling fields, as indicated by the blue and red arrows, points to a diffusionless, displacive transition type [23][24][25][26]. The inset in Fig. 1(a) for fields in the range 55-63 T, and with offsets for clarity, presents an expanded view of the first transition at temperatures between 1.6 K and 55 K. The hysteresis diminishes as T is raised. Once initiated, the transitions are observed to proceed rapidly, within ∼ 5 µs, over the entire sample (length ∼ 4 mm). Figure 1(b) shows the magnetostriction measured parallel to the field. The transition fields and ∆L/L plateaux closely match the magnetization behavior in Fig. 1(a). The red and blue arrows again represent rising and falling field transitions. The magnitude of the field-induced strain ∆L/L ∼ 10 −3 exceeds that found in structurally related perovskites [27]. For the rapidly rising fields, some rounding of the transitions is attributed to sample heating associated with latent heat released at the transitions, which occurs despite immersion in superfluid 4 He. The magnetostriction results, which show FIG. 2: (Color online) Phase diagram for LCO based on the results shown in Fig. 1(a). Blue and red symbols denote rising and falling fields respectively, while black symbols connected by the solid curve represent the thermodynamic phase boundary obtained as an average of the rising and falling field values at a given T . SSC1 and SSC2 denote spin state crystalline one and two respectively. The inset shows low lying energy states as a function of applied field for Co 3+ in LCO based on a single ion description [17].
sample expansion (∂V /∂H) T,P > 0 , are consistent with the Maxwell relations since (∂M/∂P ) T,H < 0 for T < 100 K [22]. The inset in Fig. 1(b) shows the magnetic susceptibility at 4 K measured in a pulsed single turn magnet which reaches ∼ 140 T. Figure 2 shows a phase diagram for LCO based on the inset curves in Fig. 1. Besides the non-magnetic (NM) phase at low fields, there are two ordered magnetic phases denoted spin state crystalline 1 (SSC1) and 2 (SSC2) respectively. The inset depicts low-lying single-ion states for Co 3+ in LCO [15,17].
The observed field induced phase transitions must involve interactions between the non-magnetic ions and the magnetic S z = +1 ionic configurations that are stabilized by the field. We note that the magnetic ions are both larger in volume as well as Jahn-Teller active [15,17]. Indeed, the S = +1 Co 3+ ion with a t 5 2g e 1 g electron configuration has an orbital degree of freedom since there is one electron shared between two e g orbitals. We assume that the energy gap (∆) between the S = 1 and nonmagnetic configurations is large in comparison to the exchange interactions between different magnetic configurations of Co 3+ ions. We also assume that the gap to the S = 2 magnetic configuration is much larger than 2∆ and consequently this configuration can be eliminated from the low-energy Hilbert space for fields up to 100 T. This assumption is supported by the consistency of the model predictions with experiment. Applying a magnetic field H lowers the energy of the doubly degenerate S z = +1 states (the e g orbitals being d 3z 2 −r 2 and d x 2 −y 2 ) such that their energy becomes comparable to the S = 0 nonmagnetic state. A low-energy effective model is obtained by projecting the original Hamiltonian into the subspace of these three lowest-energy states.
The orbital state of the S z = +1 ion is parameterized by an angle θ: |θ = cos(θ/2)|3z 2 −r 2 +sin(θ/2)|x 2 −y 2 , corresponding to a biaxial deformation of the local O 6 octahedra. The doublet orbital degeneracy manifested as a continuous symmetry for the angle θ is retained by linear Jahn-Teller (JT) interaction. In real compounds, however, the continuous degeneracy is lifted by the lattice anharmonicity and higher-order interactions [28]. In fact, only local elongations of the O 6 octahedra are observed in almost all JT ions with e g electrons [29]. Depending on the local elongation axis, the allowed orbitals are d 3x 2 −r 2 , d 3y 2 −r 2 and d 3z 2 −r 2 , corresponding to θ = 0 and ±2π/3. The large nonlinear JT distortion also lifts the degeneracy of the t 2g levels and quenches the orbital degrees of freedom associated with the hole left in the t 2g manifold of the S z = 1 state.
A single S = 1 ion embedded in the host matrix acquires additional energy via JT coupling with the local octahedron and exchange interaction with its nonmagnetic neighbors, amounting to a slight modification (∆) of the energy gap: Here n i = S z i = +1 such that n i = 1 for a magnetic configuration on the ion i,∆ is the renormalized spin state gap and gµ B H is the Zeeman interaction. The single-ion physics described by Eq. (1) implies a single field-induced crossover when gµ B H >∆. To account for the observed multiple transitions, we need to include the interactions between the S z = +1 sites: The leading order coupling (V ij ) is isotropic and repulsive because it results from the increased volume of the S z = 1 ions relative to the nonmagnetic ones. The second interaction terms are orbital-dependent and determine the relative orbital orientations of the magnetic ions. [21] The three-state Potts variable s i a indicates the orbital states (d 3l 2 −r 2 with l = x, y, or z) of the S z = 1 ion. The origin of the isotropic interaction V ij can be understood by using the so called "sphere-in-the-hole" model [28]. The larger ionic size of a S = 1 ion in the host matrix of S = 0 ions acts as an elastic impurity and creates a strain field which decays as ∼ 1/r 3 at large distances. A second impurity interacts with this strain field, giving rise to the elastic term in Equation (4).
We find that the key aspects of the data − namely multiple phase transitions, incremental steps in the magnetization, lattice expansions and metastability − can be captured by considering field-induced S z = +1 spin states in which we neglect orbital orientation-dependent terms. The Hamiltonian can then be mapped onto an effective Ising model on a cubic lattice with nearest neighbor and next nearest neighbor antiferromagnetic interactions J 1 = V 1 /4 and J 2 = V 2 /4, respectively, between σ i = ± 1 pseudospins. S = 0 becomes σ i = − 1 and S z = 1 becomes σ i = + 1 via the transformation n i = (1 + σ i )/2 and the adoption of an effective magnetic field h = 1 Fig. 3). We note the orbital orientation-dependent terms neglected in Eq.(3) are relevant for deciding the orbital ordering that accompanies each spin ordering at low enough temperatures.
Provided the nearest neighbor coupling dominates (i.e. J 1 > 4J 2 ), the h = 0 ground state h = 0 is a twosublattice Néel order, whose actual magnetization is half of the fully saturated value m = 1 2 corresponding to σ i = 0. We associate this ordering with the measured magnetization plateau in the region above ∼ 75 T shown in Figs. 2 and 1. Intermediate plateaux with m = 1 4 and m = 3 4 (in which σ i = ∓ 1 2 ) become viable for finite J 2 . On transforming back to the original Hamiltonian (considering both positive and negative h), we therefore realize a cascade of first order transitions between different phases with critical fields − H c1 , H c2 , H c3 and H c4 in Fig. 3 − whose relative separations depend on the nearest and next nearest neighbor repulsions V 1 and V 2 . Fig. 3 shows quantitative consistency between the measured magnetization steps and those of the model owing to the quenching of the orbital contribution to the gfactor (such that g ∼ 2) by a JT distortion of the O 6 octahedra encasing the S z = +1 sites. On associating the first two transitions (H c1 and H c2 ) with the two transitions observed experimentally, we estimate V 2 ≈ 1.2 K. The lack of further transitions in motor-generator-driven magnetic fields extending to 97 T (in Fig. 1), and in a single turn magnet system delivering fields to ≈ 140 T (shown in the inset to Fig. 1(b)), implies that V 1 25 K. On enumerating V 1 over the six nearest neighbor Co atoms in the cubic perovskite, we arrive at a lower limit 150 K for the energy scale over which collective effects influence the phase diagram of LCO.
Since the S z = +1 ions form a bipartite BCC structure, the dominant antiferro-orbital interaction gives rise to a staggered orbital order in the m = 1/4 plateau. On the other hand, antiferro-orbital interactions are frustrated in the second m = 1/2 plateau as the S z = +1 ions form a non-bipartite FCC lattice. Nonetheless, our Monte Carlo simulations found a partial orbital ordering with a layered structure.
The different crystalline structures in Fig. 3(b) (e.g. BCC and FCC) anticipated for each of the phases of LaCoO 3 provide an explanation for the discontinuous changes in the lattice and hysteretic behavior. Each magnetization plateau in Fig. 3(a) corresponds to a different optimal sublattice arrangement of S = 0 and S z = +1 spin states. Orbitals in the S z = +1 sublattice are more voluminous and directional in nature, causing the lattice to expand in response to the pressure exerted by their increased density (i.e. V 1 and V 2 ). The reduced expansion at the second transition (H c2 ) suggests the predominantly repulsive interaction is at least partially compensated by an attractive antiferro-orbital interaction, which is able to afford a more efficient arrangement of the orbitals at higher densities [21].
Finally, we reiterate that the spin value S of the first excited multiplet is still matter of debate [15,30,31]. While the generic form of the effective Ising Hamiltonian of Eq. (3) does not depend on the value of S, one would have to include further neighbor repulsions among magnetic configurations in order to explain a ∼ 0.5µ B magnetization plateau for S = 2 (∼ 0.5µ B is 1/8 of the saturated value in this case). Therefore, our original assumption of S = 1 leads to a simpler and more natural explanation of the measured plateaux. We also note that the orbital physics and magnetostrictive properties are likely to be very different for S = 2 owing to the e g orbitals no longer being JT active.
The experimental evidence therefore suggests that collective behavior involving two spin states leads to different crystalline arrangements that can be tuned with magnetic field. While magnetostriction associated with field-tuned orbital order has been reported in manganites [27], the present behavior in LaCoO 3 is different in that a strong coupling of spin, orbital and lattice degrees of freedom occurs in a Mott insulator. An entirely different type of functionality results: multiple field-tuned (diffusionless) transitions giving rise to a rapidly switchable strain. Spin state crystallization is likely to be a general property of crystalline materials subject to a spin state transition under extreme conditions [3,4,8,9], causing such materials to potentially become metastable in the vicinity of a phase transition and vulnerable to a sudden release of mechanical energy.
Acknowledgements. Work at the National High Magnetic Field Laboratory is supported by the National Science Foundation under DMR-0654118, Florida State University, the State of Florida and the U.S. DOE BES project "Science at 100 tesla". Work at Argonne National Laboratory is supported by the U.S. Department of Energy under Contract No. DEAC02-06CH211357.

Magnetization
Magnetization to ≈ 63 T is measured in a 'short-pulse' capacitor-driven pulsed magnet over temperatures ranging between 1.7 and 55 K using an extraction magnetometer [20] consisting of concentrically arranged counterwound coils of ≈ 1000 turns. Magnetization studies to ≈ 97 T are performed in a generator-driven nondestructive magnet at Los Alamos National Laboratory [19]. The magnetic susceptibility is measured inductively as an induced voltage, with the discontinuous nature of the transitions typically yielding a signal in the form of a delta function of height ∼ 10 V at each of the transitions. Susceptometer measurements in the 140 T singleturn magnet system [32] were made using a compensated coil consisting of a few turns.

Magnetostriction
Magnetostriction measurements are made using a fiber-optic Bragg grating (FBG) attached to the sample using cyanoacrylate bond [33]. The single-mode fiber contains a 1 mm long Bragg grating and is manufactured by SmartFibre, UK. The reflection of light at the Bragg wavelength shifts when the sample length changes. Optical spectra are collected with a fast 1024 pixel InGaAs infrared line array camera, which is read out at 47 kHz. Mounted on a fiber-coupled research spectrometer of focal length 500 mm with a 600 groove/mm −1 diffraction grating, the 25 µm pixel pitch results in a dispersion, δλ of 0.056 nm per pixel at 1550 nm. To illuminate the FBG, we use a superluminescent diode broadband source with an output power of 30 mW in the wavelength range 1525-1565 nm.

Dissipation effects
The hysteresis at each of the transitions evidences dissipation resulting in the unavoidable generation of heat during the experiments. This irreversible heating is likely partially compensated by cooling due to the magnetocaloric effect, which in turn is caused by the fieldinduced closing of the spin gap. The sample is partially able to impart the resultant heat to the surrounding cryogen medium in times of order 1 ms at base temperature (T ≈ 1.6 K) where the sample is immersed in superfluid 4 He. Owing to the rapid succession of the two transitions on the rising magnetic field where the rate of change of magnetic field is highest, not all of the heat generated by the first transition is able to escape the sample before the second transition, possibly contributing to the temperature-rounding of the second transition on the rising magnetic field. The slower ramp rate of the falling magnetic field enables more heat to dissipate to the cryogen between the transitions, enabling both transitions to be clearly observed with reduced rounding.

Maxwell relation and field-induced volume expansion
The Maxwell relation ∂V /∂H| T,P = −∂M/∂P | T,H , where V is the sample volume, H is the applied magnetic field, M is the sample magnetization, and P is the applied hydrostatic pressure, provides a means to check our data for self consistency. We know from prior results [34] that ∂M/∂P < 0 for temperatures T < 100K, caused by the increase of the crystal field gap with pressure. Hence, an increase in the sample volume with magnetic field is expected. This is consistent with our experimental observation of the longitudinal magnetostriction in high fields, and what we expect from our model.

Background
In the model, we assume the single-ion ground state of Co 3+ ion to be non-magnetic (S = 0) while the first excited state is magnetic with S = 1. The latter has an orbital degree of freedom since there is one electron shared between two e g orbitals. We further assume that the energy gap (∆) between the S = 1 and nonmagnetic configurations is large in comparison to the exchange interactions between different magnetic configurations of Co 3+ ions. Applying a magnetic field H lowers the energy of the doubly degenerate S z = 1 triplet states (the e g orbitals being d 3z 2 −r 2 and d x 2 −y 2 ) such that their energy becomes comparable to the S = 0 nonmagnetic state. A low-energy effective model is obtained by projecting the original Hamiltonian into the subspace of the three low-energy states.

Relevant terms in the Hamiltonian
A single S = 1 ion embedded in the host matrix acquires additional energy via JT coupling with the local octahedron and exchange interaction with its nonmagnetic neighbors, amounting to a slight modification (∆) of the energy gap (∆) entering into the single ion component (H 1 ) of the Hamiltonian given by Eq. (1) in the main text. As the density n i of S = 1 states increases, the interactions between them leads to the formation of spin-state superstructures and orbital order, necessitating a second component of the Hamiltonian describing interactions between different sites. H 2 is the same as Eq. (2) in the main text. The dominant mechanism is the elastic interaction mediated by longwavelength phonons. Short-range interactions are usually orbital dependent and have contributions from both superexchange and JT couplings. The origin of the isotropic (V ij ) term is discussed in the main text. The orbital-dependent interaction potential in Eq. (4) consists of two terms: where the ferro-and antiferro-orbital bond operators are given by O a si,sj = δ si,αij 1 − δ sj ,αij + 1 − δ si,αij δ sj ,αij . Since the orbital interaction is generally anisotropic, the variable α ij = x, y, z indicates the orbital flavor (orientation) which has distinct symmetry properties with respect to the other two on a given bond (ij). For example, α ij = z (d 3z 2 −r 2 orbital) on a nearest-neighbor bond along z-axis or a next-nearest-neighbor bond in the xy plane. Orbital configurations in which the α-orbital is empty on both sides are neglected in Eq. (7) since such terms are irrelevant for orbital ordering.
The orbital-dependent interaction potentials (V f ij and V a ij ) result mostly from the anisotropic superexchange and JT interactions, as described in more detail in the following section. Since the elongation of O 6 octahedra is the source of local anisotropic stress, the resultant phonon-mediated interaction is also orbitaldependent [35] thereby contributing to both V f ij and V a ij . The interaction potentials for two S z = +1 ions on the nearest neighbor and next-nearest neighbor bonds are summarized in Fig. 4. The electronic contribution to the interaction potentials becomes negligible beyond a few nearest neighbors as higher-order hoppings are required for such virtual exchange processes. The orbitaldependent interaction potentials (V f ij and V a ij ) result mostly from the anisotropic superexchange and JT interactions, as described in more detail in the following section. Since the elongation of O 6 octahedra is the source of local anisotropic stress, the resultant phonon-mediated interaction is also orbital-dependent [35] thereby contributing to both V f ij and V a ij . The interaction potentials for two S z = +1 ions on the nearest neighbor and next-nearest neighbor bonds are summarized in Fig. 4. The electronic contribution to the interaction potentials becomes negligible beyond a few nearest neighbors as higher-order hoppings are required for such virtual exchange processes.

The effective Hamiltonian
While nearest-neighbor ferro-and antiferro-orbital interactions are expected to be significant, these terms do not play a role in the formation of the first two fieldinduced phases. The density of S z = +1 orbital states is too low to support the existence of such bonds. It is therefore reasonable to assume that the net interaction is repulsive, with larger orbitals tending to repel each other, enabling us to construct an effective Ising model Hamiltonian (H eff ) consisting of nearest neighbor (V 1 ) and next-nearest neighbor (V 2 ) repulsive terms as given by Eq. (3) in the main text.
The magnetization curve m(H) of the original model is shown in Fig. 3a in the main text. The width of the m = 1/4 and 3/4 plateaus is ∆H 1/4 = ∆H 3/4 = 24J 2 = 6V 2 , indicating that both states are stabilized by secondnearest neighbor interactions.

Orbital ordering
In the m = 1/4 plateau phase, likely corresponding to the first field-induced phase (SSC1) in Fig. 3b in the main text, the closest distance between S z = +1 ions is r = √ 3a (i.e. third-nearest neighbors). Since superexchange and JT interactions are negligible at this distance, we resort to the anisotropic elastic interaction discussed in Ref. 35 to calculate the effective interaction potentials: where Γ = (3c 11 + 4c 44 )/(c 11 + 2c 44 ), c 11 and c 44 are elastic moduli of the crystal, and σ 0 is the strength of the local stress. Since the third-neighbor bonds are parallel to the 111 directions, the three 3l 2 − r 2 orbitals (l = x, y, z) are equivalent with respect to the bond. The corresponding bond operators are simply O f si,sj = δ si,sj and O a si,sj = 1 − δ si,sj . The elastic interaction given by Equation (7) shows an additional energy cost for ferroorbital configurations, we therefore expect the S z = +1 ions at the m = 1/4 plateau to form a bipartite bodycentered cubic (BCC) lattice with a doubled lattice constant and a staggered orbital order. Since there are three possible orbital orientations, there are six different staggered orbital configurations that have the same energy in absence of spin-orbit coupling: d 3ν 2 −r 2 , d 3ν 2 −r 2 with ν, ν = {x, y, z} and ν = ν .
In the m = 1/2 plateau phase associated with the second field-induced phase (SSC2) in Fig. 3b in the main text, the next-nearest neighbor antiferro-orbital component (V a 2 in Fig. 4) is expected to be relevant, which has a contribution from superexchange. While indirect hopping via O ions only contributes to sixth order to the exchange energy between next-nearest-neighbor Co 3+ ions, direct hopping between t 2g orbitals contributes to secondorder. On taking into account the splitting δ JT of t 2g levels due to the elongation of the local octahedron, we where t ddσ is the direct ddσ hopping between t 2g orbitals and U dd is the on-site Coulomb repulsion of d electrons.
Unlike the case of m = 1/4 plateau, The antiferro-orbital interaction V a 2 remains frustrated as the S z = 1 ions in the half-plateau form a non-bipartite face-centered cubic (FCC) lattice. Nonetheless, our Monte Carlo simulations found a partial orbital ordering with a layered structure. Detailed studies of the orbital order will be presented elsewhere.

Exchange interaction
Here we compute the exchange energies for a pair of Co 3+ ions on a nearest neighbor bond. Figure 5 shows six inequivalent configurations of Co ions on a nearestneighbor bond along the z-axis. The directional dependence of e g orbitals results in a highly anisotropic electron hopping, which in turn leads to an anisotropic exchange interaction. To simplify the calculation, we consider only the dominant pdσ hopping between e g and p orbitals of the O 2− ions. The resultant exchange energy of a z-bond thus depends on whether the d 3z 2 −r 2 orbitals are occupied. We expect E a = E b = E d = 0 , E c = E e = 1 and E f = 2 . Here the subscripts of indicate the number of occupied d 3z 2 −r 2 orbitals on the bond. For simplicity, we also neglect Hund's coupling.
We employ standard degenerate perturbation theory to calculate the bond energies up to fourth order in t pd . For configurations a, b and d in Fig. 5, the second and fourth order corrections to the bond energy are and Here, U pp and U dd are the Coulomb repulsion energies for double occupation of the p z and d 3z 2 −r 2 orbitals respectively, while ∆ pd is the diagonal energy difference between the two orbitals.

Jahn-Teller interaction
We compute the energy gain of a bond via Jahn-Teller coupling to the O 6 octahedra. For simplicity, we adopt the so-called Einstein phonon model for the displacement u of the O ion, namely E e = k|u| 2 /2. For a bond oriented along the z-axis, the JT couplings for the d 3z 2 −r 2 and d x 2 −y 2 orbitals are γu z and −γu z /2, respectively [36], where γ is a coupling constant. By integrating out the displacement u, the energy gain for the six different configurations are

Calculation of interaction potentials
Here we derive the effective interaction potentials between a pair of Co 3+ ions. In particular we will focus on the exchange contributions to the potential energy. We introduce a pseudospin 1/2, τ , to describe the orbital degrees of freedom of the S z = +1 ions such that |τ z = ±1/2 correspond to occupied orbital |3z 2 − r 2 and |x 2 − y 2 , respectively. As can be seen from the previous discussion, the exchange energy of a given bond depends only on the occupation of the active orbital (e.g. d 3z 2 −r 2 on a z-bond). We thus introduce an operator P ν for the three different nearest neighbor bonds (ν = x, y, or z) which projects onto the relevant orbital state. On introducing the following unit vectorŝ the bond projection operators become P ν = 1 2 + τ ·n ν . The exchange Hamiltonian is then given by H ex = ν ij ν n i n j K 1 (τ i ·n ν )(τ j ·n ν ) + K 2 (τ i + τ j ) ·n ν + K 3 + 0 (1 − n i )(1 − n j ) +G 1 [n i (1 − n j ) + (1 − n i )n j ] + G 2 [n i (1 − n j ) τ i ·n ν + (1 − n i )n j τ j ·n ν ] .
Here K 1 = 2 + 0 − 2 1 , K 2 = ( 2 − 0 )/2, K 3 = ( 2 + 0 + 2 1 )/4, G 1 = ( 1 + 0 )/2, and G 2 = ( 1 − 0 ). In the regime dominated by nonlinear JT interactions, quantum fluctuations of the orbitals are quenched by the lattice distortion. We thus replace the pseudospin operator by a classical O(2) vector of length |τ | = 1/2. The local elongation of the O 6 octahedron aligns the τ i vector to one of the three directionsn ν . Each corresponds to the three Potts states discussed in the main text. Noting that τ i ·n ν = 1/2 for s i = ν and −1/4 otherwise, the effective nearest neighbor interaction potentials in Equation (4) can be easily read off from H ex as where the exchange constant J is given by and the correction to the gap energy is