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 (Co2+, 3d7) in a nearly ideal octahedral oxygen environment stabilizes the Jeff = 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 Jeff = 1/2 state can potentially exhibit bond-dependent exchange interactions, including what is lately called the Kitaev exchange interaction1,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 frustration5,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 computing8,9,10. Therefore, several honeycomb cobalt-based magnets, including Na2Co2TeO611,12,13, and Na3Co2SbO614 have been investigated intensively.
Recently, another cobalt-based honeycomb-lattice material BaCo2(AsO4)2 (BCAO) has been proposed as an ideal platform for realizing a KQSL15,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 field17,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 suggested15,21. For field in the out-of-plane direction, which is theoretically more conducive to a field-induced KQSL phase22,23,24,25, a broad magnetic continuum was observed16. 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-J1–J3 model with weakly bond-dependent couplings; recent inelastic neutron scattering measurements, as well as ab initio calculations26, align more with the XXZ-J1–J3 model27. 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 behavior4,28,29. Additionally, the weak ligand-assisted exchange interactions weaken bond-dependent exchange interactions4. 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 anisotropy30. Recently, Safari et al. explored the phase diagram of BCAO using magnetotropic susceptibility measurements and theoretical models31. Their findings indicate that while classical models qualitatively describe the system’s behavior, quantum corrections are essential for precise critical field predictions31.

a Crystal structure of BaCo2(AsO4)2 showing honeycomb magnetic lattice of Co in the ab-plane layers. Blue balls denote Cobalt, green balls Arsenic and red balls oxygen. We omitted Ba as they exist between layers. The a, b, and c are the global coordinate. The a and b point along two zigzag directions of the honeycomb lattice while a* and b* points along the armchair direction. The c points along the out-of-plane direction. The x, y and z are local coordinates used in Eq. (2), where the x-axis points along the bond-direction while z points along the out-of-plane c axis. The y is perpendicular to the x and the z directions. J1 and J3 are the nearest and third nearest interactions, respectively. b Energetically competitive magnetic orders at zero and finite field, with in-plane ordering wavevectors (Q) indicated.
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 studies15,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 corrections32 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 assertions26,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 (H1) and 0.2 T (H2), 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 (H3). The saturation magnetization is about 2.3 μB/Co2+, from which we extracted gab ~ 4.6. The magnetization between the saturation field and ~ 0.15 T constitutes a one-third magnetization plateau. These findings are consistent with prior studies15,19,20. When the temperature exceeds TN ~ 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 ~ H1 and ~ H2 are consistent with the jump in magnetization in the hysteresis loop and the end of the hysteresis loop. The third peak above ~ H3 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 loop17,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.

a Magnetization per Co2+ at various temperatures with field parallel with a. For the clarity, a 0.5 μB/Co2+ offset was applied between the two curves. b Ac susceptibility as a function of field at various temperatures with field parallel with a. c Magnetic field vs. temperature phase diagrams with field along the a-axis. d Magnetization per Co2+ at various temperatures with field parallel with b*. A 0.5 μB/Co2+ offset was applied between the two curves. e Ac susceptibility as a function of field at various temperatures with field parallel with b* direction. f Magnetic field vs. temperature phase diagrams with field along the b*-axis. Data points are taken from dc magnetization M, ac susceptibility ({chi }^{{prime} }), and capacitance measurement C (see Supplementary Fig. 2) for the phase diagrams.
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 H3 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 (Mc/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 TN 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.

a Temperature dependent DC susceptibility at various magnetic field along out-of-plane direction. An offset of 0.001 μB/Co2+/T multiplied by the field is applied to each curve for clarity. b Field dependent magnetization per Co2+ measured up to 14 T at various different temperatures, with curves offset by 0.05 μB/Co2+. c Pulsed field magnetization up to 50 T at different temperatures. d Differentiated magnetization with respect to the field.
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 TN, 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 dMc/dH is a constant below ~ 20 T and above ~ 35 T. Due to the small crystal field excitation gap of Co2+, it is common to observe a linearly increasing magnetization above the saturation field, due to Van Vleck paramagnetism36,37. By extrapolating the Van Vleck contribution to the zero field, we found a saturation moment of about 0.71 μB/Co2+, from which we extract gc ~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.16 and Safari et al.31. 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 TN 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 ~ H1 and ~ H2 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.

a Magnetization per Co2+ as a function of field at various angles at T = 2 K. 0° is c-axis and 90° is b*-axis. b Magnetic field dependent capacitance with electric field along b*-axis at different angles at T = 2.5 K. For clarity, we applied the offset of 1 × 10−5 pF between two curves. The arrows indicate the phase transition fields. c Field dependent magnetic torque at various angles at T = 0.5 K. d Differentiated magnetic torque with respect to the magnetic field at various angles with an offset. The arrows indicate the peak positions. All these measurements were performed during the upsweep of the magnetic field.
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 BaCo2(AsO4)2. 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 H1 and H3. These two features are only clearly observable up to 10 T and they merge into one above 10 T.

a b*c-plane magnetic torque as a function of angle for various magnetic fields measured up to 32 T at T = 0.5 K. b Theoretical simulation of torque on the basis of second-order corrected state energies employing optimal model discussed in the text. c Coefficients (leftvert {c}_{2}rightvert) and (leftvert {c}_{4}rightvert) as a function of magnetic field extracted from 0.5 K experimental data. d The same coefficients extracted from 0.0 K theoretical simulations. e Critical magnetic fields as a function of field angle in the b*c-plane. Data points were obtained through the analysis of the magnetization, capacitance, and torque measurements. The red and black symbols correspond to the H1 and H3, respectively. The dashed lines are guides for the eye. f Theoretical phase diagram for optimized model (see text).
In general, the magnetic torque may be expanded as30:
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 c2 to c4 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 c2n coefficients have two contributions. The first contribution is field-independent, and depends exclusively on the g-anisotropy. For the implied range of g-values (gab ~4.6, gc = 1.5−2.5), the limiting ratio for this contribution leads to ∣c4/c2∣ = 0.17−0.27. The next order contribution to the c2n scales as 1/H, and is related to both the XXZ exchange anisotropy and g-anisotropy30. The findings that c2 and c4 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 (H1), and the black empty symbols indicate the saturation field (H3). 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) approach38,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 j1/2 states. This yields gab = 5.1 and gc = 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 j1/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:
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 CoO6 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 eg 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 Co2+ materials can, in principle, host strongly bond-dependent couplings1,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 Jeff = 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 eg orbitals with the oxygen p-orbitals. Unlike the couplings between t2g 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 BCAO26,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 anisotropy26,31.
At zero-field, it is instructive to compare the classical and second order corrected quantum energies (denoted ECl and EQM, 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, EQM 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 EQM. 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 ⋅ Si∣ 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: H1 = 0.16 T, H3 = 0.55 T, and for H ∣∣ b*: H1 = 0.18 T, H3 = 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.

a, b Comparison of 2nd order corrected energies per site for different magnetic states for fields H∣∣a and H∣∣b*, respectively. c, d Theoretical magnetization per site for different field directions showing 1/3 magnetization plateau. Solid lines indicate the magnetic structures depicted in Fig. 1(b), while the dashed lines represent a 60°-rotated version of these magnetic structures.
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 H3 = 36 T, while for θ = 1°, we compute H3 = 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 c2n. In Fig. 5d, we show the theoretical evolution of these coefficients with field. The finding that c2 and c4 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 anisotropy26,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 of31. 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 with33,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 J1–J3 model with J1 < 0, J3 > 0 and significant easy-plane bond-independent XXZ anisotropy (see e.g., 27). 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 BaCo2(AsO4)2 synthesized through a flux method at Oak Ridge National Laboratory in Oak Ridge, for which the detailed methods are described elsewhere27. At low temperature, dc magnetization (M) is measured up to 60 T through millisecond timescales in pulsed magnet using custom-built 3He 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 mm2 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
where α is the defected angle that is proportional to the magnetic torque, ϵ0 and ϵr are vacuum and relative permittivity, respectively, and w, l, and h0 are width, length of the BeCu plate and the gap distance between two BeCu plates at zero magnetic field. We extracted the coefficients c2n using the following equation
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:
where Egs,classical is the classical energy of the spin configuration, and Em,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 Egs,classical = −0.375 J per site, whereas the present approach provides Egs,classical + ΔE(2) = −0.563 J per site. The latter estimate (which includes local singlet fluctuations) is in much better agreement with Egs = −0.544 J per site from series expansion41 and quantum Monte Carlo42. 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 Eclassical + ΔE(2) 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]27. 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(2) 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}}}):
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 (PBE43) level using FPLO44, 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 basis46. For BCAO, this approach yields accurate crystal-field excitation energies47, 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:
In the spherically symmetric approximation48, 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 LiCoO2 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:
A significant contribution to the intersite magnetic couplings comes from ligand “charge transfer” exchange terms1,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 j1/2 doublets, which may be rendered anisotropic (bond-independent XXZ anisotropy) by the crystal-field splitting of the t2g orbitals. In a previous study of BCAO26, some of the authors accounted for this term via an estimated correction added to the computed couplings. Here, we employ the improved approach utilized in40,50,51 that captures these effects at the level of the electronic Hamiltonian. The d-orbital Wannier functions can be written:
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:
where Up 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 Up 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 oxygen52. 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.
Responses