Quantum order by disorder is a key to understanding the magnetic phases of BaCo2(AsO4)2

Introduction

The interplay between crystal field symmetry, spin-orbit coupling, and the electron configuration of cobalt ions (Co²⁺, 3d⁷) in a nearly ideal octahedral oxygen environment stabilizes the J_eff = 1/2 state, a fundamental building block of numerous intriguing physical phenomena in various configurations. One of the intriguing theoretical proposals is that cobaltates with a J_eff = 1/2 state can potentially exhibit bond-dependent exchange interactions, including what is lately called the Kitaev exchange interaction^1,2,3,4. It is well known that when the Kitaev exchange interaction dominates the spin Hamiltonian, the ground state can become a Kitaev quantum spin liquid (KQSL) due to strong magnetic frustration^5,6,7. The KQSL is highly sought after among frustrated magnets because it has the potential to host exotic emergent phenomena such as Majorana Fermions and topological order, with significant applications in quantum computing^8,9,10. Therefore, several honeycomb cobalt-based magnets, including Na₂Co₂TeO₆^11,12,13, and Na₃Co₂SbO₆¹⁴ have been investigated intensively.

Recently, another cobalt-based honeycomb-lattice material BaCo₂(AsO₄)₂ (BCAO) has been proposed as an ideal platform for realizing a KQSL^15,16 (Fig. 1a). BCAO shows a series of magnetic phase transitions under applied magnetic field (H). The zero-field order is identified as either an incommensurate spiral or double-stripe (DS), which is suppressed in favor of a collinear up-up-down (UUD) phase under field^17,18,19,20. These spin structures are illustrated in Fig. 1b. For magnetic fields in the crystallographic ab-plane, the latter phase is subsequently suppressed by relatively weak fields (0.5 T), above which a KQSL phase has been suggested^15,21. For field in the out-of-plane direction, which is theoretically more conducive to a field-induced KQSL phase^22,23,24,25, a broad magnetic continuum was observed¹⁶. However, conflicting reports exist regarding whether BCAO is best described by an extended Kitaev model with strongly bond-dependent magnetic couplings or as an XXZ-J₁–J₃ model with weakly bond-dependent couplings; recent inelastic neutron scattering measurements, as well as ab initio calculations²⁶, align more with the XXZ-J₁–J₃ model²⁷. It has also been pointed out that strong distortions of oxygen octahedra, whose the trigonal term in crystal-field splitting matrix is larger than one fifth of spin-orbit coupling, lead to Ising- and XY-like spin behavior^4,28,29. Additionally, the weak ligand-assisted exchange interactions weaken bond-dependent exchange interactions⁴. Hence, further investigation is still required to characterize the phases of BCAO comprehensively. Investigating the angular dependence of phase boundaries is an effective method to understand the origin of magnetic anisotropy³⁰. Recently, Safari et al. explored the phase diagram of BCAO using magnetotropic susceptibility measurements and theoretical models³¹. Their findings indicate that while classical models qualitatively describe the system’s behavior, quantum corrections are essential for precise critical field predictions³¹.

Quantum order by disorder is a key to understanding the magnetic phases of BaCo2(AsO4)2 — **Fig. 1: Crystal and magnetic structures.**

In this study, we explored the details of magnetic field-induced ground states in BCAO using ac and dc magnetic susceptibility, capacitance measurement, and torque magnetometry measurements at various temperatures (T). With the aim of establishing the origin of the magnetic anisotropy, we study the angular dependent T–H phase boundaries by tuning the magnetic field orientation. With the application of a magnetic field, multiple-phase transitions were identified for H along the in-plane direction, consistent with previous studies^{15,16,17,18,19,20}. In comparison, for H close to the c-axis, the critical field for suppression of ordering is drastically increased, and intermediate phase transitions disappear or are washed out at the lowest temperatures.

To understand the full details of the phase diagram, we employ first-principles-based calculations to provide a starting point for analyzing the magnetic couplings. The resulting model is refined by fitting the experimental phase diagram employing second-order energy corrections³² to account for quantum fluctuation effects on the stability of various magnetic phases. This combined approach allows for the precise reproduction of the experimental observations. Consistent with recent assertions^26,31,33,34, we conclude that quantum order-by-disorder (QOBD) effects dominate in BCAO, influencing both the critical field values and specific progression of field-induced phases.

Results

Magnetic properties and phase diagrams for in-plane fields

Figure 2a illustrates the magnetization measured at various temperatures as a function of the field applied along the a-axis. Hysteresis loops emerge between 0.1 T (H₁) and 0.2 T (H₂), which become smaller with increasing temperature and are almost invisible above 4.0 K. Beyond 0.5 T, there is a steep increase in magnetization before reaching saturation around 0.6 T (H₃). The saturation magnetization is about 2.3 μ_B/Co²⁺, from which we extracted g_ab ~ 4.6. The magnetization between the saturation field and ~ 0.15 T constitutes a one-third magnetization plateau. These findings are consistent with prior studies^15,19,20. When the temperature exceeds T_N ~ 5.5 K, all features disappear, and a monotonic increase of the magnetization with field is observed, consistent with a thermal paramagnetic state. Figure 2b shows the ac magnetic susceptibility as a function field. The temperature dependent ac magnetic susceptibility is in Supplementary Fig. 1. Three peaks are observed: The first two peaks at ~ H₁ and ~ H₂ are consistent with the jump in magnetization in the hysteresis loop and the end of the hysteresis loop. The third peak above ~ H₃ corresponds to the saturation field. ac susceptibility at 6 K does not show any anomalies. The magnetic phase diagram with field along the a-axis deduced from these measurements is presented in Fig. 2c. We observed three thermodynamic phases (I, III, and IV) and one intermediate state (II) due to the hysteresis loop^17,35. Neutron scattering experiments identified phase I as the DS phase or closely related spiral phase with in-plane wavevector close to Q = (1/4, 0)^20,27. Phase III is the UUD phase with 1/3 saturation magnetization and Q = (1/3, 0)^18,19,20. Given the differing wavevectors and the discontinuous change in wavevector through this phase transition, it is not surprising that these phases are separated by a first order phase transition. Phase II corresponds to a coexistence region of Phase I and III. Lastly, phase IV is the asymptotically polarized paramagnetic phase that is smoothly connected to high temperature paramagnetic phase.

Figure 2d–f show analogous results for fields along the b^*-axis. For this field direction, we observe the same progression of phase transitions. The hysteresis loops are wider but less pronounced, and H₃ is slightly lower in comparison with H∣∣a. Therefore, the width of the 1/3 plateau (phase III) is slightly narrower with magnetic field along the b^*-axis.

Magnetic properties and phase diagrams for out-of-plane field

Figure 3a shows the temperature dependent dc magnetic susceptibility (M_c/H (T)) measured at various constant external magnetic fields. At 1 T, a clear phase transition from the paramagnetic to DS state is observed and marked with an empty circle around 5.3 K. Raising the field to 14 T only suppresses T_N to 4 K. The phase transition is therefore very robust against out-of-plane magnetic field, which is in sharp contrast with the behavior for in-plane fields.

**Fig. 3: Magnetic properties for out-of-plane fields.**

Figure 3 b shows the magnetization versus out-of-plane field up to 14 T. At 2 K, the magnetization linearly increases without showing any anomaly, demonstrating a gradual spin canting towards the field direction without phase transitions. At 4 K, we observe a weak slope change around 14 T, denoted as a black arrow, which is consistent with the thermal phase transition between the DS and thermal paramagnet observed in Fig. 3a. With increasing temperature to 4.5 K, the slope change shifts to slightly lower field of about 11.5 T. Suppression of the phase boundary between AFM and paramagnetic phases with increasing temperature is a typical behavior of antiferromagnets.

Since 14 T is not enough to observe the magnetic saturation, we measured the magnetization with field along c-axis up to 50 T at Pulsed Field Facility in Los Alamos National Laboratory. Results from pulsed field magnetization experiments at various temperature are shown in Fig. 3c, d. The data are admittedly noisy, but the important features are clearly discernible. In the paramagnetic state at 10 K, the magnetization curve resembles the Brillouin function. Below T_N, the magnetization increases linearly and shows a slope change in the vicinity of 25 T, which we assign as the out-of-plane saturation field. Therefore dM_c/dH is a constant below ~ 20 T and above ~ 35 T. Due to the small crystal field excitation gap of Co²⁺, it is common to observe a linearly increasing magnetization above the saturation field, due to Van Vleck paramagnetism^36,37. By extrapolating the Van Vleck contribution to the zero field, we found a saturation moment of about 0.71 μ_B/Co²⁺, from which we extract g_c ~1.5. We note that the saturation field with out-of-plane configuration in this study is higher than the field claimed by Zhang et al.¹⁶ and Safari et al.³¹. As we demonstrated in Supplementary Fig. 4 and 5, a tiny misalignment of the crystal can lead to a drastically suppressed saturation field. For example, a field only 6° degrees off from the c-axis suppresses the saturation field 27 T to 4.5 T. The 25 T phase transition is smoothly connected to the zero field AFM to paramagnetic phase transition, additionally confirming that it is the saturation field. Thus, the out-of-plane magnetization measurements demonstrate that BCAO is highly anisotropic between in-plane and out-of-plane fields, with a saturation field ratio ({H}_{sat}^{c}/{H}_{sat}^{a}) of about 50.

Angular dependence of magnetic phase transitions

To investigate how the phase boundaries evolve with the field direction, we measured the magnetization, capacitance, and magnetic torque below T_N while varying the magnetic field orientation across the c–b^* plane. In Fig. 4a, we first show the evolution of the magnetization. For an out-of-plane field (θ = 0°), the magnetization does not show any features up to 14 T. However, with a small rotation of the field direction towards the b^* direction, the transitions to the intermediate UUD and polarized phases appear. At θ = 6. 3°, the phase transitions to the UUD and polarized phases occur at 1 T and 4 T, respectively. These move towards ~ H₁ and ~ H₂ as the field direction approaches b^*-axis, consistent with Fig. 2d. The magnetic moment also monotonically increases as the field is rotated towards b^* direction, consistent with the anisotropy of the g-values.

**Fig. 4: Field-angle dependence of the magnetic properties.**

Capacitance and magnetostriction measurements can be sensitive probes of magnetic phase transitions as the change in spin structure accompanies the change in electric susceptibility and lattice parameter of magnetic insulators. We found this to be the case for BaCo₂(AsO₄)₂. Figure 4b shows the change in capacitance as a function of magnetic field at various angles. We applied the electric field at 16 kHz along b^*-axis. At 27°, the capacitance change shows clear two phase transitions corresponding to the phase transition from DS to UUD, and UUD to polarized phase. They evolve to higher fields as the field becomes close to c-axis. The anomaly from the phase transition from DS to UUD phase becomes weaker while that from UUD to polarized phase is more or less similar at all angles. The capacitance change and magnetostriction data with in-plane magnetic field can be found in Supplementary Figs. 2 and 3.

Figure 4c and d show magnetic torque (τ) and its differentiation with magnetic field (dτ/dH) up to 30 T as a function of magnetic field at different fixed angles. Both τ and dτ/dH display a clear anomaly at phase transitions. When the field is only at 0.36° from the c-axis, a sharp increase in τ and a peak dτ/dH appear at about 23 T. This is reasonably well consistent with the saturation field observed in pulsed field magnetization shown in Fig. 3c. The anomaly shifts quickly to lower field as the angle becomes larger than 1° and ends up being only 2 T at 19.7°.

In Fig. 5a, we show the magnetic torque measured for continuous field-angle sweeps between the c and b^* axes at T = 0.5 K at various fixed magnitudes of H. All curves exhibit 180° periodicity and a saw-tooth shape due to the strong magnetic anisotropy. The sharp change of sign of the torque near 0° indicates the c-axis is the hard axis. Additional step-like anomalies are observed near θ = 0° marking the first order transitions H₁ and H₃. These two features are only clearly observable up to 10 T and they merge into one above 10 T.

**Fig. 5: Experimental and theoretical angle-dependence magnetic torque and phase diagrams.**

In general, the magnetic torque may be expanded as³⁰:

$$frac{tau }{H}=mathop{sum }limits_{n=1}^{infty }{c}_{2n}sin (2ntheta )={c}_{2}sin (2theta )+{c}_{4}sin (4theta )+cdots$$

(1)

The experimental field-dependence of (leftvert {c}_{2}rightvert) and (leftvert {c}_{4}rightvert) is plotted in Fig. 5c. Both coefficients increase steeply up to 4 T and then slowly decrease above 4 T up to the maximum field. The ratio of c₂ to c₄ remains more or less constant at ~0.4 across the entire field range due to the dominant saw-tooth curve shapes. In the limit of large field, such that the majority of the field angles correspond to the polarized phase, the c_2n coefficients have two contributions. The first contribution is field-independent, and depends exclusively on the g-anisotropy. For the implied range of g-values (g_ab ~4.6, g_c = 1.5−2.5), the limiting ratio for this contribution leads to ∣c₄/c₂∣ = 0.17−0.27. The next order contribution to the c_2n scales as 1/H, and is related to both the XXZ exchange anisotropy and g-anisotropy³⁰. The findings that c₂ and c₄ decay with increasing field, and that their ratio exceeds the g-only value both demonstrate significant exchange anisotropy is necessary to explain the observed torque. Below, we develop an exchange model for BCAO, and use these observations to verify the model.

From the phase transitions observed in the angular dependence of magnetic properties in Figs. 4 and 5a, an angular dependence of the magnetic phase diagram is constructed, as shown in Fig. 5e. The red empty symbols represent the phase transition from the DS to the UUD phase (H₁), and the black empty symbols indicate the saturation field (H₃). As the field direction moves from the b^*-axis to the c-axis, both phase boundaries sharply increase below 10°. The saturation field reaches above 23 T, whereas the phase boundary between the DS and the UUD merges with the saturation field within 5° above about 10 T.

Microscopic magnetic model

To first define a realistic range of magnetic couplings, we employ the des Cloizeaux effective Hamiltonian (dCEH) approach^38,39,40, based on exact diagonalization of the electronic d-orbital Hamiltonian for clusters include one and two Co sites [see Methods].

We first estimate the g-tensor by diagonalizing the electronic Hamiltonian on one Co site, and projecting the Zeeman operator into the low-energy space of j_1/2 states. This yields g_ab = 5.1 and g_c = 2.4, which further confirms accurate reproduction of the trigonal crystal field splitting in BCAO. To estimate the magnetic couplings, we exactly diagonalize the electronic Hamiltonian on pairs of sites and project the resulting low-energy Hamiltonian onto j_1/2 states. The relatively low symmetry crystalline environment allows for a large number of independent exchange constants. In particular, in the (Rbar{3}) space group, each first and third neighbor bond is constrained only by inversion symmetry, allowing six independent exchange parameters per bond:

$$begin{array}{rcl}{mathcal{H}}&=&mathop{sum}limits_{ij}{J}_{ij}^{xy}({S}_{i}^{x}{S}_{j}^{x}+{S}_{i}^{y}{S}_{j}^{y})+{J}_{ij}^{z}{S}_{i}^{z}{S}_{j}^{z}\ &&+,2{J}_{ij}^{pm pm }({S}_{i}^{x}{S}_{j}^{x}-{S}_{i}^{y}{S}_{j}^{y})+{Gamma }_{ij}^{yz}({S}_{i}^{y}{S}_{j}^{z}+{S}_{i}^{z}{S}_{j}^{y})\ &&+,{Gamma }_{ij}^{xz}({S}_{i}^{x}{S}_{j}^{z}+{S}_{i}^{z}{S}_{j}^{x})+{Gamma }_{ij}^{xy}({S}_{i}^{y}{S}_{j}^{x}+{S}_{i}^{x}{S}_{j}^{y})end{array}$$

(2)

where the bond-dependent local x-axis is oriented parallel to the bond, and the z-axis is parallel to the crystallographic c-axis. Second neighbor couplings tend to be weaker and are ignored in first approximation. By varying the strengths of the Coulomb interactions within a realistic range [see Methods], we establish reasonable ranges for the first and third neighbor couplings: ({J}_{1}^{xy} sim -8) to −12 meV, ({J}_{1}^{z}/{J}_{1}^{xy} sim 0.5) to 0.6, ({J}_{1}^{pm pm } sim -0.6) to −0.8 meV, ({Gamma }_{1}^{yz} sim 0.35) to 0.5 meV, ({J}_{3}^{xy}/{J}_{1}^{xy} sim -0.25) to −0.4, and ({J}_{3}^{z}/{J}_{3}^{xy} sim 0.2). The remaining couplings are all found to have magnitude smaller than 0.1 meV.

There are several aspects of the computed couplings that bear mentioning. In BCAO, the local trigonal distortion of the CoO₆ octahedra has a significant impact on the magnetic Hamiltonian, inducing both the large g-anisotropy, and bond-independent XXZ anisotropy of the intersite exchange. Microscopically, the third neighbor couplings are dominated by a single exchange path involving the e_g electrons, which leads to antiferromagnetic interactions that are rendered anisotropic only by the crystal field effects. As such, the off-diagonal anisotropies should be nearly vanishing and the third neighbor interactions should be parameterized by a pure XXZ anisotropy, with ({J}_{3}^{xy}) being the same sign as ({J}_{3}^{z}). For the first neighbor couplings, the dominant anisotropy is also of the bond-independent XXZ type. While Co²⁺ materials can, in principle, host strongly bond-dependent couplings^1,2,3,4, in BCAO we find that these are suppressed by a combination of factors. First, the large trigonal distortion from a pure octahedral environment significantly modifies the composition of the J_eff = 1/2 moments, and imposes significant XXZ anisotropy on all couplings instead of bond-dependent exchange interactions. Second, the largest contribution to the first neighbor exchange comes from ferromagnetic “charge-transfer” contributions arising from hybridization of the Co e_g orbitals with the oxygen p-orbitals. Unlike the couplings between t_2g electrons, these interactions do not display large bond-dependent anisotropies. As such, the bond-dependent anisotropies (({J}_{1}^{pm pm }) and ({Gamma }_{1}^{yz})) are subleading. Although some of these features are found in previously proposed models for BCAO^26,27,31,33, not all are reproduced. This raises the question of whether the phase diagram can be understood by employing a model that follows more closely these microscopic considerations.

To address this question, we made a comprehensive search of the parameter space around the predicted range, with a focus on fitting the experimental phase diagram. BCAO necessarily lies in a parameter region with multiple competing phases, making it necessary to consider quantum corrections to the classical state energies. We employ the method outlined in ref. 32 [see Methods], which considers states defined by classical spin vectors, but accounts for local fluctuation effects on the state energies at second order in perturbation theory. This naturally includes the leading quantum order-by-disorder (QOBD) effects. From this search, we propose a simplified, but extremely finely tuned, six-parameter model: ({J}_{1}^{xy}=-11.75), ({J}_{1}^{z}=-5.89), ({J}_{1}^{pm pm }=-0.33), ({Gamma }_{1}^{yz}=-0.15), ({J}_{3}^{xy}=+3.37), and ({J}_{3}^{z}=+0.79) meV. The parameters are more consistent with microscopic considerations but otherwise broadly compatible with previously proposed models for BCAO with dominant XXZ anisotropy^26,31.

At zero-field, it is instructive to compare the classical and second order corrected quantum energies (denoted E_Cl and E_QM, respectively) for the six-parameter model. For this purpose, we consider the double stripe (DS), ferromagnetic (FM), Zigzag antiferromagnetic, and Q = (1/4, 0) spiral phases pictured in Fig. 1b. We find classical energies of ({E}_{{rm{Cl}}}^{{rm{FM}}}=-3.145) meV, ({E}_{{rm{Cl}}}^{{rm{Spiral}}}=-3.131) meV, ({E}_{{rm{Cl}}}^{{rm{DS}}}=-3.022) meV, and ({E}_{{rm{Cl}}}^{{rm{Zigzag}}}=-2.899) meV per site. The spread of classical energies is small compared to the overall scale of the couplings, which provides the possibility that the ground state is selected by quantum fluctuations. Indeed, E_QM is instead minimized for the collinear DS phase, with spins oriented in the plane. The other collinear spin configurations are also competitive in energy, with ({E}_{{rm{QM}}}^{{rm{DS}}}=-3.465) meV, ({E}_{{rm{QM}}}^{{rm{Zigzag}}}=-3.459) meV, and ({E}_{{rm{QM}}}^{{rm{FM}}}=-3.408) meV per site. In contrast, the spiral phase is not energetically competitive once quantum corrections are included, and does not represent a local minimum of E_QM. This necessitates estimating the corrected energy based on the classical spin orientations, for which we find ({E}_{{rm{QM}}}^{{rm{Spiral}}}=-3.279) meV per site. In agreement with refs. ^33,34, we find including quantum fluctuations favors collinear phases, which explains the relative destabilization of the non-collinear spiral phase.

In Fig. 6, we show the computed evolution of the second order corrected state energies and magnetization per site for in-plane fields. We evaluate m = ∣∑_ig ⋅ S_i∣ from the classical spin vectors for which the second order energies correspond. Under applied field, there are two first order transitions. The first is to the magnetization plateau UUD state, and the second to the polarized state. We estimate for H ∣∣ a: H₁ = 0.16 T, H₃ = 0.55 T, and for H ∣∣ b^*: H₁ = 0.18 T, H₃ = 0.52 T. Thus, the model reproduces the slight anisotropy in the in-field critical fields observed in experiment. The predicted magnetization as a function of field is shown in Fig. 6c and d, and reproduces well the experimental magnetization curves in Fig. 2.

**Fig. 6: Energy per site calculations for the energetically competing phases.**

In Fig. 5e, f, we show a comparison of the experimental and theoretical phase diagram for fields oriented in the b^*c-plane. For strictly out-of-plane fields H∣∣c, the model reproduces the direct transition from the DS state to the polarized phase. The critical field depends very strongly on the field angle; for θ = 0°, we compute H₃ = 36 T, while for θ = 1°, we compute H₃ = 18 T. This range is comparable to the broadened saturation observed in the pulsed-field measurements in Fig. 3. The sharp increase of the critical fields in the vicinity of the c-axis is mainly due to the large bond-independent XXZ anisotropy, which leads to an easy ab-plane. Overall, the model reproduces the phase boundaries well, confirming that the specific progression of field-induced phases can be understood within a microscopically compatible model.

Finally, to further check that the model reproduces the characteristic magnetic anisotropy of BCAO, we simulated the magnetic torque data from the second-order corrected energies, where τ = ∂E/∂θ. Results are shown in Fig. 5b alongside the experimental data. The simulations reproduce very well both the positions and magnitudes of the first order jumps in τ corresponding to the various phase boundaries, as well as the rounding of τ(θ) at high field. A precise comparison can be made on the basis of the coefficients c_2n. In Fig. 5d, we show the theoretical evolution of these coefficients with field. The finding that c₂ and c₄ are well reproduced by the model implies the relative magnitudes of both the g-anisotropy and XXZ exchange anisotropy are well captured.

Overall, we conclude that the microscopically motivated model provides a very good reproduction of the experimental phase diagram of BCAO and therefore may serve as a starting point for future refinement. The present results confirm the validity of recently proposed models with primarily bond-independent XXZ anisotropy^26,27,31, but show that the complex phase diagram of BCAO illuminated by various recent studies is compatible with the microscopic expectations for this material.

Discussion

As a consequence of the exchange anisotropy and frustration between first and third neighbor couplings, BCAO lies in a region of the phase diagram where multiple classical phases are energetically competitive. This allows QOBD effects to select the ground state. Over the region of couplings we considered for BCAO, we observe several trends when comparing classical state energies with the second order corrected energies. First, primarily due to the strong XXZ anisotropy, both antiferromagnetic and polarized (ferromagnetic) states have finite energy corrections. Similar to the DMRG studies in refs. ^31,34, for the physically relevant magnitudes of the XXZ anisotropy, we find that fluctuations tend to stabilize the polarized phase over the other phases. Overall, we find that quantum fluctuations tend to reduce the critical fields by favoring the polarized state, consistent with the conclusions of³¹. Second, in the range of parameters we considered, the collinear DS and UUD phases only appear as ground states due to QOBD effects. We did not find any realistic parameters for which these are classical ground states for any magnetic field. Conversely, the spiral phase was not found to be stable once fluctuation effects were included. These findings are consistent with^33,34. Finally, we emphasize that the present model for BCAO is extremely finely tuned. As can be anticipated from the (necessarily) small differences of the state energies in Fig. 6a, b deviations of any coupling on the order of ± 0.01 meV may be sufficient to ruin the agreement with experiment. This sensitivity likely also leads to variations in the results of different numerical approaches. In this work, we opted to use the inexpensive approximation of second-order energy corrections to facilitate a fine parameter search. Future refinement of the model could employ more sophisticated numerical methods, provided those approaches adequately treat the QOBD effects.

These findings raise the question: which features of our proposed model are necessary to reproduce the experimental phenomenology of BCAO? Both microscopically and experimentally, there is strong justification for a J₁–J₃ model with J₁ < 0, J₃ > 0 and significant easy-plane bond-independent XXZ anisotropy (see e.g., ²⁷). Provided leading order effects of quantum fluctuations are treated, we find that these considerations are sufficient to reproduce much of the phenomenology of the BCAO phase diagram. This includes the 1/3 magnetization plateau for in-plane fields, and the lack of this plateau for out-of-plane fields. However, the additional variations of the critical fields with in-plane field direction, highlighted in the present work, point to additional bond-dependent anisotropic couplings. Considering the largest of such bond-dependent couplings, as identified by first principles calculations, allows for full reconciliation of the phase diagram and critical field values.

In conclusion, this work is a comprehensive study of the magnetic field-induced ground states of BCAO, employing various techniques such as dc magnetization, capacitance, ac magnetic susceptibility at low temperatures. We mapped out the full T–H phase diagram for different magnetic field orientations, including H along a, b^*, c, and ({c}^{{prime} }). To understand the specific progression of field-induced phases, we employed first-principles-based calculations to derive a model that reproduces all essential features of the low-temperature phase diagram. The stabilization of two colinear phases (I and III) in both zero and finite fields is attributed to the significant influence of QOBD effects. BCAO presents an intriguing material where quantum fluctuation effects dominate the phase diagram at both zero and finite field. This investigation highlights new avenues for understanding the complex behavior of BCAO under external magnetic fields.

Methods

Experimental details

High-quality single crystals of BaCo₂(AsO₄)₂ synthesized through a flux method at Oak Ridge National Laboratory in Oak Ridge, for which the detailed methods are described elsewhere²⁷. At low temperature, dc magnetization (M) is measured up to 60 T through millisecond timescales in pulsed magnet using custom-built ³He pumping (400 mK-3 K), and ac magnetic susceptibility (χ_ac) and capacitance (C) measurements are shown here performed in PPMS for 2 K-30 K under the magnetic fields up to 14 T at Pulsed Field Facility in Los Alamos. For dc magnetization and ac magnetic susceptibility measurements, the sample with 1 mm² was aligned along a, b^*, c, and ({c}^{{prime} }) and mounted on the sample holder inside the compensation and pick-up coil, where ({c}^{{prime} })-axis is 6° from c-axis on the b^*–c plane. The electrical contacts were made by painting silver epoxy on the two parallel opposite surfaces for the capacitance measurement and Andeen-Hagerling AH-2700A capacitance bridge was employed to measure the capacitance.

Magnetic torque was measured with the 31 T resistive magnet at the National High Magnetic Field Laboratory in Tallahassee, Florida using a home-made torque magnetometer. The torque magnetometer was formed from two parallel BeCu plates. A tiny aligned crystal with the mass of 0.13 mg was attached using epoxy on the cantilever plates. The capacitance change was monitored while the field or angle was swept using an AH-2700A capacitance bridge. More details can be found in https://nationalmaglab.org/user-facilities/dc-field/measurement-techniques/torque-magnetometry-dc/.

We then converted the measured capacitance to the deflected angle (α) using the following equation

$$C(alpha )=frac{{epsilon }_{0}{epsilon }_{r}w}{tan alpha }ln frac{ltan alpha -{h}_{0}}{{h}_{0}},$$

(3)

where α is the defected angle that is proportional to the magnetic torque, ϵ₀ and ϵ_r are vacuum and relative permittivity, respectively, and w, l, and h₀ are width, length of the BeCu plate and the gap distance between two BeCu plates at zero magnetic field. We extracted the coefficients c_2n using the following equation

$${c}_{2n}=frac{2}{pi }{int_{{theta }_{1}}^{{theta }_{i}+pi }}frac{tau }{H}sin (2ntheta )dtheta ,$$

(4)

where θ_i is the initial angle, τ is the measured magnetic torque signal.

Theoretical phase diagrams

In order to account for quantum fluctuation effects, we employ the method of second-order energy corrections outlined in ref. 32. At each site i, the local quantization axis (denoted ({hat{z}}_{i}^{{prime} })) is along the local spin vector, such that (langle {S}_{i}^{z{prime} }rangle =+frac{1}{2}) for all sites. Fluctuations are then driven by terms acting like ({J}_{ij}^{–}{S}_{i}^{-{prime} }{S}_{j}^{-{prime} }). The second-order energy correction due to such fluctuations is:

$$Delta {E}^{(2)}=sum _{{rm{bonds}},ij}frac{{J}_{ij}^{++}{J}_{ij}^{–}| langle m| {S}_{i}^{-{prime} }{S}_{j}^{-{prime} }| gsrangle {| }^{2}}{{E}_{gs,{rm{classical}}}-{E}_{m,{rm{classical}}}}$$

(5)

where E_gs,classical is the classical energy of the spin configuration, and E_m,classical is the classical energy of the state with the two spins along the bond ij flipped. For the example of a Heisenberg honeycomb antiferromagnet, the classical ground state energy is E_gs,classical = −0.375 J per site, whereas the present approach provides E_gs,classical + ΔE⁽²⁾ = −0.563 J per site. The latter estimate (which includes local singlet fluctuations) is in much better agreement with E_gs = −0.544 J per site from series expansion⁴¹ and quantum Monte Carlo⁴². In order to capture the energetically competitive magnetic orders depicted in Fig. 1b, we employ a combination of 4-, 6-, and 8-site clusters of spins with appropriate boundary conditions, and numerically minimize E_classical + ΔE⁽²⁾ with respect to the spin orientations. We consider the commensurate spiral state with Q = [1/4, 0] as a proxy for the proposed incommensurate spiral with Q = [0.27, 0]²⁷. We find the Q = [1/4, 0] spiral to be the classical ground state for a range of parameters, but it is completely replaced by the colinear ferromagnetic and Q = [1/4, 0] DS phases once ΔE⁽²⁾ is included.

First principles calculations

To estimate the magnetic couplings, we consider an electronic model for the Co 3d-orbitals, which is a sum of, respectively, one- and two-particle terms: ({{mathcal{H}}}_{{rm{el}}}={{mathcal{H}}}_{1p}+{{mathcal{H}}}_{2p}). The one-particle terms include intersite hopping, intrasite crystal field, and spin-orbit coupling, ({{mathcal{H}}}_{1p}={{mathcal{H}}}_{hop}+{{mathcal{H}}}_{{rm{CF}}}+{{mathcal{H}}}_{{rm{SO}}}):

$${{mathcal{H}}}_{hop}=sum _{ijalpha beta sigma }{t}_{ij}^{alpha beta }{c}_{i,alpha ,sigma }^{dagger }{c}_{j,beta ,sigma }$$

(6)

$${{mathcal{H}}}_{{rm{CF}}}=sum _{ialpha beta sigma }{d}_{i}^{alpha beta }{c}_{i,alpha ,sigma }^{dagger }{c}_{i,beta ,sigma }$$

(7)

$${{mathcal{H}}}_{{rm{SO}}}=sum _{ialpha beta sigma {sigma }^{{prime} }}{lambda }_{{rm{Co}}}leftlangle {phi }_{i}^{alpha }(sigma )| {bf{L}}cdot {bf{S}}| {phi }_{i}^{beta }({sigma }^{{prime} })rightrangle {c}_{i,alpha ,sigma }^{dagger }{c}_{i,beta ,{sigma }^{{prime} }}$$

(8)

where ({c}_{i,alpha ,sigma }^{dagger }) creates an electron at Co site i, in d-orbital α, with spin σ. To parameterize the one-particle terms, we performed fully relativistic density functional theory (DFT) calculations at the GGA (PBE⁴³) level using FPLO⁴⁴, a 12 × 12 × 12 k-grid, and the structure from ref. 45. The resulting Kohn-Sham Bloch functions were projected onto Co d-orbitals to yield single-particle terms from the Kohn-Sham Fock matrix written in the Wannier basis⁴⁶. For BCAO, this approach yields accurate crystal-field excitation energies⁴⁷, confirming the validity of the computed single-particle parameters.

For the two-particle terms, we consider a combination of on-site and nearest neighbor intersite contributions: ({{mathcal{H}}}_{2p}={{mathcal{H}}}_{U}+{{mathcal{H}}}_{J}^{{rm{nn}}}), respectively. The on-site contributions are given by:

$${{mathcal{H}}}_{U}=sum _{ialpha beta delta gamma }sum _{sigma {sigma }^{{prime} }}{U}_{alpha beta gamma delta },{c}_{i,alpha ,sigma }^{dagger }{c}_{i,beta ,{sigma }^{{prime} }}^{dagger }{c}_{i,gamma ,{sigma }^{{prime} }}{c}_{i,delta ,sigma }$$

(9)

In the spherically symmetric approximation⁴⁸, U_αβγδ are parameterized by the Slater parameters ({F}_{0}^{dd},{F}_{2}^{dd},{F}_{4}^{dd}). For the (screened) on-site metal Coulomb interactions, we use ({F}_{2}^{dd}=11.55) eV, and ({F}_{4}^{dd}=7.68) eV, following the fitted parameters for LiCoO₂ from ref. 49. We consider ({U}_{t2g}={F}_{0}^{dd}+(4/49)({F}_{2}^{dd}+{F}_{4}^{dd})) within the realistic range of 5.0–5.5 eV.

The nearest neighbor intersite Hund’s coupling is given by:

$${{mathcal{H}}}_{J}^{{rm{nn}}}=sum _{ijalpha beta sigma {sigma }^{{prime} }}{J}_{H,ij}^{alpha beta },{c}_{i,alpha ,sigma }^{dagger }{c}_{j,beta ,{sigma }^{{prime} }}^{dagger }{c}_{i,alpha ,{sigma }^{{prime} }}{c}_{j,beta ,sigma }$$

(10)

A significant contribution to the intersite magnetic couplings comes from ligand “charge transfer” exchange terms^1,4, which account for perturbative processes where two holes meet on a single oxygen ligand. When downfolded into the d-orbital Wannier basis, these appear as intersite Hund’s coupling terms. They give an overall ferromagnetic contribution to the intersite magnetic couplings between the low-energy j_1/2 doublets, which may be rendered anisotropic (bond-independent XXZ anisotropy) by the crystal-field splitting of the t_2g orbitals. In a previous study of BCAO²⁶, some of the authors accounted for this term via an estimated correction added to the computed couplings. Here, we employ the improved approach utilized in^40,50,51 that captures these effects at the level of the electronic Hamiltonian. The d-orbital Wannier functions can be written:

$$leftvert {tilde{d}}_{i,alpha ,sigma }rightrangle ={phi }_{i,alpha }^{0}leftvert {d}_{i,alpha ,sigma }rightrangle +sum _{n,beta }{phi }_{i,alpha }^{n,a}leftvert {p}_{n,a,sigma }rightrangle$$

(11)

where i labels a given Co metal site, α is a d-orbital, n labels an oxygen ligand site, and a ∈ {x, y, z} is a p-orbital. The tilde on the left side of the equation indicates a Wannier function, whereas the absence of a tilde indicates a bare atomic orbital. Rotating the p-orbital Coulomb terms into the (tilde{d}) Wannier basis results in residual intersite Hund’s coupling given by:

$$begin{array}{rcl}{J}_{H,ij}^{alpha beta }&=&mathop{sum}limits_{n}mathop{sum}limits_{alpha beta }left({U}^{p}-{J}_{H}^{p}right){phi }_{i,alpha }^{n,a}{phi }_{i,alpha }^{n,b}{phi }_{j,beta }^{n,a}{phi }_{j,beta }^{n,b}\ &&+,{J}_{H}^{p},{phi }_{i,alpha }^{n,a}{phi }_{i,alpha }^{n,a}{phi }_{j,beta }^{n,b}{phi }_{j,beta }^{n,b}end{array}$$

(12)

where U^p and ({J}_{H}^{p}) are the atomic ligand Hubbard and Hund’s couplings, respectively. The computed magnetic couplings depend most strongly on the overall magnitude of ({J}_{ij}^{alpha beta }), rather than the specific orbital dependence, which provides for some inconsequential flexibility in the choice of U^p and ({J}_{H}^{p}). We take ({U}^{p}=2{J}_{H}^{p}), and consider ({J}_{H}^{p}) within the range of 0.4 to 0.6 eV. This choice corresponds to roughly 60% to 80% of the atomic value for oxygen⁵². The sum over oxygen sites in Eq. (12) is then evaluated from the ab-initio Wannier function coefficients for the bridging oxygen atoms for each bond.