Field-tunable toroidal moment in a chiral-lattice magnet

Ferrotoroidal order, which represents a spontaneous arrangement of toroidal moments, has recently been found in a few linear magnetoelectric materials. However, tuning toroidal moments in these materials is challenging. Here, we report switching between ferritoroidal and ferrotoroidal phases by a small magnetic field, in a chiral triangular-lattice magnet BaCoSiO4 with tri-spin vortices. Upon applying a magnetic field, we observe multi-stair metamagnetic transitions, characterized by equidistant steps in the net magnetic and toroidal moments. This highly unusual ferri-ferroic order appears to come as a result of an unusual hierarchy of frustrated isotropic exchange couplings revealed by first principle calculations, and the antisymmetric exchange interactions driven by the structural chirality. In contrast to the previously known toroidal materials identified via a linear magnetoelectric effect, BaCoSiO4 is a qualitatively new multiferroic with an unusual coupling between several different orders, and opens up new avenues for realizing easily tunable toroidal orders.

I n localized spin systems, a toroidal moment, which violates space inversion and time-reversal symmetry, can be generated by a head-to-tail arrangement of magnetic moments 1,2 . Such a magnetic vortex can have two senses, corresponding to the opposite directions of the toroidal moment. Ferrotoroidicity, that is, the uniform arrangement of toroidal moments has been actively discussed as the fourth primary ferroic order, in addition to ferromagnetism, ferroelectricity, and ferroelasticity [2][3][4][5][6][7][8][9] . Ferrotoroidal order is thus of great fundamental and technological interest in condensed matter physics and spintronics 3,6,10,11 . From the symmetry-allowed terms in the free energy expression, toroidal moments should lead to antisymmetric components in the linear magnetoelectric effect. Indeed, discovered so far ferrotoroidal systems have been limited to a couple of linear magnetoelectric materials 5,12,13 .
Ferrotoroidal phases in linear magnetoelectric materials have been studied using either polarized neutron analysis or secondharmonic generation technique 5,12,13 . The pioneering work on LiCoPO 4 using the latter method has revealed the presence of toroidal domains 5 . However, by the very nature of the toroidal ordering in this compound, a hysteretic poling of ferrotoroidic domains was only realized by crossed magnetic and electric fields 14 . It is worth noting that the observed signal in LiCoPO 4 results from a staggered arrangement of toroidal moments, i.e., it is a ferri-(not ferro)toroidal order, related to the collinear magnetic structure with bi-spin vortices 15 . Another case is pyroxene LiFeSi 2 O 6 , which also shows a finite off-diagonal magnetoelectric effect 13 . The control of toroidal domains in this compound was only possible by applying crossed magnetic and electric fields, as revealed by the polarized neutron analysis. This appears to be a common problem in these toroidal materials, which require simultaneous application of crossed magnetic and electric field, and is a consequence of symmetry constraints [12][13][14] . This restriction is directly related to their linear magnetoelectric effect, rendering these materials curious, but impractical.
A natural question to ask is whether it is possible at all to control toroidal moments using a single field. Our work presented below answers the question affirmatively, and clarifies the underpinning physical picture. The direct coupling between ferromagnetism and ferroelectricity in such chiral multiferroic materials as BaCoSiO 4 means that toroidal moments can be readily controlled using a single field, in a natural environment for the simultaneous existence of multi-dipole orders.
In a chiral structure, an object subjected to spatial inversion cannot be superimposed upon itself by any combination of rotations and translations 16,17 . Owing to this special symmetry, chirality has been found to be instrumental in stabilizing unusual magnetic orders such as multiferroicity 18,19 , skyrmionic order [20][21][22] , helicity 23,24 , and chiral magnetic soliton lattice 25 . Chirality combined with the magnetic frustration characteristic of antiferromagnetic interactions in an equilateral triangle tends to generate a 120°vortex-like configuration 26,27 , as shown in Fig. 1, giving rise to a nonzero toroidal moment breaking both the spatial inversion and time-reversal symmetry [26][27][28] . Depending on the sense of the in-plane spin rotations, the toroidal moment is either positive (+) or negative (−). Manipulating toroidal moments directly by a magnetic field becomes possible in such a chiral magnetic vortex, if an out-of-plane spin component is present and coupled with the in-plane spin texture through Dzyaloshinskii-Moriya (DM) interactions 29,30 .
Following this strategy, in this work, we find a ferritoroidal order in a unique vortex-like spin configuration in the chiral magnet BaCoSiO 4 . By applying a small magnetic field, the toroidal moments are uniformly aligned, thereby leading to a ferrito ferrotoroidal transition. This toroidal transition, as well as the simultaneously scalar ferri-to ferrochiral transition, is fully explained within a magnetic Hamiltonian accounting for the magnetic frustration and antisymmetric DM interactions. A key property of this Hamiltonian, as derived from first-principles calculations, is a rather special and not immediately obvious hierarchy of Heisenberg exchange parameters, which does not correlate with the length of the corresponding Co-Co bonds.

Results
Crystal structure and multi-stair metamagnetic transitions. The stuffed tridymite BaCoSiO 4 crystallizes in the polar space group P6 3 . Its crystal structure is therefore chiral and adopts only one enantiomorph 31 . Co atoms are tetrahedrally coordinated by oxygen with a large off-center distortion and the nearest Co atoms form spin trimers in the ab plane (Fig. S1). In this structure, magnetic interactions between Co 2+ ions are expected to be small due to long and indirect exchange paths through adjacent SiO 4 tetrahedra. The temperature-dependent magnetic susceptibility shows an anomaly at T N~3 .2 K (Fig. S5), reflecting a longrange antiferromagnetic magnetic order as discussed in the following. The Curie-Weiss law describes well the high-temperature (150 ≲ T ≲ 300 K) magnetic susceptibility, with negative Weiss temperatures θ ab CW ¼ À10ð2Þ K for H∥ab and θ c CW ¼ À26:2ð4Þ K for H∥c. The fitted Curie constants correspond to effective moments μ ab eff ¼ 4:6ð4Þμ B and μ c eff ¼ 4:4ð2Þμ B , consistent with the high-spin state of the Co 2+ cation with S = 3/2 and a nearly isotropic g-tensor of ≈2.3 (note that this deviates considerably from the nonrelativistic value g = 2, indicating sizeable spin-orbit effects). As shown in Fig. 2a, the magnetization hysteresis loop of BaCoSiO 4 at 2 K exhibits a sequence of metamagnetic transitions. Starting with a zero-field cooled sample, the first transition tõ 0.1μ B is observed at low fields (≤150 Oe) for H∥c, stemming from an alignment of weak ferromagnetic domains. After a slow and linear ramp, a second transition to~0.4μ B occurs at a critical field μ 0 H C~1 .2 T. This corresponds to a spin flip in one of the ferrotoroidal sublattices, as we will elaborate further below. Similar transitions occur with small hysteresis loops for reversed fields. Using a pulsed magnetic field, we measured the magnetization up to 60 T along different crystallographic directions (Fig. 2a, inset). A much higher field is needed to reach saturation with the magnetic field applied along the c-axis, which implies the presence of a considerable easy-plane anisotropy, formally not expected for Co 2+ in a tetrahedral environment, but consistent with g > 2. The slope changes around 4 and 7 T before the saturation suggest additional transitions of a completely different nature (Fig. S5).
Spin vortex revealed by neutron diffraction. To further characterize the magnetic ground state, we have measured the powder neutron diffraction of BaCoSiO 4 in a zero magnetic field. The diffraction pattern at 1.8 K shows a set of satellite reflections that can be indexed with a propagation vector k = (1/3, 1/3, 0) with respect to the crystallographic unit cell. The thermal evolution of the reflection (2/3 2/3 0) confirms the magnetic origin of the satellite reflections (Fig. 2b, inset). A power-law fit of the integrated intensity as a function of temperature gives a critical exponent β = 0.37 (1) and T N = 3.24(1) K, in agreement with the magnetization results. The symmetry analysis 32 and Rietveld refinement 33 based on the neutron data yield a complex magnetic structure where the in-plane components of magnetic moments form a vortex-like configuration (Fig. 2d). The structure is consistent with the magnetic space group P6 3 (No. 173.129) with a ffiffi ffi 3 p ffiffi ffi 3 p magnetic supercell. To test the reliability of the refined in-plane spin orientation, we evaluate the profile factor R p of the fit as a function of uniform rotation (ϕ) of spins within the ab plane. Figure 2b reveals clear minima at ϕ~0°(the structure in Fig. 2d) and ϕ~120°, indicating the global orientation of the spin structure with respect to the lattice is strongly constrained. The bulk magnetization data imply the existence of a weak ferromagnetic canting along the c-axis, which is allowed by the magnetic space group symmetry; however, it is too small to be unequivocally determined from current neutron data. A satisfactory fit can be achieved by setting canting angles to zero, yielding an ordered moment m Co (0 T) = 2.71(5) μ B at 1.8 K and~3.67μ B at zero temperature from extrapolating the powerlaw fitting (Fig. S4). To appreciate the toroidal nature of the magnetic ground in BaCoSiO 4 , it is instructive to decompose the structure into three interpenetrating sublattices (red, cyan, and blue) (Fig. 2d), each of which is a network of trimers (up-and down-triangles). Spins on every trimer form a 120°configuration that resembles a vortex and generates a nonzero toroidal moment. All trimers belonging to a single sublattice have an identical toroidal moment, giving rise to three ferrotoroidal sublattices. In zero field, two of them (red and blue) have the same total moment t, while the remaining one (cyan) has the opposite moment, leading to a net moment of −1t or +1t within a macroscopic magnetic domain. We dub this structure ferritoroidal.
The field dependence of the magnetic ground state in BaCoSiO 4 is investigated using single-crystal neutron diffraction. Figure 2c shows the integrated intensity of the magnetic reflection (2/3 2/3 0) and the nuclear reflection (−1 1 0) as a function of the magnetic field along the c-axis at 1.5 K. The former is increasingly suppressed by fields and eventually disappears at μ 0 H C~1 .2 T, while the latter gains significant extra intensity, indicating a fieldinduced transition to a k = 0 magnetic structure. The refined k = 0 magnetic structure has the same magnetic space group symmetry as the zero-field structure. The key difference is that the sublattice (cyan) with the opposite toroidal moment is flipped by 180°in the field, yielding a uniformly aligned toroidal moment for all three sublattices with a total toroidal moment +3t (Fig. 2e), termed ferrotoroidal. This remarkable field-induced ferri-to ferrotoroidal transition directly correlates with the metamagnetic transition observed in the bulk magnetization measurements at the same critical field.
Density functional theory calculations. While the magnetic configuration, at first glance, seems nearly incomprehensibly complicated, it actually has a straightforward microscopic explanation. To show that, we first calculated the isotropic Heisenberg magnetic interactions using the density functional theory (DFT). Inspecting the Co t 2g bands crossing the Fermi level, we found that the five shortest Co-Co bonds, with d Co-Co between 5.11 and 5.41 Å, all have hopping integrals of roughly the same order (Fig. S6), so they should be included in the minimal model (Fig. 3a). Next, we performed total energy calculations for selected spin configurations and determined exchange parameters J by fitting the results to the mean-field energies of a Heisenberg model. Figure 3b shows the fitted model parameters as a function of onsite interaction U, using the room-temperature crystal structure as input in DFT calculations. All five exchange couplings are antiferromagnetic. Most interestingly, despite all bond lengths being similar, two interactions stand out as dominant, independent of U: the intralayer coupling J t and the interlayer coupling J z . The Co atoms connected by these two bonds form the three interpenetrating sublattices shown in Fig. 2d. Within each sublattice, we expect that frustrated J t interactions impose 120°s pin configurations on trimers that are antiferromagnetically coupled by J z . This explains the major part of the experimentally determined magnetic structure. In addition, relativistic DFT calculations indicate a strong easy-plane single-ion anisotropy (~2 K, comparable to the dominant exchange interactions). This is nontrivial, since Co 2+ in a tetrahedral environment features a full e g shell and half-filled t 2g , and formally is not supposed to have a sizeable orbital moment. The key is that the CoO 4 tetrahedra are considerably distorted and moreover Co 2+ is strongly off-centered (Fig. 3a), so that the e g and t 2g orbitals are actually mixed. This is also consistent with the pulsed-field magnetization measurements. Three relevant physical quantities are defined with spins {S i = 1, 2, 3} numbered anticlockwise: toroidal moment t = ∑ i r i × S i , where r i is the vector from the center of the triangle to spin S i ; vector chirality ϵ = S 1 × S 2 + S 2 × S 3 + S 3 × S 1 ; scalar chirality κ = (S 1 × S 2 ) ⋅ S 3 . Green symbols + and − for a toroidal moment and vector chirality denote the direction of these quantities with respect to the net magnetic moment, + for parallel and − for antiparallel. The magnetic vector chirality characterizes the sense of spin rotation along an oriented loop (or line), while the toroidal moment is associated with that around a center. Scalar spin chirality is a measure of non-coplanarity that does not necessarily have a sense of rotation. In the current example, toroidal moment and scalar chirality have a one-to-one correspondence to the net magnetic moment for a given handedness, since all of them are odd under time reversal. In a crystalline material with these chiral spin trimers as basic units, a ferro alignment of net magnetic moments of trimers will lead to a ferro ordering of toroidal moments and scalar chirality. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-25657-6 ARTICLE Of the three sublattices in Fig. 2d, one (cyan) has its toroidal moment opposite to the others. This results from the subleading exchange interactions fJ 0 t ; J 00 t ; J c g (Fig. 3a). The first two connect the trimers within the ab plane, while the last one connects those along the c-axis. A close inspection of the lattice connectivity reveals that all three subleading interactions with the help of J t form various triangular units. The total energy associated with subleading interactions is the lowest if spins on all these triangular units have 120°arrangements (neglecting the weak ferromagnetic canting). Yet, this condition cannot be satisfied. Figure 3d colors all "frustrated" triangles that do not have 120°c onfigurations in a ferri-and ferrotoroidal state. We see clear switching of colored triangles from one state to another, indicating a direct competition between the intralayer J 0 t and J 00 t couplings and the interlayer J c couplings. In BaCoSiO 4 , this competition favors a ferritoroidal state by a small margin of energy in zero field.
Antisymmetric DM interactions. However, this Heisenberg model, with (or without) the single-ion anisotropy, is insufficient in explaining the weak ferromagnetic canting and the spin-space anisotropy evidenced from our experimental data. To this end, we introduce the antisymmetric DM interactions 28-30 within the structural trimers (Fig. 3a). All three components of a DM vector on a nearest-neighbor bond are allowed due to the lack of symmetry constraints. It is rather difficult to calculate the DMI from the first principles, so for the moment we are assuming that all three components are present. The DM vectors on different bonds of the trimer are related by the 3-fold rotations. Assuming a uniform canting along c and 120°configurations in the ab plane for a trimer, the out-of-plane DM component D zẑ contributes À 3 ffiffi 3 p 2 jD z jS 2 xy in energy, where S xy is the length of in-plane spin component. Therefore, this term always favors coplanar spin configurations instead of canting. Taking into account that the 120°configuration has a 3-fold vortex symmetry, we find that the energy associated with the in-plane DM component, D xy , is 3 ffiffi ffi 3 p ðD xy Á S xy Þ S z Á b z À Á , where S z is the out-of-plane spin component. There are two important ramifications. First, the two components of spins are locked together, namely if S z flips, so does S xy , to keep the DM energy gain. This is actually the essential physical factor that couples the net component of magnetization p supercell solved from powder neutron diffraction data, showing three interpenetrating ferrotoroidal sublattices (red, blue, and cyan) formed by the dominant exchange interactions J t (intralayer) and J z (interlayer). The direction of the toroidal moment for each sublattice is denoted + if it is parallel to the c-axis and − if antiparallel. The red and blue sublattices have the same toroidal moment, while the cyan has the opposite moment, leading to a ferritoroidal state with a total moment +1t. The primitive crystallographic unit cell is indicated by the dotted lines. e Magnetic structure of BaCoSiO 4 at 2 T solved from single-crystal neutron diffraction data. All spins on the cyan sublattice are reversed, leading to a ferrotoroidal state with a total toroidal moment +3t. Triangles in panels (d, e) lie in two adjacent layers, which are bridged by the interlayer interaction J z . M z to the toroidal moment and allows to control toroidal moments by altering M z , e.g., by external fields. Second, to minimize the DM energy for a fixed spin canting, S xy has to be antiparallel to D xy if S z k b z or parallel to D xy if S z k Àb z, which means that we can read off the direction of the DM vector directly from the experimental spin structure, assuming D z = 0. In Fig. 3c, we show by green arrows the total DM vector for each bond, which is the vector sum of D ⊥ and D ∥ components and makes ã 30°angle with the bond, as determined from the magnetic structure. For each sublattice connected by J t and J z bonds, the DM interaction creates a uniform c-axis canting and corresponding in-plane toroidal spin texture. However, different sublattices are independent of each other, hence ferri-and ferrotoroidal states would have had the same energy, if not for the subleading Heisenberg interactions fJ 0 t ; J 00 t ; J c g.

Discussion
We now fully understand the metamagnetic ferrotoroidal transition. Indeed, because of the DMI-induced ferromagnetic canting, the small fJ 0 t ; J 00 t ; J c g-driven energy gain associated with the ferritoroidal arrangement competes directly with the Zeeman interaction favoring the ferrotoroidal phase. This leads to the spin flip (also a toroidal flip) manifested through the sudden increase of magnetization along the c-axis by a factor of three at the metamagnetic transition (Fig. 2e). A direct numerical calculation of magnetization using the complete model with all interactions discussed thus far is given in Fig. 2a and shows good agreement with the data. See Supplementary information and "Method" section for more details and model parameters.
At a fundamental level, the physics emerging in BaCoSiO 4 originates from its chiral crystal structure. For a single triangle in three dimensions, there is a set of three mirror planes perpendicular to the triangle and one extra mirror plane containing the triangle. The former allows both D z and D ⊥ components of the DM interaction, while the latter allows only D z . When D ⊥ is absent, D z with the correct sign could create a vortex configuration; however, the ferromagnetic canting is unfavored in this case. Thus, the intimate coupling between magnetization and toroidal moment is lost. In BaCoSiO 4 , the triangles are decorated by distorted CoO 4 tetrahedra that break all the mirror planes while still preserving the 3-fold symmetry. The D ∥ term is then allowed and plays an essential role in generating chiral vortex structures wherein the c-axis canting is coupled with the sense of the spin rotation in the ab plane. A direct control of toroidal moments using only magnetic fields is therefore possible through controlling the bulk magnetization by a uniform magnetic field, instead of using a conjugate field such as the curl of magnetic fields. The same arguments show that at this transition the scalar chirality κ (Fig. 1), existing in Co triangles, also changes from the ferrichiral  (Table S2). c Minimal energy configurations for three spins {S i , i = 1, 2, 3} on a triangle with an antiferromagnetic Heisenberg interaction and an in-plane DM interaction. The DM vectors {D i , i = 1, 2, 3} (green arrows) are related by 3-fold rotation symmetry and have the same sense of rotation as the tilting of apical oxygens shown in panel (a). Each vector makes an~30°angle with the bond, see main text for details. For this set of DM vectors, the resulting spin structure (pink arrows) generates a toroidal moment t (black arrows) that is always parallel to the magnetization M. d Energy balance between the ferri-to ferrotoroidal state in magnetic fields. The "frustrated" triangular units that do not have 120°configurations (and cost more energy) are highlighted in pink for both states. The ferrotoroidal state has less colored triangles in the ab plane; therefore, it is energetically favored by the interactions fJ 0 t ; J 00 t g, similarly the ferritoroidal state is favored by the interaction J c . Competition between these subleading interactions results in the ferritoroidal structure with a lower energy in zero field. The transition to the ferrotoroidal state occurs when the energy difference is compensated by the Zeeman energy in magnetic fields.
to ferrochiral state, similar to the three-sublattice description for the toroidal transition. The ferrochiral state with a noncoplanar spin configuration acquires a considerable Berry curvature, which can lead to a variety of exotic physical phenomena such as topological magnon excitations 34,35 .
In summary, we studied a rare chiral triangular-lattice magnet BaCoSiO 4 through bulk magnetization measurements and neutron diffraction experiments. We uncovered a novel vortex-like spin texture and a magnetic field-induced ferri-to ferrotoroidal transition for the first time. Combining ab initio DFT calculations and theoretical modeling, we have derived the microscopic energy balance and were able to explain quantitatively the complex magnetic structure and the field-induced metamagnetic/toroidal phase transition, neither of which had been observed before in any compound. Our work shows that BaCoSiO 4 is an excellent platform to study field-tunable toroidal moments and to explore their interplay with the structural and magnetic chirality. Further studies on the magnetoelectric effects and dynamical responses of toroidal spin textures are liable to bring up further new physics and potential applications.

Methods
Sample preparation and characterization. A powder sample of BaCoSiO 4 was prepared by the direct solid-state reaction from stoichiometric mixtures of BaCO 3 , Co 3 O 4 , and SiO 2 powders all from Alfa Aesar (99.99%). The mixture was calcined at 900°C in the air for 12 h and then re-grounded, pelletized, and heated at 1200°C for 20 h and at 1250°C for 15 h with intermediate grindings to ensure a total reaction. Finally, the sample was rapidly quenched from 1250°C to room temperature to avoid the decomposition at intermediate temperature. The resulting powder sample is fine and bright blue in color. Large single crystals were grown using a laser-diode heated floating zone technique. The optimal growth conditions were growth speed of 2-4 mm/h, atmospheric air flow of 0.1 L/min and counterrotation of the feed and seed rods at 15 and 30 r.p.m., respectively. Single-crystal xray diffraction data were collected at 95 K using a Rigaku XtaLAB PRO diffractometer with the graphite monochromated Mo K α radiation (λ = 0.71073 Å) equipped with a HyPix-6000HE detector and an Oxford N-HeliX cryocooler. Peak indexing and integration were done using the Rigaku Oxford Diffraction CrysA-lisPro software. An empirical absorption correction was applied using the SCALE3 ABSPACK algorithm as implemented in CrysAlisPro. Structure refinement was done using the FullProf suite 33 .
Magnetization measurement. Temperature dependence of magnetization M(T) was measured under a field of 0.1 T in a commercial magnetic property measurement system (MPMS-XL7, Quantum Design). Magnetic hysteresis loops with the field along the a-and c-axis were measured at 2 K using the same setup. A 2.6 mg piece oriented BaCoSiO 4 single crystal was glued on a 3 mm × 3 mm piece of weighing paper by GE varnish first, and then they were mechanically fixed on a straw attaching to the MPMS sample rod. This mounting setup was tested without placing any sample and it shows a diamagnetic background signal around −10 −6 emu at ground temperature, whereas the BaCoSiO 4 sample shows >10 −3 emu magnetization. Thus, a significant influence from sample mounting can be ruled out. The core diamagnetism has a typical magnitude of 10 −6 emu/mole, which is also negligible, compared with the observed susceptibility signal magnitude. Magnetization up to 60 T was measured by an induction magnetometry technique 36 using a capacitor-bank-driven pulsed magnet at the National High Magnetic Field Laboratory pulsed-field facility at Los Alamos. The pulsed-field magnetization values are calibrated against measurements in a 7 T direct current magnet (MPMS-XL7, Quantum Design).
Neutron diffraction. Neutron powder diffraction experiments were performed on the time-of-flight powder diffractometer POWGEN at the Spallation Neutron Source at Oak Ridge National Laboratory (ORNL). A powder sample of~2 g was loaded in a vanadium cylinder can and measured in the temperature range of 1.8-10 K with neutron wavelength band centered at λ = 1.5 and 2.665 Å, covering the d-space range 0.5-9.0 and 1.1-15.4 Å, respectively. Single-crystal neutron diffraction experiments were carried out on the single-crystal neutron diffractometer HB-3A DEMAND equipped with a 2D detector at the High Flux Isotope Reactor, ORNL. The measurement used the neutron wavelength of 1.553 Å selected by a bent perfect Si-220 monochromator 37,38 . The single crystal (~0.2 g) was mounted in a vertical field superconducting cryomagnet with a magnetic field up to 5.5 T and measured over the temperature range of 1.5-10 K with a magnetic field applied along the c-axis. The data refinements were performed by the FullProf suite 33 .
Spin-polarized DFT calculations. Electronic structure calculations were performed using the full potential local orbital basis 39 and generalized gradient approximation exchange and correlation functional 40 . The crystal structure from ref. 31 was used. We correct for the strong electronic correlations on Co 3d orbitals using a DFT+U method 41 with a fixed value of the Hund's rule coupling J H = 0.84 eV as suggested in ref. 42 . The Heisenberg Hamiltonian parameters were determined by an energy mapping technique 43 .
Full spin model calculations. The experimental magnetization data shown in Fig. 2a was modeled using a full spin Hamiltonian, including exchange interactions, single-ion anisotropy, and external fields, H ¼ 1 2 ∑ ij J ij S i Á S j þ A∑ i ðS z i Þ 2 À gμ B μ 0 ∑ i H Á S i , with g = 2 and S = 3/2. Starting with the Heisenberg interactions strengths produced by DFT calculations, by trial and error, we found the following set of parameters reasonably reproducing the experimental magnetization data, J t = 2.22 K, J z = 1.87 K, J 0 t ¼ 0:21 K, J 00 t ¼ 0:61 K, J c = 0.48 K, D ∥ = 1.71 K, D ? ¼ D k = ffiffi ffi 3 p , D z = 0 K, and A = 7.5 K. Direct numerical minimization of the classical energy was performed for this model using a ffiffi ffi 3 p ffiffi ffi 3 p supercell. The magnetization is obtained by averaging all spins of the lowest energy configuration at each field value, shown as the red line in Fig. 2a.

Data availability
The data that support the findings of this work are included in the article and the Supplementary information file. All raw data in the current work are available from the corresponding authors on reasonable request.