Geometrically frustrated magnetic systems offer a platform to explore magnetic states with suppressed long-range ordering and unconventional excitations [1,2]. Fractional magnetization plateaus are one of the most remarkable consequences of geometrical frustration [3], but in spite of over three decades of intense research, the discovery of plateaus in a new compound often poses a challenge to explain. The most famous example is the Shastry-Sutherland Heisenberg compound SrCu 2 (BO 3 ) 2 and its improbable sequence of plateaus at 1/8, 2/15, 1/6, 1/4, 1/3, 2/5 and 1/2 [4][5][6][7][8][9] that resisted for 15 years until a plausible explanation was put forward [10]. These magnetization plateaus have been associated with a distinct set of exotic super-lattice spin structures originating from the bosonic crystallization of the triplets or of boundstates of triplets. However, due to the high magnetic fields required for these transitions, such as 27.2 T for the first plateau, the direct evidence for the actual spin structure remains to be measured using probes such as neutron scattering. The difficulty in reaching the magnetic fields required for the higher magnetization plateaus in SrCu 2 (BO 3 ) 2 and the desire to study the implications of the geometric frustration in a different setting than that of the spin-1/2 Heisenberg model motivates the search for other magnetic compounds hosting the Shastry-Sutherland lattice (SSL).
The SSL [11] can be viewed as an orthogonal arrangement of dimers, and this structure has already been realized in several rare-earth compounds where the spins are Ising-like, but a clean realization of the Ising model with just interand intra-dimer couplings is still missing. The interest in that model lies in the very solid prediction of a unique 1/3 magnetization plateau [12,13], a prediction that has not been verified so far. In the ReB 4 family, long-range Ruderman-Kittel-Kasuya-Yosida (RKKY) interactions give rise to various fractional magnetization plateaus: a 1/5 plateau in NdB 4 ; a 1/3 plateau with Up-Up-Down (UUD) ferrimagnetic ordering, as well as narrow 1/2 and 3/5 plateaus in HoB 4 [14]; and multiple plateaus, including a 1/2 plateau, in TmB 4 [15,16]. Conversely, the insulating compounds BaNd 2 ZnO 5 [17] and BaNd 2 ZnS 5 [18] exhibit ferromagnetic intra-dimer interactions stabilizing a double-Q magnetic structure without fractional plateaus.
Recently, a novel family of insulating rare-earth-based SSL compounds, the melilites Re 2 Be 2 GeO 7 , has been reported, crystallizing in the tetragonal P 42 1 m space group [19]. In these compounds, Re magnetic ions form SSL orthogonal dimer planes (see FIG. 1(a)). Recently, this family has garnered significant interest, sparking active research efforts. For instance, work performed on Nd 2 Be 2 GeO 7 has revealed both short-range spin correlations and long-range magnetic ordering, while studies on Pr 2 Be 2 GeO 7 have identified dynamic spin-freezing behavior [20]. Similarly, investigations of Yb 2 Be 2 GeO 7 indicate the absence of long-range order, pointing to a possible quantum spin liquid state [21]. Despite these advances, comprehensive magnetic property measurements across this family remain scarce, leaving many fundamental aspects yet to be explored.
In this paper, we present Er 2 Be 2 GeO 7 as the closest realization of the Ising model on the SSL to date. As an insulating system, it lacks long-range RKKY interactions, and its dimers adopt an antiferromagnetic configuration in zero field, consistent with the expectations for an Ising model with strong intradimer interactions. However, when a magnetic field is applied along the [001] direction-parallel to the Ising spins-the magnetization curve reveals two well-defined plateaus at 1/4 and 1/2 of the saturation magnetization, with no sign of the conventionally expected 1/3 plateau. This finding is entirely unexpected, as all known variations of the Ising model on the SSL predict either a single 1/3 plateau for short-range interactions or multiple plateaus for long-range interactions, but never exclusively the 1/4 and 1/2 plateaus. Notably, these unexpected plateaus emerge at experimentally accessible magnetic fields, making their detailed investigation both feasible and compelling.
Intrigued by this unexpected result, we conducted a theoretical exploration of several new model variants and revisited the sample's structure using single-crystal neutron scattering. This combined investigation revealed a clear and elegant solution: a subtle orthorhombic distortion in the lattice. This distortion induces a slight asymmetry in intra-dimer interactions along the two orthogonal directions, destabilizing the expected 1/3 magnetization plateau and giving rise to the observed 1/4 and 1/2 plateaus.
The crystal structure and phase purity of Er 2 Be 2 GeO 7 were thoroughly analyzed using powder X-ray diffraction and Rietveld refinement (see Appendix B). No sign of an impurity phase was detected. Detailed results, including the single crystal synthesis, can be found in the Appendix B. Temperature dependent magnetic susceptibility and isothermal magnetization measurements were conducted on high-quality single crystals of Er 2 Be 2 GeO 7 for both H ∥ [001] and H ∥ [010] directions.
In FIG. 1 (b,c), we present isothermal magnetization measurements along both crystallographic directions, which confirm the presence of strong directional magnetic anisotropy. A pronounced Ising-like anisotropy is evident, with a domi-nant spin component aligned along the [001] direction. This significant anisotropy highlights the restricted spin dynamics characteristic of an Ising system. At 2.5 K, M [001] approaches 7.39 µ B at 7 T, while M [010] reaches 4.5 µ B , confirming that the direction [001] has the largest projection of the easy axis. A more quantitative analysis will be performed in the theory section, after the Van Vleck contribution to the magnetization has been calculated and subtracted. We take the saturation for the [001] direction as 7.39 µ B , which can be used to label the field-emergent magnetization plateaus as fractions of M S . Additionally, Tunnel Diode Oscillator (TDO) measurements, detailed in Appendix C, showed no evidence of any further magnetic transitions up to 35 T.
In FIG. 1 (d), the low-temperature magnetic susceptibility χ [001] initially increases to 2 K, then decreases to around 0.85 K, indicating an antiferromagnetic transition. In contrast, χ [010] increases sharply from 2 K to 0.85 K and remains constant below this transition temperature (see FIG. 1(e)), indicating ferromagnetic behavior in the [010] direction. In FIG. 1 (f), we present M [001] (H) below the ordering temperature at 0.3 K and the corresponding dM/dH curve. Interestingly, we observe the emergence of multiple fractional magnetization plateaus as a function of field. These plateaus correspond to distinct magnetic phases that are stable within different field ranges and appear at 1/4 and 1/2 of the saturation magnetization. The critical fields associated with these plateaus are identified by differentiating M [001] (H) with respect to H and using the peaks of the dM/dH curve to determine the start and end of the plateaus. The critical fields associated with the start and end of the 1/4 plateau are 0.23 T to 0.33 T, and for the 1/2 plateau, they are at 0.33 T and 0.56 T. It should be noted that the critical fields for these plateaus are about two orders of magnitude smaller than those required for SrCu 2 (BO 3 ) 2 , making them more accessible for investigation using standard experimental techniques. For a magnetic field applied along the [010] direction at 0.3 K, M [010] (H) does not exhibit any magnetization plateaus. However, as shown in the inset of FIG. 1(e), M [010] (H) undergoes a rapid increase at a rate of 21.8 µ B /T for fields between 0 T and 0.07 T, reaching approximately 1.5 µ B and following a sigmoid-like curve. Remarkably, within the narrow field range H ∈ [-0.07, 0.07] T, the moment shifts linearly from -1.5 µ B to +1.5 µ B without any observable hysteresis. Beyond 0.1 T, M [010] (H) gradually approaches 4 µ B , with no further transitions detected.
The specific heat of the Er 2 Be 2 GeO 7 single crystal was measured under magnetic fields applied parallel to the [001] direction. Measurements were conducted at multiple fields ranging from 0 to 1 T, with fine steps in both field and temperature to capture the behavior across the narrow plateau phases. To optimize measurement time in the sub-Kelvin range, data collection primarily focused on the 0.4 K to 0.9 K range for most fields, while a few representative fields-one from each plateau-were measured over the full 0.055 K to 4 K range to ensure no significant features were overlooked. Panels in FIG. 2(a-d) display the magnetic heat capacity C m , grouped by field values corresponding to different plateau phases. In the m = 0 phase at 0 T, a sharp peak is observed at 0.88 K, consistent with the antiferromagnetic transition identified in the magnetic susceptibility. As the field increases to 0.23 T, this peak shifts to lower temperatures (down to ≈ 0.75 K), decreases in intensity, and broadens significantly. In the m = 1/4 regime, the broad peak continues to shift to lower temperatures while becoming even wider. Additionally, a new sharp feature emerges around ≈ 0.55 K, leading to a two-peak structure: a sharp peak at ≈ 0.55 K and a broader peak at ≈ 0.67 K. As the system transitions from the m = 1/4 phase to the m = 1/2 phase, these two peaks merge into a single, sharper peak centered at 0.55 K. Within the m = 1/2 phase, the sharp peak gradually flattens into a plateau-like feature at 0.5 T. Finally, in the fully polarized state, these features are completely suppressed. To better illustrate the field and temperature dependence of the specific heat, a color plot of C m (T, H) is presented in FIG. 2(e). Additionally, the magnetization M (H) at 0.3 K is shown to correlate the C m /T features with the different plateau phases. The magnetic en-tropy S m was determined by integrating the C m /T curves up to the highest available temperature. As depicted in FIG. 2(f), the entropy loss remains uniform across all plateau regimes. Furthermore, the magnetic entropy S m exceeds R ln(2), indicating additional contributions from the population of a lowlying crystal electric field (CEF) excited level above 3 K. This result aligns with the CEF analysis discussed in the following section.
To investigate the CEF levels of Er 2 Be 2 GeO 7 , inelastic neutron scattering (INS) measurements were performed on highpurity powder samples. According to Hund's rules, the ground-state multiplet of the Er 3+ ion is 4 I 15/2 , with a 2J + 1 = 16-fold degeneracy that splits into 8 doublets under the influence of the crystal field. The INS experiments revealed energy bands corresponding to transitions between these split CEF levels. FIG. 3(a,b) shows representative INS spectra at 6 K using E i = 30 meV and 60 meV. Six out of seven expected transition bands were identified, with the lowest CEF excitation observed at 1.58 meV. Measurements with higher energy coverage (E i = 150 meV) did not reveal any additional CEF levels. The CEF scheme was determined by fitting the single-ion CEF Hamiltonian as detailed in the Methods section. FIG. 3(c) and (d) present constant-Q cuts of the experimental data in the energy ranges ∆E = 0-25 meV (E i = 30 meV) and ∆E = 0-50 meV (E i = 60 meV), respectively. The experimental data are overlaid with the CEF fit (red line), which closely aligns with the observed spectra. Minor discrepancies in intensity may arise from unaccounted phonon contributions due to imperfect background subtraction. The insets of FIG. 3(c,d) show the experimental M (H) and χ(T ) along the [001] and [010] directions, overlaid with the CEF-calculated curves. The CEF calculations closely match the experimental data, capturing the anisotropy. Minor deviations may result from diamagnetic effects at higher temperatures and exchange interactions at lower temperatures. The fitted B n m parameters, energy levels, eigenvectors, and gtensor components are provided in the Appendix E. The analysis reveals that the easy axis of the magnetic moment lies in the plane defined by the c-direction and the dimer bond.
To determine the magnetic structure of Er 2 Be 2 GeO 7 below the transition at T N = 0.9 K, two neutron powder diffraction (NPD) spectra were obtained: one at 0.3 K below the transition and another in the paramagnetic phase at 100 K. Initial Rietveld refinement of the Bragg peaks at 100 K successfully refined all peaks within the space group P 42 1 m, yielding a low R wp of 5.02 % as shown in FIG. 4 (a). At 0.3 K, magnetic Bragg peaks became evident on top of the nuclear Bragg peaks, as illustrated in FIG. 4 (b), pointing to a magnetic propagation vector k = (0, 0, 0). High-intensity peaks at 2θ = 18.8 • and 26.7 • , indexed as (1 0 0) and (1 1 0), hint at a significant magnetic moment outside the a-b plane in Er 2 Be 2 GeO 7 . A magnetic symmetry analysis was performed to investigate the maximal magnetic subgroups compatible with the space group P 42 1 m with propagation vector k = (0, 0, 0) using the Bilbao Crystallographic Server (MAXMAGN program [22]). There are six k-maximal magnetic subgroups identified. The best fitting model is given by P 2 1 2 1 2 ′ (No. 18.19), where all Er atoms and the corresponding ordered moments are described by a single Wyckoff site. In this subgroup, the magnetic moment components m(a), m(b), and m(c) are independent. Allowing these components to vary freely during the refinement process resulted in m(a) and m(b) being comparable. The CEF analysis, as discussed in the previous section, indicates that the easy-axis anisotropy lies in the plane defined by the c axis and dimer bond direction, and thus we introduced the constraint |m(a)| = |m(b)| into the model. The model aligns well with the data, showing a magnetic R-factor R M value of 2.17 as shown in FIG. 4(b). The magnetic structure is characterized by antiferromagnetic dimer pairs with a dominant moment in the crystallographic c direction as illustrated in FIG. 4 (c). The ordered moment tilts away from the c axis by 24 • towards the dimer bond direction. Each Er ion exhibits an ordered moment of µ ≈ 7µ B . The model predicts a net ferromagnetic component in the ab plane. This agrees with the magnetic susceptibility behavior where χ [001] shows an antiferromagnetic transition but χ [010] shows a ferromagnetic transition at 0.83 K. The deduced magnetic moment components from the fit for all four atoms within the unit cell are documented in FIG. 4(d). We compare the magnetic moment determined in this way to the one predicted in the ground-state doublet of the CEF Hamiltonian (see appendix A 3 b):
for the Er 3+ site labeled as 1 in FIG. 4, where g J = 1.2 is the Landé g-factor (the magnetic moment of the remaining sites of the unit cell are related to that one by symmetry). NPD and CEF both predict magnetic moments pointing in the dimer plane and with dominant component along the z direction, but show a significant difference in the size of the moment (7µ B versus 5µ B ). A possible solution for this discrepancy is the enhancement of the ordered magnetic moment by exchangeinduced mixing with excited CEF levels. This phenomena has previously been reported in non-Kramers rare earth compounds such as Tb 2 Ti 2 O 7 [23] and PrRu 2 Si 2 [24], where the rare-earth ions develop a magnetic moment in spite of having a singlet CEF ground state.
Single-crystal neutron scattering measurements were performed using the HYSPEC spectrometer [25] at Oak Ridge National Laboratory (ORNL) to investigate the emergence of the magnetic structure under a magnetic field. The sample was oriented in the [h, k, 0] scattering plane with the field along the crystallographic c direction. In the paramagnetic phase at 12 K, weak (1, 0, 0) and (0, 1, 0) reflections were observed. For the tetragonal space group P 42 1 m, reflections such as (h, 0, 0) and (0, k, 0) with odd h and k are forbidden by the 2 1 symmetry element. The (1, 0, 0) reflection is approximately 5% of the intensity of the allowed (1, 1, 0) reflection, explaining why it was undetectable in NPD due to background noise.
To rule out the possibility of multiple scattering or higherorder wavelength contamination, a separate neutron singlecrystal diffraction experiment was conducted (see Appendix D), confirming that the forbidden peaks are intrinsic to the crystal structure. Two possible subgroups of P 42 1 m that allow these reflections were identified: P 4 (No. 81) and Cmm2 (No. 35). Among these, only Cmm2 was found to be consistent with the observed magnetization plateaus of this compound (see below). The limited number of observed reflections made it impractical to refine the magnetic structure fully. The low intensity of these forbidden peaks indicates minimal distortion, and the lattice constants satisfy a = b within the experimental resolution, supporting the validity of the highsymmetry P 42 1 m approximation. To approximate the magnitude of the distortion, the NPD data at 2 K was refined using the Cmm2 space group under the condition that it reproduces the observed intensity ratio I(100)/I(110) = 0.05, thus constraining the distortion within the experimental resolution. In the refined structure (see FIG. 5), Er 3+ ions were found to displace along the dimer bonds, splitting the average intra-dimer bond length of 3.3468 Å into 3.142 Å and 3.554 Å, leading to anisotropy in the intra-dimer bonds. In FIG. 5, we present the diffraction channel obtained by integrating the energy around the elastic line. The paramagnetic background at 12 K has been subtracted from the 60 mK data to isolate the magnetic contribution. In zero field, the diffraction pattern in the [h, k, 0] scattering plane displays magnetic peaks indexed by the magnetic propagation vector k = (0, 0, 0), in agreement with NPD data. Upon applying a field of H = 0.275 T, corresponding to 1/4 plateau, additional magnetic Bragg peaks emerge, indexable by k = (0.5, 0, 0) and (0, 0.5, 0). This indicates a cell doubling in the a and b directions with possibility of domain formation. Increasing the field to 0.45 T, corresponding to 1/2 plateau, results in a loss of intensity of k = (0.5, 0, 0) and (0, 0.5, 0) peaks with the emergence of faint (1/3, 0, 0) and (2/3, 0, 0) satellite peaks. Finally, in the polarized state (m = 1) at 7 T, only peaks at k = (0, 0, 0) remain. Although a full magnetic refinement could not be performed due to limited access to the peaks in the [h, k, 0] plane, theoretical calculations discussed below propose plausible magnetic structures for the m = 1/4 and m = 1/2 phases, which agree with the observed peaks.
The low-energy excitations observed in the inelastic channel are presented in Fig. 6. In the paramagnetic phase at 11 K and 3 K, a CEF band centered around ≈ 1.6 meV is observed, consistent with previous CEF analyses. At 3 K, as the temperature approaches the ordering temperature, this 1.6 meV band begins to split, possibly due to the molecular field generated by the onset of long-range spin correlations. At 0.06 K, a new flat-band excitation emerges at 0.35 meV, showing no detectable dispersion within the energy resolution. This energy band appears to be well-gapped. Under applied magnetic fields of 0.275 T and 0.45 T, corresponding to the m = 1/4 and m = 1/2 phases, the flat energy band broadens, with its intensity redistributing and shifting to 0.5 meV and 0.25 meV, respectively. Notably, magnetic Bragg peaks from the elastic line extend into the inelastic channel, creating continuum-like features at the corresponding (H, K, L) positions. Although low-energy excitation remains confined below 0.5 meV, the CEF Kramers doublet at 1.6 meV evolves across these phases (see Appendix E for line cuts), likely influenced by the development of a local ordered moment. Finally, at 7 T, the excitation shifts upward to approximately 1-1.25 meV, exhibiting weak dispersion. This behavior likely arises from a combina- tion of ferromagnetic spin order and field-dependent modifications to the CEF levels.
The two ground states of the CEF Hamiltonian have their largest components along the fully polarized states |±15/2⟩ (see table II in SI), suggesting that a good starting point to describe the magnetic properties of Er 2 Be 2 GeO 7 is the Ising model:
This is further supported by the absence of significant intermediate phases between the fractional magnetization plateaus (FIG. 1) and by the dispersionless mode in inelastic neutron scattering (FIG. 6). Recent studies on the theory of Rare-Earth Shastry-Sutherland compounds [23] have also shown that intra-dimer quantum fluctuations allow the existence of a m = 0 plateau at small fields only if the zero-field ground state is a singlet dimer product state. Given that both an m = 0 plateau and long-range magnetic order are observed in Er 2 Be 2 GeO 7 , the zero-field ground state may not be a dimer product state and instead is more likely to be an Ising order stabilized by farther-range interactions. The possible effect of residual off-diagonal interactions will be discussed later in this section.
The most remarkable property of Er 2 Be 2 GeO 7 is the existence of 1/4 and 1/2 magnetization plateaus in a field H ∥ [001]. This is in apparent contradiction with the Ising description since the Ising model on the SSL with intra-dimer coupling J 1 and inter-dimer coupling J 2 , is known to have a single fractional magnetization plateau at m = 1/3 [13]. To resolve this discrepancy within the Ising approximation, the model may be modified in two ways: by considering farther-range interactions or by allowing the Ising couplings to not be fully symmetric under the space group P 42 1 m.
We first consider the effect of long-range dipolar interactions on top of intra and inter-dimer exchange interactions:
The ground state of H Ising,LRI has been determined by systematic enumeration of magnetic unit cells and constrained energy minimization within each unit cell [26,27]. For methodological details, see section A 5. A typical result for the magnetization curve of this model is shown in Fig. 13. This approach shows that dipolar interactions lead to the stabilization of an infinite sequence of fractional magnetization plateaus (a devil's staircase), and that the 1/3 plateau remains the dominant one, in stark contrast with the properties of Er 2 Be 2 GeO 7 . These two aspects persist independently of the specific choice of couplings in H Ising,LRI .
The other possibility is to allow for the existence of a spatial anisotropy in the Ising couplings J zz ij . This hypothesis is supported by the single crystal neutron diffraction (FIG. 5-b), where a weak but noticeable intensity is observed at two reflections -(100) and (010) -which are forbidden in the P 42 1 m space group. These forbidden reflections persist at temperatures significantly higher than the energy scale of magnetic interactions and hence are likely of nuclear origin. The observed reflections are compatible with two of the maximal subgroups of P 42 1 m: P 4 and Cmm2. The distortions compatible with each of these subgroups induce a spatial anisotropy in the SSL Ising couplings, as illustrated in FIG. 7.
In order to test these anisotropic Ising models, the groundstate energy and magnetization have been determined by comparing ground-state energy lower bounds (obtained as described in section A 4) with the energy of the candidate states illustrated in FIG. 9 a magnetization curve compatible with that of Er 2 Be 2 GeO 7 with magnetization plateaus at m = 1/4 and m = 1/2. An arbitrarily small anisotropy in the intra-dimer interaction is enough to stabilize narrow m = 1/4 and m = 1/2 plateaus on the boundaries of the m = 1/3 plateau, while a relative anisotropy of 33% leads to the complete suppression of the m = 1/3 plateau (see the ground-state phase diagram in Fig. 8).
The role of intra-dimer anisotropy in destabilizing the 1/3 plateau can be intuitively understood by the following simple FIG. 8. Ground state phase diagram of the ASSLIM as a function of magnetic field and relative anisotropy (δJ1 = J ′ 1 -J1, J1 = (J1 + J ′ 1 )/2). The isotropic intra-dimer interaction J1 = 0.079 meV and inter-dimer interaction J2 = 0.013 meV were chosen to be consistent with the parameter estimates written in Fig. 11, while the no third-neighbor interaction J3 was considered.
argument. The 1/3 plateau structure is highly favorable in any Shastry-Sutherland system dominated by antiferromagnetic intra-and inter-dimer interactions because, in this structure, all up-down dimers gain inter-dimer exchange energy through their interactions with a nearby polarized dimer. By contrast, in the 1/4 resp. 1/2 plateau, some resp. all up-down dimers can be flipped and do not gain any energy through the interdimer exchange. However, the 1/3 plateau requires putting half the polarized dimers in each direction and is penalized by the intra-dimer anisotropy (J 1 -J ′ 1 ), while in the 1/4 and 1/2 plateaus all polarized dimers are in the same direction and can take advantage of the anisotropy. And since the 1/4 and 1/2 plateaus are degenerate with the 1/3 plateau at the lower and upper critical fields, even a small anisotropy immediately stabilizes 1/4 and 1/2 plateaus at the boundaries of the 1/3 plateau.
The identification of the 1/4 and 1/2 plateaus relies on having found one state that saturates the lower bound, but if one only includes J 1 , J ′ 1 and J 2 , the phases with magnetization m = 0, m = 1/4 and m = 1/2 are macroscopically degenerate (their residual entropies are listed in table I): In all plateaus, there are dimers that can flipped independently (shown in green in FIGS. 9-a and e). Besides, in the 1/4 plateau some rows can be shifted independently (see FIG. 10). However, the macroscopic degeneracy of the plateaus is incompatible with the experimental data: (i) With the specific heat since no hint of a residual entropy was found (FIG. 2-e); (ii) With neutron scattering, which revealed long-range order in all the plateaus. This indicates that residual interactions, either in the form of longer-range interactions or quantum fluctuations, must be considered.
Let us first look at the effect of longer-range interactions. Third-neighbor interactions (J 3 shown in FIG. 7-a) are enough to lift all dimer flip degeneracies, stabilizing the ordered states shown in FIG. 9 for the m = 0 and m = 1/2 plateaus. Note that the stripe ordered phase stabilized by J 3 in the m = 0 plateau leads to magnetic peaks either in the [100] or [010] reflections while both were observed in Er 2 Be 2 GeO 7 . This points to the coexistence of magnetic domains with net moments in each of the four allowed directions, which also explains why no net magnetization was observed at zero field in Er 2 Be 2 GeO 7 .
By contrast, the row shift degeneracy of the m = 1/4 plateau is remarkably resilient to longer-range interactions and subsists even if couplings up to the sixth nearest neighbor are considered. The row shift degeneracies eventually get lifted if full long-range dipolar interactions are taken into account in a distorted version of the SSL with Cmm2 symmetry, but in that case the ordered zig-zag state shown in FIG. 10-b is selected. This order is incompatible with neutron scattering since it would induce reflections at (h ± 1/4, k ± 1/4, 0) and (h±1/4, k∓1/4, 0), whereas magnetic peaks at (h±1/2, k, 0) and (h, k ± 1/2, 0) are found in the m = 1/4 plateau of Er 2 Be 2 GeO 7 . So another explanation has to be found for the 1/4 plateau. Interestingly enough, in the spin-1/2 Heisenberg Shastry-Sutherland compound SCBO the 1/4 plateau has the structure of Fig. 10-c [8], which would be compatible with the peaks observed in single-crystal neutron scattering in Er 2 Be 2 GeO 7 . So it is plausible that residual quantum fluctuations are responsible for the stabilization of this plateau in Er 2 Be 2 GeO 7 , but a proof that realistic off-diagonal interactions are sufficient remains to be seen. Note that for the 1/2 plateau the absence of a phase transition in the specific heat also suggests that quantum fluctuations, rather than longrange interactions, lift the degeneracy. Indeed, in both cases quantum fluctuations are expected to stabilize an entangled state (singlet or triplet) on the unpolarized dimers that would not be seen in neutron scattering [23].
Next we try to be more quantitative. Assuming, as usual in the context of rare earths, that only the nearest couplings J 1 , J ′ 1 , and J 2 are influenced by exchange, all other couplings are taken to be dipolar, which leads in particular to J 3 = 0.006 meV for the magnetic moments deduced from the CEF. The three critical fields that separate the plateaus at 0, 1/4, 1/2 and saturation can then be used to extract the nearest coupling, leading to J 1 = 0.092 meV, J ′ 1 = 0.066 meV and J 2 = 0.013 meV (see inset of Fig. 11). They have been obtained using the conversion h = (g zz /2)µ B H z = H z × 0.287 meV/T, where we make use of the effective g-tensor calculated in section A 3 b.
To check the consistency of our theory, we now discuss three other experimental sets of data: the flat mode in inelastic neutron scattering, the temperature dependence of the magnetization curve, and the temperature of the phase transition in zero field. In all cases, it turns out to be important to keep further dipolar couplings. The largest ones are actually the ferromagnetic nearest and next-nearest neighbor inter-layer couplings J 1 ⊥ = -0.021 meV and J 2 ⊥ = -0.006 meV. The flat mode observed in inelastic neutron scattering has an Model m = 0 m = 1/4 m = 1/2 m = 1
Dipolar interactions
Ground-state degeneracy (Ω), residual entropy (s0) and magnetic order (defined as the quotient between the high-and lowtemperature symmetry groups, G(T → ∞)/G(T = 0)) of the ASSLIM with and without J3 couplings. The row concerning the effect of quantum fluctuations results from qualitative considerations about the formation of dimer product states and the stabilization of the stripe order shown in FIG. 10 c.
energy given by 2J 1 +4J 3 -4J 1 ⊥ +4J 2 ⊥ = 0.24 meV, which is smaller than the experimental value of 0.35 meV but already much closer than if the inter-layer interactions were neglected.
Since J 1 , J ′ 1 ≫ J 3 , the melting temperature of the long-range stripe order at zero field in the absence of inter-layer coupling can be estimated by an effective model where the two antiferromagnetic states of each dimer are described by an Ising variable which interacts with its four neighbors due to J 3 , leading to the critical temperature:
significantly lower than the observed temperature of 0.88 K.
If one includes the dominant inter-layer couplings, the resulting effective model is an anisotropic 3D Ising model with anisotropy parameter ∆ = J 3 /(2J 2 ⊥ -2J 1 ⊥ ) = 0.22, for which Monte-Carlo simulations [28] determined the critical temperature:
The transition temperature corrected with J ⊥ is much closer to the experimentally measured temperature of 0.88 K. The persisting mismatch is probably due to other residual interactions.
Next, we calculated the temperature dependence of the magnetization for the purely 2D model using the Corner Transfer Matrix Renormalization Group (CTMRG), a numerical method briefly described in section A 6. To compare with the experimental data, we have subtracted the Van Vleck contribution to the magnetization (M
VV in eq. A19). The overall shape of the curve at low but non-zero temperature is compatible with the experimental data, further supporting the hypothesis that the system is well described by an Ising model, and that the broadening of the jumps between plateaus is thermal. Note also that, after the background is subtracted, the saturation magnetization of Er 2 Be 2 GeO 7 is compatible with the size of the magnetic moment of the CEF ground state doublets |±⟩ given in Eq.1. There is however a mismatch between the theoretical temperature, 0.15 K, needed to reproduce the experimental curve obtained at 0.3 K. This mismatch is again attributed to further couplings.
Finally, we comment on the most prominent feature of the specific heat data, a very high and narrow peak at the boundary between the 1/4 and 1/2 plateau. The interpretation of this peak, which has all the characteristics of a phase transition, is based on the following key observations : (i) There is no phase transition in the 1/2 plateau; (ii) There is a jump in magnetization at low temperature between 1/4 and 1/2, hence a first-order transition. (iii) The thermal phase transition in the 1/4 plateau takes place at 0.54 K, and it gives rise to a small peak that does not seem to connect to the very pronounced peak at 0.56 K and 0.34 T. Now, a first order transition can either terminate at a critical point, or be connected to a continuous phase transition through a tricritical point. In the fieldtemperature phase diagram (Fig. 2 (e)), the specific heat data are more consistent with a horizontal line at 0.54 K for the thermal transition of the 1/4 plateau phase cut by a first-order transition line terminating at a critical point that gives rise to a much stronger peak. So the most natural interpretation is that this prominent peak in the specific heat is a critical point terminating a first-order transition line starting at the transition between the 1/4 and 1/2 plateaus at zero temperature. A similar interpretation has been proposed for the peak observed in SrCu 2 (BO 3 ) 2 at the transition between the dimer phase and the plaquette phase [29].
Through a comprehensive combination of cutting-edge experimental techniques and advanced theoretical modeling, we have discovered and elucidated the unprecedented magnetic properties of Er 2 Be 2 GeO 7 , establishing it as a unique and powerful platform for exploring exotic phases of matter in Shastry-Sutherland lattice (SSL) systems. Our results mark the first observation of 1/4 and 1/2 fractional magnetization plateaus in an Ising SSL compound -an outcome that defies long-standing theoretical expectations of a 1/3 plateau and signals a fundamental shift in our understanding of such systems.
We demonstrated that this remarkable plateau sequence arises from a subtle orthorhombic distortion that introduces spatial anisotropy in dimer interactions, captured effectively by the anisotropic Shastry-Sutherland Ising model (ASSLIM). This model successfully reproduces the observed magnetization behavior and is validated by detailed structural and neutron scattering data, including the detection of symmetryforbidden reflections and anisotropic bond lengths. Notably, our experimental results reveal a complete absence of residual entropy across all phases, despite theoretical expectations of macroscopic degeneracy. This striking result is fully accounted for by incorporating two key factors: longerrange Ising interactions extending beyond second neighbors, and quantum fluctuations arising from off-diagonal terms in the effective Hamiltonian. Impressively, we find that degeneracy is lifted at zero field through third-neighbor interactions, while in the 1/4 and 1/2 magnetization plateaus, it is suppressed by the emergence of entangled dimer states on nonpolarized dimers-a direct manifestation of quantum fluctuations stabilizing unique ordered phases.
Our work provides a clear roadmap for the design and investigation of related rare-earth SSL systems, positioning Er 2 Be 2 GeO 7 and the broader melilite family as exceptional platforms to uncover and understand new phenomena in frustrated magnetism. The accessible energy scales and field ranges in Er 2 Be 2 GeO 7 , compared to other SSL compounds, open avenues for high-precision measurements and controlled theoretical investigations. FIG. 11. Ground-state energy (a) and magnetization (b) of Er 2 Be 2 GeO 7 after subtracting the calculated background contribution from the population of excited CEF levels, compared with the theoretical prediction of the ASSLIM model. The theoretical prediction was numerically computed using the Corner Transfer Matrix Renormalization Group method. Inset of (b) shows the values of interaction parameters required for agreement between of the critical fields predicted by the ASSLIM and those measured in Er 2 Be 2 GeO 7 , in comparison with the magnitude of dipole-dipole interactions. The third-neighbor interaction J3 was assumed to be entirely due to dipole-dipole interaction.
Looking ahead, our findings open several compelling avenues for future research that can deepen the understanding of Er 2 Be 2 GeO 7 and related systems. In particular, the distinct mechanisms responsible for lifting degeneracy across the various plateau phases-ranging from longer-range Ising interactions to quantum fluctuations-present a rich landscape for further theoretical and experimental investigation. Additionally, our work suggests that incorporating long-range interactions is essential to accurately capture the full energy and temperature scales, pointing to the value of refined modeling approaches. The observation of satellite peaks indicative of incommensurate correlations at the boundary between the 1/4 and 1/2 plateaus raises intriguing questions about emergent spin textures and competing orders in this regime. Moreover, determining the universality class of the phase transitions at zero field and within the 1/4 plateau remains an open challenge that could shed light on the critical behavior in frustrated Ising systems. Addressing these questions will benefit from a deeper understanding of the underlying microscopic interactions, potentially guided by ab-initio calculations and advanced numerical methods, paving the way for new discoveries in rare-earth-based Shastry-Sutherland systems.
Beyond this compound, our study underscores the enduring potential of the SSL geometry to promote unexpected phases and transitions, and it highlights rare-earth-based SSL materials as fertile ground for discovering and controlling exotic quantum states.
The research performed at Duke University was supported by NSF award number DMR-2327555. A portion of this research used resources at the High Flux Isotope Reactor and Spallation Neutron Source, a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory. A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by National Science Foundation Cooperative Agreement No. DMR-2128556 and the State of Florida. The research performed at EPFL has been supported by the Swiss National Science Foundation Grant No. 212082. The research performed at FAU was supported by the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus. For the numerical calculations of the unit-cell based ground-state search, we acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg. Y., M.E., and D.G. performed magneto-transport measurements; A.R., J.A.K., K.P.S., and F.M. provided theoretical interpretations; L.Y., A.R., J.A.K., K.P.S., F.M., and S.H. wrote the manuscript with comments from all authors.
Appendix A: Methods
Polycrystalline Er 2 Be 2 GeO 7 powder was prepared by a solidstate reaction route using the starting precursors of Er 2 O 3 (99.9 %, Alfa Aesar) with BeO (99.99 %, Alfa Aesar) and GeO 2 (99.99 %, Alfa Aesar). The starting precursors were weighed in the 1 : 2 : 1 molar ratio and the mixture was sintered at 1285 °C for 48 hours with intermediate grindings. Powder X-ray diffraction (PXRD) data were analyzed by performing Rietveld refinement using FullProf Suite. Once the phase purity of the powder sample was confirmed, the single crystals of Er 2 Be 2 GeO 7 were grown using four mirrors optical floating zone furnace (Model: FZ-T-12000-X-VII-VPO-PC, Crystal System Corporation, Japan). The grown crystals were analyzed and oriented using the Laue diffractometer (MULTI-WIRE LABS MWL120) and subsequently cut to the required dimensions using a wire saw. Heat capacity was measured on oriented single crystal samples of Er 2 Be 2 GeO 7 and powder samples of Lu 2 Be 2 GeO 7 (non-magnetic) using Helium-4 (1.8 K ≤ T ≤ 300 K) and dilution refrigerator (0.06 K ≤ T ≤ 2 K) set up attached to the Physical Properties Measurement Systems, Quantum design (PPMS Dynacool, QD, USA) accompanied by 14 T magnets. Magnetic measurements from 300 K to 2 K were performed using the Cryogenic Ltd SQUID (superconducting quantum interference device) magnetometer. Additional measurements from 1.8 K to 0.3 K were performed using the Helium-3 probe in SQUID magnetometer.
Neutron powder diffraction (NPD) measurements to study the magnetic structure of Er 2 Be 2 GeO 7 were performed using the HB-2A [30] diffractometer at the High Flux Isotope Reactor (HFIR) in Oak Ridge National Laboratory (ORNL). For the HB-2A experiment, 2.7 grams of the powder sample were loaded in an Al can and sealed under a Helium environment, and subsequently loaded onto a 3 He cryo-stick in a Cryostat. The measurements were performed using a collimation of open-21 ′′ -12 ′′ and a Ge monochromator to select incident neutrons of wavelength λ = 2.41 Å. Diffraction patterns were collected at 0.3 K and 100 K by counting for 120 sec at each point ranging between 2θ = 5 • and 130 • . Single crystal neutron scattering measurements were performed at the HYSPEC spectrometer [31] at ORNL. Approximately 1 gram of the single crystal sample of Er 2 Be 2 GeO 7 was oriented in the [h,k,0] scattering plane, with a magnetic field of up to 7 T applied along the crystallographic c-axis. The sample environment included a dilution refrigerator capable of reaching a base temperature of 60 mK. An incident neutron energy of 3.8 meV was selected using a Fermi chopper operating at 180 Hz. The paramagnetic phase dataset, collected at 12 K, was employed to subtract the nuclear component from the scattering. Inelastic neutron scattering experiments were performed on powder samples of Er 2 Be 2 GeO 7 and Lu 2 Be 2 GeO 7 (for phonon contributions) at SEQUOIA [32] spectrometer at ORNL to probe crystal electric field (CEF) levels. Powder samples were loaded in a Al can and sealed under 4 He atmosphere to ensure good temperature coupling. Measurements were collected at incident energy E i = 150 meV in high flux mode and at E i = 60 meV and 30 meV in high-resolution mode at 6 K, 50 K and 250 K. Data analysis was performed using the DAVE MSlice [33], and the PyCrystalField [34] packages, to determine the crystal electric field (CEF) parameters.
The crystal field effect on the magnetic ion can be parameterized using the single-ion crystal field Hamiltonian:
Here, B m n are the CEF parameters that parameterize the effects of the ligand environment on the magnetic ion, and O m n are the Stevens operators [34][35][36].
The fits to these observed CEF levels in the INS spectrum were carried out using PyCrystalField software [34] where B m n were optimized during the fitting procedure. The number of CEF parameters in the Hamiltonian was reduced to fifteen by rotating the coordinates so that the quantization axis z coincided with [001], while x and y were assigned to [110] and [1][2][3][4][5][6][7][8][9][10], respectively, corresponding to the Er1 site. This coordinate system ensures that the mirror axis [-0.707 0.707 0] lies in the xy plane and is parallel to y, ensuring that all imaginary CEF parameters are zero. The initial B m n were estimated using an electrostatic point charge model based on the 8 O 2- ligands surrounding the Er 3+ ion. Starting with these CEF parameters, 2D Q -∆E slices at all three temperatures (6 K, 50 K and 250 K) for both E i = 30 meV and E i = 60 meV were fitted. To further constrain the fit, magnetic susceptibilities (χ) along [001] and [010] and isothermal magnetization data at 5 K, 10 K, and 20 K in the same directions were included, ensuring the solutions captured the pronounced anisotropy observed in this system. The initial and the optimized B m n parameters are provided in the table in Appendix E. A comprehensive list of eigenvalues and their corresponding eigenvectors can also be found in the Appendix E.
a. Effective spin-1/2 Hamiltonian from the Crystal Electric Field At temperatures much smaller than the smallest CEF gap, T ≪ 1.63 meV/k B = 18.7 K only the lowest lying CEF levels, which we label as |+⟩ and |-⟩, will be populated and the system can be described by an effective spin-1/2 Hamiltonian. In order to derive this effective Hamiltonian, one must start from a Hamiltonian including the most general two-body interaction [37]
where
q (J i ) are spherical tensors of rank k, built as polynomials of order k in the angular momentum J i . They include the usual bilinear exchanges such as J i • J j as well as interactions between higher-order multipoles, up to the maximum rank k max = 2j. The components of the J α i are represented in the local reference frame xi , ŷi , ẑi , related to each-other by the symmetry operations of the space group of the lattice.
The effective spin-1/2 Hamiltonian can be obtained to leading order in J /∆ CEF by projecting H int in the ground-state manifold of H CEF
where P is the projector onto the ground-state manifold of H CEF and S α i are the effective spin-1/2 spin operators, defined as
where τ α σσ ′ are the Pauli matrices and the operators S α i act non-trivially only on site i. The spin-1/2 couplings J α,β i,j can be calculated from the physical Hamiltonian H int as
In some rare-earth compounds, such as the pyrochlore Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 [38]), H eff is dominated by diagonal terms proportional to J zz ij , which means that their lowtemperature magnetic properties can be well approximated by a classical Ising model. The dominance of the diagonal terms is a consequence of the fact that the strength of multipolemultipole interactions J kq,k ′ q ′ i,j -which result from a combination of physical processes such as magnetic dipole-dipole interactions, direct exchange, superexchange and lattice mediated interactions -are highly suppressed at rank k = 8 or higher. This fact, combined with the Wigner-Eckart theorem [39] which states that the only non-zero matrix elements in
implies that interactions can only connect states |m J ⟩ and |m ′ J ⟩ in first-order perturbation theory if |m J -m ′ J | ≤ 7. In the cases of Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 , the CEF ground state doublets are almost collinear with the fully polarized states, |±15/2⟩ and |±8⟩ respectively [40]. This means that, to first order in J /∆ CEF , the effective spin-1/2 Hamiltonian is an Ising model with interaction constants given by J zz i,j . Similarly, in Er 2 Be 2 GeO 7 the ground state doublets of the CEF have their largest component collinear with the fully polarized states |±15/2⟩ (see SM Table II), albeit to a lesser extent than in Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 . This justifies us in using an Ising model as a starting point to understand the magnetic properties of Er 2 Be 2 GeO 7 , and in invoking quantum fluctuations due to transverse exchange in the discussion of some properties.
Due to the Kramers degeneracy, there is some freedom in the choice of the basis of the ground-state manifold of H CEF . In the absence of a detailed model for the physical Er 3+ moments, we choose the basis of the CEF ground state doublet such that the projection of intra-dimer Heisenberg interaction leads to a diagonal coupling matrix:
The CEF states presented in table III reflect this choice of basis. If the generating tile and the cluster translation vectors (u 1 , u 2 ) are chosen such that every pair of interacting spins is contained in at least one cluster, the Hamiltonian can be decomposed as a sum over clusters. If different clusters overlap, some spins or bonds will be contained in several clusters and the decomposition is not unique. The possible decompositions are parameterized by the weights 0 ≤ α i,j ≤ 1 and 0 ≤ β i ≤ 1, responsible for splitting the energy of every site and interacting bond among the clusters that share it:
H loc (σ| tx,y , α, β) (A24)
For the decomposition to be consistent, the weights must satisfy the constraints:
x,y i,i ′ ∈t0,0
x,y i,i ′ ∈t0,0 j,j ′ ∈t0,0 δ i+xa1+ya2,i ′ δ j+xa1+ya2,j ′ α i,j = 1.
In the case of the Shastry-Sutherland Ising model, we take as generating cluster the set of 12 spins highlighted in fig. 12. Independent minimisation of the energy of each cluster leads to a lower bound on the ground state energy
Each value of (α, β) compatible with the constraints A26 and A27 defines a lower bound on the ground state energy. Therefore, the most restrictive lower bound can be obtained by optimising over the weights as
H loc (σ| t0,0 , α, β) .
If the ground state energy lower bound obtained in this way matches the energy of any of the candidate states, the bound is said to be saturated and we can be sure to have found the correct ground-state energy.
A recently proposed way to determine ground states for Hamiltonians of the form of Eq. ( 3) is to consider all possible magnetic unit cells up to a certain limitation and optimize an effective Hamiltonian on each of the unit cells with appropriately resummed couplings [26,27]. These couplings are chosen so that the energy per site of the pattern on the unit cell corresponds to the thermodynamic limit being tiled by this unit cell. The treatment of finite-range interactions in this framework is well-established [43]. The key insight that makes this approach possible for long-range interactions is that one can rewrite a diagonal long-range interaction for a periodic pattern with a K-site unit cell and translational vectors a 1 and a 2 i̸ =j
as sums over the unit cell of the pattern. The appropriately resummed couplings
with δ i,j being the Kronecker delta, L(a 1 , a 2 ) the lattice spanned by a 1 and a 2 and ζ L,α (x, y) being the Epstein ζfunction. Inverting the logic, this means that finding the ground state of (A32) on a unit cell gives the lowest energy state for the thermodynamic limit that fits the considered cell.
Therefore, the workflow is as follows: First, we determine the unit cells for the respective (non-)distorted SSL lattice with atomic positions provided by experimental measurements following the procedure described in [26]. Then we evaluate the resummed couplings for each unit cell using efficient implementations of the Epstein ζ-function [44]. In the end, a parameter study is conducted for varying the h values. Here, for each h value, an optimization of the spin state is performed on each unit cell and the configuration with the overall lowest energy in the thermodynamic limit is selected. From the resulting state, observables such as the magnetization or static structure factors can be computed.
This method is limited by the number and shape of the unit cells considered, as well as the optimization algorithm on the respective unit cells. For this study, we choose unit cells with no more than 16 elementary four site unit cells of the respective SSL lattice, with no more than six elementary unit cells in one linear direction.
... ... ... ... FIG. 14. Illustration of the Corner Transfer Matrix Renormalisation Group (CTMRG) algorithm.
The Corner Transfer Matrix Renormalization Group (CTMRG) algorithm [41,45,46] is a numerical method for 2D classical lattice models which relies on the approximate contraction of a tensor network formulation of the partition function. Since the cost of exact contraction of a 2D tensor network grows exponentially with system size, the CTMRG algorithm probes the thermodynamic limit by approximating the environment around a tensor by finite-dimensional corner (C i ) and edge (E i ) tensors accounting respectively for a quadrant or half-column of the infinite two-dimensional tensor network, as illustrated in Fig. 14. After the environment tensors are initialized with the desired boundary conditions, sites are iteratively added to the lattice until convergence has been reached. Once the environment tensors C i and E i have converged, the expectation value of a local operator can be calculated by performing a contraction similar to that of FIG. 14 where the tensor a in the center is replaced by a tensor representation of the operator.
The cut-off bond dimension χ to which the corner and edge tensors are truncated after each CTMRG step controls the precision of the algorithm. FIG. 16. Tunnel diode oscillator measurement showing the frequency response of the sample as a field up to 35 T is applied using resistive magnets at 300 mK. In the higher field range, specifically from 7 T to 35 T, no additional peaks are observed, indicating that no further magnetic transitions occur at higher fields.
ing a floating zone furnace. The phase purity of the grown crystals was confirmed using a powder x-ray diffraction pattern on crushed Er 2 Be 2 GeO 7 single crystal samples.
TABLE II. Crystal field parameters B m n calculated using the point charge model and after fitting the experimental data. The values are reported in units of meV. B m n Point Charge Model Fit B 0 2 × 10 3 -229 -2.611 B 1 2 × 10 1 -6.56 -4.736 B 2 2 × 10 1 2.19 1.325 B 0 4 × 10 4 -5.47 -0.670 B 1 4 × 10 3 -1.57 2.273 B 2 4 × 10 3 0.489 -2.988 B 3 4 × 10 4 41.05 -13.240 B 4 4 × 10 5 -96.54 -57.425 B 0 6 × 10 6 2.88 -0.889 B 1 6 × 10 4 0.214 1.101 B 2 6 × 10 5 0.070 -6.066 B 3 6 × 10 5 -1.21 7.818 B 4 6 × 10 5 -0.289 0.557 B 5 6 × 10 5 -1.42 7.547 B 6 6 × 10 6 was expected to be of the same order as the background. The full list of energy levels and the corresponding eigenvectors derived from the fitting is listed in Table III. The initial and the optimized B m n parameters are provided in Table II. Note that the the crystal-field parameters presented in table III were calculated in a coordinate system which has x parallel to [110] and z parallel to [001], as discussed in the main text. The basis of the CEF ground state was chosen such that the projection of intra-dimer Heisenberg interactions leads to a diagonal effective coupling matrix. The g-tensor calculated from the ground state eigenvectors is as follows : tion. The results, shown in Fig. 18, reveal a splitting pattern that is only partially captured by the single-ion CEF model. While the calculated spectrum correctly predicts the linear Zeeman splitting of the Kramers doublet, the experimental data show additional field-induced features that are likely due to internal molecular fields originating from the long-range magnetic order.
To verify that the forbidden (1,0,0) peak observed in the HYSPEC experiment is not due to multiple scattering or higher-order wavelength contamination, we conducted a neutron diffraction experiment at the WAND 2 diffractometer at
TABLE III. Eigenvalues and their corresponding eigenvectors for the single-ion CEF Hamiltonian, derived from fitting the Er 2 Be 2 GeO 7 INS data as discussed in the main text. Wavefunctions are presented in the |mJ ⟩ basis. E (meV) 15 2 13 2 11 2 9 2 7 2 5 2 3 2 1 2 -1 2 -3 2 -5 2 -7 2 -9 2 -11 2 -13 2 -15 2 0.000 -0.094 -0.109 0.08 -0.206 0.128 0.165 -0.039 -0.257 -0.117 0.251 -0.198 -0.1 0.283 0.513 0.061 0.589 0.000 0.589 -0.061 0.513 -0.283 -0.1 0.198 0.251 0.117 -0.257 0.039 0.165 -0.128 -0.206 -0.08 -0.109 0.094 1.674 -0.650 -0.334 -0.022 -0.322 0.047 0.189 0.149 -0.023 -0.254 0.081 0.175 -0.173 -0.305 -0.228 -0.155 -0.034 1.674 0.034 -0.155 0.228 -0.305 0.173 0.175 -0.081 -0.254 0.023 0.149 -0.189 0.047 0.322 -0.022 0.334 -0.650 6.240 0.208 0.246 -0.145 0.093 0.291 0.075 -0.320 -0.481 -0.009 0.342 -0.117 -0.272 -0.298 -0.293 -0.252 0.009 6.240 0.009 0.252 -0.293 0.298 -0.272 0.117 0.342 0.009 -0.481 0.320 0.075 -0.291 0.093 0.145 0.246 -0.208 15.604 0.001 0.169 -0.156 -0.005 -0.189 0.538 0.440 -0.369 0.199 -0.212 -0.132 0.348 0.080 -0.130 -0.209 0.050 15.604 -0.050 -0.209 0.130 0.080 -0.348 -0.132 0.212 0.199 0.369 0.440 -0.538 -0.189 0.005 -0.156 -0.169 0.001 23.191 -0.034 0.450 -0.137 -0.329 0.565 -0.070 0.396 0.374 0.155 0.130 -0.083 -0.025 0.025 -0.002 -0.008 0.000 23.191 0.000 0.008 -0.002 -0.025 -0.025 0.083 0.130 -0.155 0.374 -0.396 -0.070 -0.565 -0.329 0.137 0.450 0.034 28.002 0.001 -0.071 0.108 -0.011 0.076 -0.691 0.446 -0.497 -0.167 -0.131 -0.043 0.043 0.029 -0.030 -0.066 -0.006 28.002 0.006 -0.066 0.030 0.029 -0.043 -0.043 0.131 -0.167 0.497 0.446 0.691 0.076 0.011 0.108 0.071 0.001 45.056 -0.265 0.631 0.265 -0.384 -0.475 -0.157 -0.229 -0.097 0.012 0.017 0.028 -0.005 -0.005 0.001 0.003 0.000 45.056 0.000 -0.003 0.001 0.005 -0.005 -0.028 0.017 -0.012 -0.097 0.229 -0.157 0.475 -0.384 -0.265 0.631 0.265 62.195 0.000 -0.000 0.000 0.000 -0.004 -0.004 0.009 -0.011 -0.003 0.083 -0.157 0.269 -0.565 0.655 -0.206 -0.323 62.195 0.323 -0.206 -0.655 -0.565 -0.269 -0.157 -0.083 -0.003 0.011 0.009 0.004 -0.004 -0.000 0.000 0.000 0.000 ORNL at 3 K in the paramagnetic phase of the Er 2 Be 2 GeO 7 .
This experiment utilized a Ge(113) monochromator with an incident wavelength of 1.49 Å and employed a different single crystal sample aligned in the [h, k, 0] plane. Artifacts such as multiple scattering or higher-order contamination can produce spurious forbidden peaks in single crystal diffraction. However, in this WAND experiment, we observed peaks cor-responding to the (2h + 1, 0, 0) and (0, 2k + 1, 0) families (Fig. 17), extending to higher h and k, effectively ruling out accidental multiple reflections. Since forbidden peaks were observed in both experiments, which used different incident wavelengths and different single crystal samples, we conclude that the forbidden peak arises from intrinsic structural changes, not experimental artifacts.