Introduction

In 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 moments1,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 ferroelasticity2,3,4,5,6,7,8,9. Ferrotoroidal order is thus of great fundamental and technological interest in condensed matter physics and spintronics3,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 materials5,12,13.

Ferrotoroidal phases in linear magnetoelectric materials have been studied using either polarized neutron analysis or second-harmonic generation technique5,12,13. The pioneering work on LiCoPO4 using the latter method has revealed the presence of toroidal domains5. 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 fields14. It is worth noting that the observed signal in LiCoPO4 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 vortices15. Another case is pyroxene LiFeSi2O6, which also shows a finite off-diagonal magnetoelectric effect13. 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 constraints12,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 BaCoSiO4 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 translations16,17. Owing to this special symmetry, chirality has been found to be instrumental in stabilizing unusual magnetic orders such as multiferroicity18,19, skyrmionic order20,21,22, helicity23,24, and chiral magnetic soliton lattice25. Chirality combined with the magnetic frustration characteristic of antiferromagnetic interactions in an equilateral triangle tends to generate a 120° vortex-like configuration26,27, as shown in Fig. 1, giving rise to a nonzero toroidal moment breaking both the spatial inversion and time-reversal symmetry26,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) interactions29,30.

Fig. 1: Toroidal moment and magnetic chirality of spin vortex configurations with 3-fold rotation symmetry.
figure 1

Four chiral noncoplanar structures are divided into groups of two left-handed ones and two right-handed ones. Those within each group are connected by 2-fold rotations and those between two groups by a mirror operation. The spins (axial vector) are represented as arrows. The dashed line with an ellipsoid indicates one of the three 2-fold axes. Three relevant physical quantities are defined with spins {Si = 1, 2, 3} numbered anticlockwise: toroidal moment t = ∑iri × Si, where ri is the vector from the center of the triangle to spin Si; vector chirality ϵ = S1 × S2 + S2 × S3 + S3 × S1; scalar chirality κ = (S1 × S2) ⋅ S3. 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.

Full size image

Following this strategy, in this work, we find a ferritoroidal order in a unique vortex-like spin configuration in the chiral magnet BaCoSiO4. By applying a small magnetic field, the toroidal moments are uniformly aligned, thereby leading to a ferri- to 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 BaCoSiO4 crystallizes in the polar space group P63. Its crystal structure is therefore chiral and adopts only one enantiomorph31. 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 Co2+ ions are expected to be small due to long and indirect exchange paths through adjacent SiO4 tetrahedra. The temperature-dependent magnetic susceptibility shows an anomaly at TN ~ 3.2 K (Fig. S5), reflecting a long-range 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 \({\theta }_{{{{{{{{\rm{CW}}}}}}}}}^{ab}\,=-10(2)\,{{{{{\mathrm{K}}}}}}\) for H∥ab and \({\theta }_{{{{{{{{\rm{CW}}}}}}}}}^{c}\,=-26.2(4)\,{{{{{\mathrm{K}}}}}}\) for H∥c. The fitted Curie constants correspond to effective moments \({\mu }_{{{{{{{{\rm{eff}}}}}}}}}^{ab}=4.6(4){\mu }_{{{{{{{{\rm{B}}}}}}}}}\) and \({\mu }_{{{{{{{{\rm{eff}}}}}}}}}^{c}=4.4(2){\mu }_{{{{{{{{\rm{B}}}}}}}}}\), consistent with the high-spin state of the Co2+ 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 BaCoSiO4 at 2 K exhibits a sequence of metamagnetic transitions. Starting with a zero-field cooled sample, the first transition to ~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 μ0HC ~ 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 Co2+ 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).

Fig. 2: Zero-field magnetic structure and field-induced ferri- to ferrotoroidal transition in BaCoSiO4.
figure 2

a The isothermal bulk magnetization data (open circles) measured at 2 K with the field parallel to the c-axis and the theoretical calculations (red line) from the full spin Hamiltonian, showing a good agreement between data and calculations. The magnetization up to 60 T in a pulsed magnetic field at 1.5 K is shown in the lower inset, and the magnetization hysteresis at low fields in the upper inset. See caption of panel (d) for the meaning of symbols + and −. b The refined agreement factor for the powder neutron diffraction data as a function of uniform rotation in the ab plane. The dashed line marks the best refinement. Inset shows the integrated intensity of the magnetic reflection (2/3 2/3 0) as a function of temperature with an order parameter fit \(\sim {(1-T/{T}_{{{{{{{{\rm{N}}}}}}}}})}^{2\beta }\) (solid line), where β is the critical exponent. The error bars are used to show the standard deviation given by the square root of the number of neutron counts. c The integrated intensities of reflection (2/3 2/3 0) and (−1 1 0) as a function of the magnetic field with H∥c at 1.5 K. d Zero-field magnetic structure of BaCoSiO4 in a \(\sqrt{3}\,\times \sqrt{3}\) supercell solved from powder neutron diffraction data, showing three interpenetrating ferrotoroidal sublattices (red, blue, and cyan) formed by the dominant exchange interactions Jt (intralayer) and Jz (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 BaCoSiO4 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 Jz.

Full size image

Spin vortex revealed by neutron diffraction

To further characterize the magnetic ground state, we have measured the powder neutron diffraction of BaCoSiO4 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 TN = 3.24(1) K, in agreement with the magnetization results. The symmetry analysis32 and Rietveld refinement33 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 P63 (No. 173.129) with a \(\sqrt{3}\,\times \sqrt{3}\) magnetic supercell. To test the reliability of the refined in-plane spin orientation, we evaluate the profile factor Rp 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 mCo(0 T) = 2.71(5) μB at 1.8 K and ~3.67μB at zero temperature from extrapolating the power-law fitting (Fig. S4). To appreciate the toroidal nature of the magnetic ground in BaCoSiO4, 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 BaCoSiO4 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 μ0HC ~ 1.2 T, while the latter gains significant extra intensity, indicating a field-induced 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 t2g bands crossing the Fermi level, we found that the five shortest Co–Co bonds, with dCo–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 Jt and the interlayer coupling Jz. The Co atoms connected by these two bonds form the three interpenetrating sublattices shown in Fig. 2d. Within each sublattice, we expect that frustrated Jt interactions impose 120° spin configurations on trimers that are antiferromagnetically coupled by Jz. 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 Co2+ in a tetrahedral environment features a full eg shell and half-filled t2g, and formally is not supposed to have a sizeable orbital moment. The key is that the CoO4 tetrahedra are considerably distorted and moreover Co2+ is strongly off-centered (Fig. 3a), so that the eg and t2g orbitals are actually mixed. This is also consistent with the pulsed-field magnetization measurements.

Fig. 3: Microscopic magnetic model and the underlying mechanism for the ferritoroidal to ferrotoroidal transition.
figure 3

a Magnetic exchange pathways of BaCoSiO4, showing three intralayer couplings \(\{{J}_{t},{J}_{t}^{\prime},{J}_{t}^{^{\prime\prime} }\}\) and two interlayer couplings {Jz, Jc}. Three components of a DM vector on the nearest-neighbor Jt bonds are indicated by red arrows in a local reference frame. b Density functional theory calculation of exchange interaction strengths as a function of onsite interaction U. The dashed line marks the set of couplings with U = 4.41 eV that matches the Weiss temperature from the magnetic susceptibility (Table S2). c Minimal energy configurations for three spins {Si, i = 1, 2, 3} on a triangle with an antiferromagnetic Heisenberg interaction and an in-plane DM interaction. The DM vectors {Di, 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 \(\{{J}_{t}^{\prime},{J}_{t}^{^{\prime\prime} }\}\), similarly the ferritoroidal state is favored by the interaction Jc. 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.

Full size image

Of the three sublattices in Fig. 2d, one (cyan) has its toroidal moment opposite to the others. This results from the subleading exchange interactions \(\{{J}_{t}^{\prime},{J}_{t}^{^{\prime\prime} },{J}_{c}\}\) (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 Jt 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° configurations 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}_{t}^{\prime}\) and \({J}_{t}^{^{\prime\prime} }\) couplings and the interlayer Jc couplings. In BaCoSiO4, 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 interactions28,29,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}\hat{{{{{{{{\bf{z}}}}}}}}}\) contributes \(-\!\frac{3\sqrt{3}}{2}| {D}_{z}| {S}_{xy}^{2}\) in energy, where Sxy 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, Dxy, is \(3\sqrt{3}({{{{{{{{\bf{D}}}}}}}}}_{xy}\cdot {{{{{{{{\bf{S}}}}}}}}}_{xy})\left({{{{{{{{\bf{S}}}}}}}}}_{z}\cdot \widehat{{{{{{{{\bf{z}}}}}}}}}\right)\), where Sz is the out-of-plane spin component. There are two important ramifications. First, the two components of spins are locked together, namely if Sz flips, so does Sxy, to keep the DM energy gain. This is actually the essential physical factor that couples the net component of magnetization Mz to the toroidal moment and allows to control toroidal moments by altering Mz, e.g., by external fields. Second, to minimize the DM energy for a fixed spin canting, Sxy has to be antiparallel to Dxy if \({{{{{{{{\bf{S}}}}}}}}}_{z}\parallel \widehat{{{{{{{{\bf{z}}}}}}}}}\) or parallel to Dxy if \({{{{{{{{\bf{S}}}}}}}}}_{z}\parallel -\widehat{{{{{{{{\bf{z}}}}}}}}}\), which means that we can read off the direction of the DM vector directly from the experimental spin structure, assuming Dz = 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 a ~30° angle with the bond, as determined from the magnetic structure. For each sublattice connected by Jt and Jz 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 \(\{{J}_{t}^{\prime},{J}_{t}^{^{\prime\prime} },{J}_{c}\}\).

Discussion

We now fully understand the metamagnetic ferrotoroidal transition. Indeed, because of the DMI-induced ferromagnetic canting, the small \(\{{J}_{t}^{\prime},{J}_{t}^{^{\prime\prime} },{J}_{c}\}\)-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 BaCoSiO4 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 Dz and D⊥ components of the DM interaction, while the latter allows only Dz. When D⊥ is absent, Dz 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 BaCoSiO4, the triangles are decorated by distorted CoO4 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 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 excitations34,35.

In summary, we studied a rare chiral triangular-lattice magnet BaCoSiO4 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 BaCoSiO4 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 BaCoSiO4 was prepared by the direct solid-state reaction from stoichiometric mixtures of BaCO3, Co3O4, and SiO2 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 counter-rotation of the feed and seed rods at 15 and 30 r.p.m., respectively. Single-crystal x-ray 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 CrysAlisPro software. An empirical absorption correction was applied using the SCALE3 ABSPACK algorithm as implemented in CrysAlisPro. Structure refinement was done using the FullProf suite33.

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 BaCoSiO4 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 BaCoSiO4 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 technique36 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 monochromator37,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 suite33.

Spin-polarized DFT calculations

Electronic structure calculations were performed using the full potential local orbital basis39 and generalized gradient approximation exchange and correlation functional40. The crystal structure from ref. 31 was used. We correct for the strong electronic correlations on Co 3d orbitals using a DFT+U method41 with a fixed value of the Hund’s rule coupling JH = 0.84 eV as suggested in ref. 42. The Heisenberg Hamiltonian parameters were determined by an energy mapping technique43.

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, \({{{{{{{\mathcal{H}}}}}}}}=\frac{1}{2}{\sum }_{ij}{J}_{ij}{{{{{{{{\bf{S}}}}}}}}}_{i}\cdot {{{{{{{{\bf{S}}}}}}}}}_{j}+A{\sum }_{i}{({S}_{i}^{z})}^{2}-g{\mu }_{{{\mbox{B}}}}{\mu }_{0}{\sum }_{i}{{{{{{{\bf{H}}}}}}}}\cdot {{{{{{{{\bf{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, Jt = 2.22 K, Jz = 1.87 K, \({J}_{t}^{\prime}=0.21\) K, \({J}_{t}^{^{\prime\prime} }=0.61\) K, Jc = 0.48 K, D∥ = 1.71 K, \({D}_{\perp }={D}_{\parallel }/\sqrt{3}\), Dz = 0 K, and A = 7.5 K. Direct numerical minimization of the classical energy was performed for this model using a \(\sqrt{3}\,\times \sqrt{3}\) 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.