Introduction

In recent years, there has been growing interest in materials exhibiting geometric frustration due to their potential to support a quantum spin liquid (QSL) state1,2. A quantum spin liquid is a unique state of matter where the spins remain highly disordered and fluctuate even at very low temperatures. This frustration commonly occurs in materials where magnetic rare-earth ions are arranged on a two-dimensional triangular lattice. Of these, YbMgGaO4 has received considerable attention as being a potential host of a quantum spin liquid state3,4,5,6,7. However, YbMgGaO4 has random site mixing between the Mg and Ga atoms, which can make interpretation of the measurements difficult, and has been the source of much debate about the true nature of the ground state8,9,10,11. This has led to further study of other ytterbium-based triangular materials without site disorder, such as the rare-earth chalcogenides AYbX2 (A = Na, K, Rb, Cs; X = O, S, Se)12,13,14,15,16,17,18,19,20 or the orthoborate materials (Na, K)(Sr, Ba)Yb(BO3)221,22,23,24,25. We have recently reported the single crystal study of another triangular lattice borate compound without site mixing: Ba3Yb(BO3)3, and have found it to exhibit a spin-\(\frac{1}{2}\) quantum dipole ground state26.

Due to the chemical similarity of the rare-earth elements, different rare-earth elements can be substituted for Yb in the lattice, changing the spin of the magnetic ion without significantly altering the crystal structure. Many of the extensively studied Yb-based triangular materials have also been reported with other rare-earths substituted into the lattice. Some examples are ErMgGaO427 and TmMgGaO428, the rare-earth chalcogenides ARECh2 (A = Na, K, Rb; RE = rare-earth; Ch = O, S, Se)17,29,30,31,32,33,34 and several examples in the orthoborate family (Na, K)(Sr, Ba)RE(BO3)2 (RE = rare-earth)22,25,35,36,37.

A natural next step in the study of these materials is substituting spin-\(\frac{1}{2}\) Yb with spin-\(\frac{3}{2}\) Er. Since Er3+ is also a Kramers ion, it is expected to have a doublet ground state like Yb, so these compounds would provide an opportunity to study the effect of higher spin on the ground states. Both powder and single crystal samples of ErMgGaO4 have been produced27 which show similar site mixing to YbMgGaO4. Muon spin resonance (μSR) measurements show no signs of ordering down to 25 mK in ErMgGaO4 and rule out the presence of a glassy state38. Inelastic neutron scattering studies have been performed on NaErS2 and (K, Cs)ErSe2 to extract the crystal electric field (CEF) parameters; these materials show a large \({J}_{z}=\pm \!\frac{1}{2}\) component to the ground state, which allows for quantum transitions between the doublet states39,40. KErSe2 was also found to order into a stripe antiferromagnetic state at low temperatures41 and to host low energy spin-wave-like excitations42. Motivated by these reported results and our recent work on Ba3Yb(BO3)326, we have begun studies of the analogous compound Ba3Er(BO3)3.

In this study, we detail our efforts in synthesizing both powder and single crystal forms of Ba3Er(BO3)3. We conducted comprehensive magnetic and thermal characterization, complemented by inelastic neutron scattering measurements on this material. These measurements allowed us to derive the CEF spectrum and propose a potential magnetic ground state for Ba3Er(BO3)3. Our examination of single crystals at low temperatures uncovered a potential phase transition occurring at approximately 100 millikelvin. Through symmetry analysis and classical simulation, we propose that Ba3Er(BO3)3 exhibits two-sublattice exchange interactions: In this scenario, the spins on a honeycomb lattice exhibit ferromagnetic tendencies due to the presence of additional spins at the hexagonal centers, ultimately leading to an ordered antiferromagnetic phase.

Results and discussion

Crystal structure and symmetry considerations

The double borates form two different crystal structures—trigonal with space group R\(\bar{3}\) or hexagonal with space group P63cm—depending on which rare-earth is used43. Ba3Er(BO3)3 forms the hexagonal structure with space group P63cm (space group no. 185)44,45,46. The crystal structure is illustrated in Fig. 1a, b. The phase purity of Ba3Er(BO3)3 was confirmed using powder x-ray diffraction. Rietveld refinement to the P63cm structure is shown in Supplementary Fig. 1. The corresponding point group is C6v with 12 elements, generated by a three-fold rotation C3, a two-fold screw 21 along z, and a mirror M. The lattice structure and symmetries are shown in Fig. 1c. The unit cell is spanned by the lattice vectors a, b, and c (see Fig. 1) and contains six Er3+ ions (two layers × three Er3+ ions in each layer). These Er3+ ions occupy two distinct Wyckoff positions and have different symmetry properties, as detailed in Table 1. As we will show, the two environments are important in the understanding of magnetic properties of Ba3Er(BO3)3.

Fig. 1: Crystal structure, basis vectors and symmetry of Ba3Er(BO3)3.
figure 1

a Crystal structure of Ba3Er(BO3)3 viewed along the c-axis, showing the superexchange pathways between the Er3+ ions. The magenta and purple atoms are the 2a and 4b Er sites, large light green atoms are Ba, small red atoms are O, and small dark green atoms are B. The colored arcs indicate the different O-B-O angles in the two exchanges (black arcs for J2a and red/blue arcs for J4b. b Crystal structure viewed along the a-axis, showing the hexagonal layers of Er atoms, separated by layers of Ba, B, and O atoms). c Symmetry operations of the two Er sites. The symmetry operations 21, C3, and M correspond to a 180-degree screw axis in the z-direction, a three-fold rotation about the z-axis, and a mirror plane, respectively. The two colors indicate the two Wyckoff positions occupied by the Er3+ atoms in Ba3Er(BO3)3.

Full size image
Table 1 Wyckoff positions.
Full size table

Magnetic susceptibility and magnetization

Magnetic susceptibility measurements were performed on a single crystal sample with field applied both parallel to the c-axis and parallel to the ab-plane down to 300 mK using a helium-3 insert, as shown in Fig. 2a. As shown in the inset of Fig. 2a, in both orientations we observe no sharp features in the susceptibility down to 300 mK that would indicate magnetic ordering. The inverse susceptibility data have been fit to the Curie–Weiss law: \({\chi }^{-1}=\frac{T-\theta }{C}\), where θ is the Curie–Weiss temperature and C is the Curie constant. The inverse susceptibility data and Curie–Weiss fits are shown in Fig. 2b. We performed fits in the high-temperature region (100–300 K) and in the low-temperature region (10–30 K). The fit parameters were used to compute the effective magnetic moment (μeff) of the Er3+ ions. The Curie–Weiss fit in the high-temperature region was also used to calculate the Landé g-factor.

Fig. 2: Magnetic susceptibility data for Ba3Er(BO3)3.
figure 2

a Temperature-dependent magnetic susceptibility data for a single crystal sample of Ba3Er(BO3)3. Black points show data with field applied in-plane (χab) and blue points show data with field applied along the c-axis (χc). The inset shows the low-temperature region from 300 mK to 1.6 K, showing the anisotropy in the system. b Inverse magnetic susceptibility data are shown. Black points show data with field applied in plane, and blue points show data with field applied along the c-axis. Red lines show the fit to the Curie–Weiss law.

Full size image

In Ba3Er(BO3)3 we observe strong directional anisotropy in the susceptibility data. We do not observe any discrepancy between the field-cooled and zero field-cooled data in either orientation. In previously reported triangular compounds Ba3Yb(BO3)3, YbMgGaO4, and NaYbBa(BO3)2, large anisotropy is also observed4,21,26. In these compounds the intralayer Yb–Yb distance is less the interlayer distance, and the susceptibility and magnetization is larger with field applied along the c-axis, indicating Ising-like spins with an easy axis along the c-direction. In Ba3Er(BO3)3 the trend is reversed: the susceptibility χc is smaller than χab, particularly at low temperatures, as shown in Fig. 2a. This is despite Ba3Er(BO3)3 having nearly the same crystal structure as Ba3Yb(BO3)3, with the intralayer Er–Er distance being smaller than the interlayer Er–Er distance (see Table 2). A similar anisotropy (χab > χc) is reported the analogous Er-based system ErMgGaO4 at low temperatures27,38. The change in the nature of the anisotropy indicates that the Er spins take on XY-like behavior with an easy-plane perpendicular to the c-axis, even though there is no significant change in the crystal structure compared to the Yb compounds. This change in the anisotropic behavior does not occur in all layered triangular systems: in the chalcogenides ARSe2 (A = Na, K, Cs; R = Yb, Er) χab > χc regardless of the magnetic ion12,19,20,30,40,42, while for NaBaR(BO3)2 (R = Yb, Er) χc > χab for both Yb and Er21,36.

Table 2 Distances and exchange parameters.
Full size table

As shown in Fig. 2b, the inverse susceptibility curves for both orientations are nearly parallel from 100 K to 300 K, with θab = −3.75 K and θc = −15.46 K. These values are smaller than those observed in ErMgGaO427, but are similar to the values observed in the Er chalcogenides NaErCh2 (Ch = O, S, Se)17,29,30. The calculated effective moments and Landé g-factors in this region are μeff(ab) = 9.02 μB and gab = 1.13 for H∥ab and μeff(c) = 9.07 μB and gc = 1.13 for H∥c, in good agreement with the μ = 9.58 μB and g = 1.2 expected for a free Er3+ ion47.

With field applied in the ab-plane, the inverse susceptibility remains nearly linear down all the way to 300 mK. In contrast, the perpendicular direction is linear down to around 10 K, where it levels off slightly before turning down sharply at 5 K. This downturn behavior is also observed in Yb-based triangular lattice compounds, such as Ba3Yb(BO3)3, NaYbO2, and NaBaYb(BO3)316,17,21,26, and is believed to be caused by thermal population of CEF levels. The Curie–Weiss fits from 10 K to 30 K give θab = −0.1876 K and μeff(ab) = 8.51μB for H∥ab, showing a slight decrease in both values for both parameters, and θc = −24.16 K and μeff(c) = 9.50 μB for H∥c. Below 5 K the parameters sharply decrease. The inverse susceptibility remaining nearly linear down to very low temperatures is observed in several other Er-based triangular lattice systems22,27,29, and the discrepancy between the low temperature behavior for the two orientations is also seen in single crystals of NaErSe230.

To further investigate the magnetic behavior of this material, isothermal magnetization measurements were carried out at different temperatures ranging from 300 mK to 2 K with field up to 7 T applied parallel to the c-axis and parallel to the ab-plane, as shown in Fig. 3a. For both orientations there is a nonlinear region at low field where the magnetization rapidly increases before the curve becomes nearly linear above 3 T. With field applied in the ab-plane, the magnetization increases quickly at first and then slowly tapers off as the field increases. In contrast, with field applied along the c-axis the magnetization begins to level off around 1 T before increasing again. This behavior can be seen more clearly in the plot of dM/dH shown in Fig. 3b. The dM/dH curve goes down to a local minimum at 1 T, then rises to a local maximum at 3 T before going back down at higher fields. As with susceptibility, the magnetization data shows strong anisotropy in this material. Similar anisotropy in the saturated magnetization and curve shape are observed in ErMgGaO4 (where a crossing of the magnetization curves is also observed) and AErSe2 (A = Na and K)27,30.

Fig. 3: Magnetization data for Ba3Er(BO3)3.
figure 3

a Isothermal magnetization data of Ba3Er(BO3)3 taken at various temperatures. Closed points show data with field applied in-plane and open points show data with field applied along the c-axis. The different colors indicate the temperatures of the measurements. b First derivative of isothermal magnetization with field applied along the c-axis showing anomaly between 1 T and 3 T, with colors representing the same temperatures as in a.

Full size image

Heat capacity

To further characterize this material, we also performed heat capacity measurements at various fields from room temperature down to 60 mK using helium-4 and dilution refrigerator (DR) setups. The measurements were carried out on a single crystal sample of Ba3Er(BO3)3 with field applied parallel to the c-axis. The heat capacity data up to 15 K with several applied fields (0 T, 0.25 T, 0.5 T, 1 T, 3 T) are shown in Fig. 4a. Our heat capacity setup requires thin, flat samples where field is applied perpendicular to the surface of the sample. Ba3Er(BO3)3 cleaves easily between the ab-planes, so we were unable to produce a sample with the correct dimensions and orientation to measure the heat capacity with field applied in-plane.

Fig. 4: Field-dependent thermal measurements of Ba3Er(BO3)3.
figure 4

Colors represent different fields applied during measurements. a Heat capacity of Ba3Er(BO3)3 at various fields applied along c-axis, up to 15 K. b Low-temperature heat capacity data with field applied in smaller increments. c Magnetic heat capacity of Ba3Er(BO3)3 obtained by subtracting T3 fit to the 0 T data. d Magnetic entropy up to 20 K. The dashed line indicates entropy of \(R\ln 2\). Note that c and d are plotted on a log scale to better show the low temperature behavior.

Full size image

The zero field data shows a sharp peak at 0.1 K, possibly indicating the onset of magnetic order. The peak broadens and shifts to higher temperature with increasing field strength, and appears more like a two-level Schottky anomaly. In order to further study the transition between the sharp peak and broad features, we measured the heat capacity while increasing the applied field in small steps. For clarity, the fields below 1 T are shown separately in Fig. 4b for the low-temperature region up to 1.5 K. As the field is slowly increased up to 0.1 T, the intensity of the peak decreases but the peak’s position at 0.1 K does not change. As the field is further increased beyond 0.15 T, the peak begins to broaden and shift to higher temperatures, matching the behavior we observed for the higher fields. We are planning a follow-up neutron powder diffraction measurement to confirm the nature of the ordering and determine the magnetic structure at low temperatures. We consider this to be outside the scope of the current paper and plan to present those results in a future publication.

The magnetic heat capacity and corresponding magnetic entropy are shown in Fig. 4c, d. The lattice contribution to the heat capacity was estimated by fitting the zero field data to a cubic model (T3), which was subtracted to obtain the magnetic heat capacity (Cmag). The magnetic heat capacity of Ba3Er(BO3)3 reveals a low-temperature peak that shifts to higher temperature with increasing field, which lines up with the Schottky-like feature we observe. We also observe a peak centered at 6.5 K which is present for all fields measured.

The entropy as a function of temperature [see Fig. 4d] shows a plateau at around 2 K, with a value of \(R\ln 2 \sim 5.76\,\,{{\mbox{J}}}\,\,{{{\mbox{mol}}}}^{-1}{{{\mbox{K}}}}^{-1}\). This indicates that there is an effective two-level state at very low energy, agreeing with the prediction of a Kramers’ doublet ground state39,40. Furthermore, another plateau is observed at about 10 K, with an entropy value about \(\frac{4}{3}R\ln 2 \sim 7.68\,{{\mbox{J}}}{{{\mbox{mol}}}}^{-1}{{{\mbox{K}}}}^{-1}\). If we assume that this peak is due to another CEF doublet, we should expect to see the entropy increase by \(R\ln 2\) above the previous plateau. However, the entropy only increases by \(\frac{1}{3}R\ln 2\) over the ground state doublet, so this doublet must only be associated with one-third of the Er ions. We associate this CEF excitation with the Er atoms at Wyckoff position 2a, which indeed comprise one-third of the total Er atoms.

Neutron scattering measurements

To get a better understanding of the CEF spectrum of Ba3Er(BO3)3, inelastic neutron scattering measurements were carried out on a powder sample at the SEQUOIA time-of-flight spectrometer at Oak Ridge National Laboratory48. Measurements were performed with incident energies 150 meV, 50 meV, and 12 meV, at temperatures 5 K, 30 K and 150 K. The neutron spectra for all data sets are given in Supplementary Fig. 2. Er 3+ has \(S=\frac{3}{2}\)—so all of its energy levels are Kramers doublets—and has \(J=\frac{15}{2}\)—so one would expect to see a total of 7 excitations above the ground state. Our data show two low-energy excitations at 0.8 and 1.9 meV and two higher-energy excitations around 10 meV as shown in Fig. 5a. Additionally, we see signatures of two closely lying excitations at higher energies of approximately 12.5 and 13.5 meV—see Supplementary Fig. 2. The neutron spectrum with Ei = 12 meV taken at 5 K and a line cut to show the intensity as a function of energy are shown in Fig. 5. The first excitation at 0.8 meV corresponds to a temperature of 9.3 K which coincides with the broad feature we observe in the magnetic heat capacity measurements [see Fig. 4c]. The next excitation occurs at 1.9 meV which corresponds to a temperature of 22 K, which matches with the small peak in the magnetic heat capacity we see at the highest temperatures of our measurement. The CEF spectrum we observe is similar to the values calculated for ErMgGaO4 and Ba3ErB9O18, (although the excitations tend to appear at slightly lower energies in Ba3Er(BO3)3 than in either of these materials)38,49, and for those observed in NaErS2 and KErSe239,40.

Fig. 5: Neutron scattering data of Ba3Er(BO3)3.
figure 5

a Inelastic neutron scattering (INS) data taken on a powder sample of Ba3Er(BO3)3 at 5 K. b Blue points show integrated intensity from the boxed region in a .The red line shows the Stevens operator fit to the data.

Full size image

The low energy levels observed in the magnetic heat capacity and inelastic neutron scattering are captured by the CEF Hamiltonian H. The two Er sites have different point symmetries, so depending on the local symmetry being C3 (for 4b positions) or C3v (for 2a positions), we have

$${H}_{{{{{{{{\rm{4b}}}}}}}}}= \, {B}_{2}^{0}{O}_{2,0}+{B}_{4}^{0}{O}_{4,0}+{B}_{4}^{3}{O}_{4,3}+{B}_{4}^{-3}{O}_{4,-3}+{B}_{6}^{0}{O}_{6,0}+{B}_{6}^{3}{O}_{6,3}\\ +{B}_{6}^{-3}{O}_{6,-3}+{B}_{6}^{6}{O}_{6,6}+{B}_{6}^{-6}{O}_{6,-6}$$
(1a)
$${H}_{{{{{{{{\rm{2a}}}}}}}}}={B}_{2}^{0}{O}_{2,0}+{B}_{4}^{0}{O}_{4,0}+{B}_{4}^{3}{O}_{4,3}+{B}_{6}^{0}{O}_{6,0}+{B}_{6}^{3}{O}_{6,3}+{B}_{6}^{6}{O}_{6,6}$$
(1b)

Here On,m are the standard Stevens operators50, and \({B}_{n}^{m}\) are the Stevens parameters to be determined. We perform a point charge calculation to obtain the Stevens parameters and the CEF levels51 (for further details on the calculation, see Supplementary Note 3 and Supplementary Fig. 4). See ref. 52 for derivation. For the ion positions we use the data from the Materials Project53,54. We consider the following two models:

  • The Ba-Er-B-O model: We calculate the CEF contribution from all types of ions.

  • The O-only model: We only include the CEF contribution from the oxygen ions O2−. This describes the scenario in which the CEF is mainly due to the ligands.

The results are shown in Table 3. The Ba-Er-B-O model gives comparable CEF levels for Er3+ in both the 2a and 4b environments. On the other hand, the O-only model has the first CEF excitation for the 2a position at E1 = 0.2322 meV (2.7 K) and for the 4b position at 1.156 meV (13.4 K); as a result, it provides one scenario where the 2a has a significantly lower first CEF excitation, in agreement with the CEF result from the entropy analysis. We point out that the point charge model so far involves no fitting parameters, yet already shows qualitative agreement with the experimental result.

Table 3 Experimental and theoretical comparison.
Full size table

To further quantitatively improve our result, we include a radial factor γ in the O-only model to control the ion’s radius. As one can see, the O-only model with a radial factor γ = 1.1 reproduces the low energy CEF levels extracted experimentally quite well. For comparison, we also fit CEF Hamiltonian directly to our neutron data, using the point charge calculations as a guide for which peaks to fit to which ion. The result of that fit is also given in Table 3 for comparison and the fit is shown in Fig. 5b. However, we note that this fit is likely unreliable: since we only observe 4 distinct peaks in our data, the fit is largely unconstrained and gives several additional transitions at higher energies which we do not observe in our neutron data (details on the data fitting and comparison of the three models to the neutron data are given in Supplementary Note 2, and the results are shown in Supplementary Fig. 3). While the O-only model seems to give a good estimate of the CEF levels, we point out two potential issues: (1) it is yet to find a reason why keeping only the anions (O2−) improves the modeling. Naively this describes the situation where the positive charges from the cations are maximally screened; such screening effect is beyond a point charge description. (2) We note that the Stevens parameters from the point charge model show a large difference from the experimental values (particularly in the values of the 4a Er site). More experimental data are needed in order to better understand the crystal field physics. To summarize, a point charge model for the CEF levels from only O2− contributions qualitatively reproduces the magnetic heat capacity at low temperature due to the ground state Kramers pair and low energy first excitation Kramers pair only for Wyckoff position 2a.

Exchange physics

The low energy physics in the temperature range 0 K < T < 11 K happens in the ground state doublet and the first excited doublet for the 2a position ions. If we restrict to a smaller temperature range 0 K < T < 2 K, the physics happens purely in the ground state doublet. In this section we focus on the smaller temperature range 0 K < T < 2 K. Our goal is to analyze the possible spin Hamiltonian consistent with the experimentally observed magnetic heat capacity. At nearest-neighbor level there are two types of exchanges: those connecting two 4b sites, which we call J4b, and those connecting 2a and 4b, which we call J2a. See Table 2 for more information. The triangular lattice can be viewed as a honeycomb lattice of 4b ions decorated by a (larger) triangular lattice of 2a ions located at the hexagon centers; see Fig. 1c for an illustration. The exchange pathways for J2a and J4b are both Er-O-B-O-Er, although the structure of those pathways are different (see Fig. 1a). In addition to the different symmetries of the two Er sites, the O-B-O bond angles are different along the two exchange pathways: between the 2a and 4b sites both angles are 125∘, while the two angles between the two 4b sites are 105∘ and 109∘. Because the 2a and 4b sites have different crystallographic environments, the exchanges J4b and J2a can have different energy scales. This phenomenon of different exchange parameters for different Wyckoff sites was calculated for another hexagonal rare-earth compound GdInO3, where it was hypothesized that the different interactions between the sites lead to reduced frustration55. Based on the magnetization measurements, we estimate that the exchange interactions should be smaller than 0.1 meV, but further measurements would be needed to confirm this.

2D magnetic order

Two pieces of information can be extracted from the magnetic heat capacity data in Fig. 4b. First, the peak at ~0.1 K does not immediately vanish at small magnetic field, resembling an antiferromagnetic ordering transition. Second, there is a lower, broader peak at ~0.3 K, evident at small magnetic field, that moves towards higher temperatures as the field increases. This trend of moving to higher temperature with increasing field suggests that the field is assisting rather than competing with the ordering. Therefore, we associate this broader peak with ferromagnetic fluctuations of the system. Based on this analysis, we propose that Ba3Er(BO3)3 realizes two-sublattice exchange physics, in which the honeycomb lattice spins develop ferromagnetic correlations due to the additional spins at the hexagon centers but eventually order in an antiferromagnetic phase.

3D magnetic order

Even if the magnetic order is determined within each layer, nontrivial ordering can still happen between the layers. To study the possible 3D magnetic order we now determine the magnetic symmetry that is consistent with the 2D magnetic order. Within each layer, the magnetic order has three-fold rotation symmetry C3 and the magnetic mirror symmetry \({M}_{x}^{{\prime} }=M{{{{{{{\mathcal{T}}}}}}}}\) (product of the ordinary mirror reflection M given in Table 1 and time reversal \({{{{{{{\mathcal{T}}}}}}}}\)). Furthermore, the two layers of Er3+ in one unit cell can still have different order configurations: the order on the two layers can be either related by 21 (the ordinary screw in Table 1) or \({2}_{1}^{{\prime} }={2}_{1}{{{{{{{\mathcal{T}}}}}}}}\). Physically, these two cases correspond to antiferromagnetic and ferromagnetic (collinear) order between the layers. The corresponding magnetic space groups are P6\({}_{3}^{{\prime} }\) cm and \({{{{{{{\rm{P}}}}}}}}{6}_{3}{{{{{{{{\rm{c}}}}}}}}}^{{\prime} }{{{{{{{{\rm{m}}}}}}}}}^{{\prime} }\)56, respectively. Note that both cases allow an effective Zeeman field on the 2a spin:

$${H}_{{{{{{{{\rm{Zeeman}}}}}}}}}={h}^{z}\mathop{\sum}\limits_{i}{S}_{i,{{{{{{{\rm{2a}}}}}}}}}^{z},$$
(2)

where hz is the effective Zeeman field allowed by the magnetic space group symmetry at the 2a site. The isothermal magnetization measurement (Fig. 3) observes zero magnetization at zero field for temperatures all the way down to 0.3 K (the lowest temperature reached). One scenario consistent with this is that the two layers order antiferromagnetically, corresponding to the magnetic space group P6\({}_{3}^{{\prime} }\) cm.

Exchange interactions

In this section we discuss the influence of exchange interactions in the context of the experimental data. Based on the structure and symmetry, a minimal exchange model consists of the Ising hamiltonian

$${H}_{{{{{{{{\rm{Ising}}}}}}}}}={J}_{{{{{{{{\rm{2a}}}}}}}}}\mathop{\sum}\limits_{i\in {{{{{{{\rm{2a}}}}}}}},j\in {{{{{{{\rm{4b}}}}}}}}}{S}_{i}^{z}{S}_{j}^{z}+{J}_{{{{{{{{\rm{4b}}}}}}}}}\mathop{\sum}\limits_{i,j\in {{{{{{{\rm{4b}}}}}}}}}{S}_{i}^{z}{S}_{j}^{z}+{H}^{z}\mathop{\sum}\limits_{i}{S}_{i}^{z},$$
(3)

with antiferromagnetic interactions J2a > 0, J4b > 0. (We note that a related model with Heisenberg exchange interaction57 may also be relevant to our compound Ba3Er(BO3)3). This model contains several simpler limits. First, for J2a = J4b it reduces to the triangular lattice Ising model, which is a canonical problem in statistical physics, famous for its frustration and lack of magnetic order. Second, if we include only J4b > 0 exchange, it reduces to a honeycomb lattice Ising model on the 4b sites, which is unfrustrated and shows an antiferromagnetic ordering transition. The phase diagram of the model (3) at zero-field (Hz = 0) can be obtained in closed form and contains three phases58. At nonzero field Hz > 0, the phase diagrams of the two limiting cases above both contain an ordered phase and a disordered phase59,60.

We now discuss the experimental heat capacity shown in Fig. 4 in context of exchange physics. Ideally, one would like to explain the following key features of the data: At low field there is a singular peak at about 100 mK presumably reflecting long range ordering. This peak decreases in magnitude with weak fields, and becomes noticeably broader above 0.1 T. Above this field, the now broadened low-temperature peak moves upward in temperature with increasing field. In addition, a second distinct broad high temperature (at around 6.5 K in zero field) peak is present at all fields up to about 2 T, when the peak originating at low temperature merges with it.

The persistent sharpness of the low-T peak in non-zero fields below 0.1 T, shown in Fig. 4b, suggests that it is antiferromagnetic in nature, distinct from the peak in an Ising ferromagnet which broadens immediately upon application of a field. By contrast, in an antiferromagnet, the ordering peak typically weakens and initially moves to lower temperature upon application of a field, due to the competition of antiferromagnetism, which wishes to anti-align the moments and the field, which tries to align them. At somewhat larger fields, the broadened peak typically reverses direction and increases its temperature, as the Zeeman term increases the thermal energy required to release the magnetic entropy. These behaviors are common to many antiferromagnetic models: The low field suppression of ordering temperature happens in the honeycomb lattice Ising model61, and the high field enhancement of the heat capacity peak and peak temperature happens in both the honeycomb lattice and the triangular lattice Ising models62.

In the experiment, we indeed observe an initial weaking of the ordering peak, but not a shift to lower temperatures. This is perhaps due to the very low temperature of the measured ordering transition. A broadened low-temperature peak moving to higher temperatures for H > 0.15 T is consistent with the expectations just described. As discussed earlier, the high-temperature 6.5 K peak may originate from the crystal field excitations of the 2a sites.

Given the rough consistency of the picture of weak antiferromagnetism, we briefly consider the nature of the possible ordered state. We further introduce an Ising exchange coupling between the two layers, and the total Hamiltonian is

$${H}_{{{{{{{{\rm{total}}}}}}}}}={H}_{{{{{{{{\rm{Ising}}}}}}}},u}+{H}_{{{{{{{{\rm{Ising}}}}}}}},d}+{J}_{u-d}\mathop{\sum}\limits_{i}{S}_{i}^{z}{S}_{j}^{z}.$$
(4)

In general this model results in three-dimensional ordering. In the case of zero field, its phase diagram is a simple extension of the 2D phase diagram58. We conclude that for J4b > J2a, Ju−d > 0, it supports the following candidate three-dimensional magnetic structure below the zero field transition: Néel order on the honeycomb lattice decorated with extra spins aligned ferromagnetically on the hexagon centers, while the layers order antiferromagnetically with a two-layer periodicity. Further magnetic measurements below the 100 mK ordering transition will be required to test and refine this proposed structure.

Conclusions

We successfully synthesized high-quality single crystal samples of Ba3Er(BO3)3, a compound with an erbium-based triangular lattice. Our low-temperature magnetic measurements unveiled significant anisotropy within this system. The magnetic entropy at zero field, as a function of temperature, displayed two distinct plateaus-one at 2 K and another at 10 K, corresponding to \(R\ln 2\) and \(\frac{4}{3}R\ln 2\) respectively. We’ve linked these observations to the unique environments occupied by the magnetic ions Er3+ (specifically, their positions at Wyckoff positions 2a and 4b), which possess distinct symmetry properties and crystal field environments. Our calculations align with this interpretation. Additionally, we delved into the exchange interactions occurring within the ground state doublet. Moreover, our research highlights the presence of a zero-field ordering transition, which we contend is of an antiferromagnetic nature. Several characteristics of the heat capacity can be elucidated within this context.

Methods

Sample synthesis and single crystal growth

Polycrystalline Ba3Er(BO3)3 samples were synthesized using the solid-state reaction technique. BaCO3 (99.95% metal basis, Alfa Aesar), H3 11BO3 (>99 atom %, Sigma Aldrich), and Er2O3 (99.95%, Alfa Aesar) were used as starting precursor materials. The  11B isotope of boron was used in the synthesis to allow for neutron measurements with these samples, due to its lower neutron absorption cross section compared to  10B. Using stoichiometric ratios of the starting materials resulted in some leftover Er2O3 after sintering, so 10% excess of H3BO3 and BaCO3 by weight were added to account for their loss during synthesis. These were mixed thoroughly in a mortar and pestle and then pressed into a pellet. The pellets were sintered at 500 ∘C for 12 h to remove any water the powders could have absorbed, and then at 1000 ∘C for 48 h to complete the reaction. As needed, samples were resintered at 1000 ∘C with intermediate grindings to achieve the pure phase. Once the pure phase was obtained (see Supplementary Note 1, Supplementary Fig. 1), the powders were pressed into cylindrical rods using a hydrostatic press under 700 bar pressure, in preparation for single crystal growths. The cylindrical rods (known as feed and seed) were sintered at 1130 ∘C under O2 atmosphere in a vertical tube furnace to achieve higher density. These rods were then grown into single crystals using a four-mirror optical floating zone furnace. Growth parameters were refined over the course of several growth attempts and eventually high purity centimeter-sized crystals were grown.

Magnetic and heat capacity measurements

Magnetic susceptibility and magnetization were measured from 300 K down to 300 mK using a superconducting quantum interference device (SQUID) magnetometer with  3He attachment. Further measurements were carried out using a Quantum Design Dynacool Physical Property Measurement System (QD PPMS). Heat capacity data were taken from room temperature down to 60 mK using helium-4 and dilution refrigerator (DR) options attached to the PPMS.

Inelastic neutron scattering

Inelastic neutron scattering (INS) experiments were performed on a powder sample of Ba3Er(BO3)3 on the SEQUOIA spectrometer at Oak Ridge National Laboratory48 to determine the crystal electric field (CEF) ground state of the compound. Measurements were taken at 5 K, 30 K, and 150 K using incident energy 12 meV, 50 meV, and 150 meV at each temperature. The powder was sealed with helium exchange gas to ensure good temperature coupling. A powder sample of the nonmagnetic analog Ba3Lu(BO3)3 was also measured to use for nonmagnetic background subtraction. The neutron data were analyzed using DAVE63 and the CEF levels were fit using the Python package PyCrystalField64.