Quantifying uncertainties in crystal electric field Hamiltonian fits to neutron data

We systematically examine uncertainties from fitting rare earth single-ion crystal electric field (CEF) Hamiltonians to inelastic neutron scattering data. Using pyrochlore and delafossite structures as test cases, we find that uncertainty in CEF parameters can be large despite visually excellent fits. These results show Yb$^{3+}$ compounds have particularly large $g$-tensor uncertainty because of the few available peaks. In such cases, additional constraints are necessary for meaningful fits.

ion magnet stability [6,7], and the exchange interactions between ions [8]. For rare earth magnetic ions, where magnetic anisotropy is strong, a common way to experimentally measure magnetic anisotropy is by fitting the crystal electric field (CEF) Hamiltonian to measured CEF excited levels. Often, this is done with neutron scattering, where the low-energy excited levels are clearly resolved [9]. However, fits to CEF neutron scattering peaks can sometimes be underdetermined (c.f. Yb 2 Ge 2 O 7 [10,11]) and CEF-derived anisotropy does not always match bulk measures of anisotropy (c.f. YbCl 3 [12]). Because neutron CEF studies rarely report uncertainties for the fitted CEF parameters, it is unclear how serious the discrepancies are. If CEF-derived quantities are to be useful for other studies (for example, using the gtensor to fit exchange constants), the error bars for the fitted quantities must be accurately known.
In this study, we propose a method for quantifying uncertainties of CEF fits by using a stochastic search method to map out the χ 2 contour. We test the method by fitting to simulated neutron scattering data for various rare earth ions. We find that the g-tensor uncertainties are strongly ion-dependent, with Yb 3+ often having extremely large uncertainty. These results not only demonstrate a method for rigorously defining error bars on CEF fits, but also reveal which ions are most in need of additional constraints and which quantities are most susceptible to error when fitting CEF levels.

Method
The CEF Hamiltonian can be written as: where O m n are the Stevens Operators [13,14] and B m n are scalar CEF parameters. At a single wavevector Q, the neutron cross section of a crystal field Hamiltonian is written where A is a normalization factor, p n is the Boltzmann weight, and | Γ m |Ĵ ⊥ |Γ n | 2 is computed from the inner product of the matrix element of magnetic moment with the CEF eigenstates |Γ n . In general, because of magnetic form factors and potentially anisotropic g-tensors, this leads to a wavevector dependence of the intensity. However, it is common practice to fit to a constant wavevector Q, allowing for the Q-dependent terms to be neglected. In a neutron experiment, the CEF parameters B m n are fitted to the observed intensities in Eq. 2. To explain the method by which we determine the CEF uncertainties, we consider the example of Yb 2 Ti 2 O 7 .
Example: Pyrochlore Yb 2 Ti 2 O 7 Yb 2 Ti 2 O 7 is a pyrochlore material with magnetic Yb 3+ ions in a D 3 scalenohedron ligand environment, with a three-fold rotation axis along the local [111] direction shown in Fig.  1(a) [15]. In the Stevens Operator formalism [13,14], the D 3 symmetry gives six symmetryallowed CEF parameters:  best fit = .
200 K (c) show the point-charge model simulated scattering at 10 K and 200 K, respectively. The three curves show three fits to the simulated scattering data, offset along the y axis for clarity. The bottom (grey) shows the original model and best fit, the middle (green) shows the minimum g zz to within uncertainty, and the top (red) shows the maximum g zz to within uncertainty.
Using PyCrystalField software [16], we simulated CEF Hamiltonian using a point charge model based on the structure reported in Ref. [17] and 10 nearest oxygen ions. We calculated the neutron spectra at T = 10 K and T = 200 K to simulate intensities at realistic experimental temperatures. To simulate counting statistics of a real neutron experiment, we added intensity-dependent noise to the simulated data based on the Poisson counting statistics of neutron experiments, plus an intensity-independent Gaussian background noise. This method allows us to precisely define the error bar of each simulated data point. For the peak widths, we use a linear energy dependent Gaussian resolution function to define the Gaussian widths of the peaks, plus an energy independent Lorentzian broadening contribution which varies with temperature to account for finite-lifetimes at nonzero temperatures. These two broadning contributions were simulated with a Voigt profile for computational efficiency. This gave a realistic simulated neutron scattering data where the "correct" CEF Hamiltonian is exactly known.
After generating this data set, we defined a global χ 2 fit function based on nine fitted parameters: the six CEF parameters, an overall scale factor, and the two Lorentzian broadening factors. The Gaussian width resolution function was fixed to the simulated values, and thus treated as precisely known. The T = 10 K simulation shows three peaks with three intensities, giving a total of six observable quantities related to the CEF Hamiltonian. Including a second higher temperature reveals two additional transitions: |Γ 1 → |Γ 2 and |Γ 1 → |Γ 3 , bringing the total observable quantities to eight (the energies of these transitions are determined by the low-temperature peak energies, and thus contain no new information). Adding to this the two Lorentzian broadening parameters (which are not related to the CEF Hamiltonian), the total independent observed quantities in Fig. 5 comes to ten. Fitting nine parameters to ten observable quantities is a fully constrained fit. The model which best fits this simulated data, obviously, uses the nine parameters used to generate the data. This has a reduced χ 2 ∆χ 2 red = 1 of the best fit can be considered a valid solution to within one standard deviation uncertainty [18]. Thus, to calculate uncertainties in the CEF Hamiltonian, we must determine the allowed variation of the nine parameters such that ∆χ 2 < 1.
To calculate this, we use a two-step method: incremental search, and then Monte Carlo search. First we select a fitted parameter, fix it to a slightly increased value from the optimum fit, and fit the remaining eight parameters. If the fitted solution is less than ∆χ 2 red = 1, we save the solution as valid and increase the fixed parameter again, using the last fit for starting parameters. If the new solution has greater than ∆χ 2 red = 1 from the original optimum, we return to the optimum solution and repeat the process decreasing the fixed value from the optimum. We then repeat the process for each CEF parameter. As a second step, after looping through each CEF parameter, we then employ a series of Monte Carlo Markov Chains using each valid solution as a starting point, keeping all solutions within ∆χ 2 red = 1. In this way, we effectively map out the allowed variations in each parameter. The distributions of various χ 2 red solution for Yb 2 Ti 2 O 7 are shown in Fig. 2. This family of ∆χ 2 < 1 solutions in Fig. 2 reveals the uncertainties in both the CEF parameters and the CEF derived quantities. The CEF parameter uncertainties are straightforward, defined by the range of parameter fit values. For derived quantities like the ground state eigenkets or the g tensor, we calculate these quantities for each solution and then take the range of calculated values to be the uncertainty bounds. In this way the uncertainties are propagated through the CEF calculations.
The resulting CEF parameters with uncertainty are in Table 1, the eigenvectors and eigenvalues with uncertainty are in Table 2, and the g tensor is g xx = g yy = 4.10 +0.14 −0.15 , g zz = 2.05 +0.06 −1.3 . Several things are worth noting about the Yb 2 Ti 2 O 7 CEF uncertainty calculations. First, some quantities-like g zz -can vary quite a lot even though neutron scattering signal barely changes. To illustrate this, the maximal and minimal g zz models are plotted in Fig. 1(b)-(c). All these solutions would be considered "good fits" to the data (the fitted energy eigenvalues in Table 2 have tiny uncertainties), but g zz = 0.7 is far from the true value g zz = 2.0. Second, the uncertainties in both the fitted CEF values and the resulting quantities can be highly asymmetric, evidenced both in Fig. 2 and the g-tensor variation.

Pyrochlores
The substantial variation in CEF solutions for Yb 2 Ti 2 O 7 is not so surprising given that only three peaks are visible in the low-temperature neutron spectrum. Such a fit is poorly constrained. Other ions with larger effective J values would have more visible peaks, and would thus fare much better. To test this, we repeated the above method but replaced the Yb 3+ ion in Yb 2 Ti 2 O 7 point charge model with other rare earth ions: Sm 3+ , Nd 3+ , Ce 3+ , Dy 3+ , Ho 3+ , Tm 3+ , Pr 3+ , Er 3+ , and Tb 3+ . The ligand environment is exactly the same for each fit-the only thing that changes is the magnetic ion. (Note that not all these materials exist as cubic pyrochlores; the point here is to compare the relative uncertainties for different ions in identical ligand environments.) The uncertainties in the ground state eigenkets and g tensor values from the χ 2 red contours are shown in Table 3. As expected, most other ions have smaller uncertainty in the ground state CEF wavefunction than Yb 3+ . The presence of more CEF levels constrains the fit much better. (One exception to this is Pm 3+ , not listed in Table 3: this ion only gives two visible peaks in its neutron spectrum, and the range of possible solutions is so great that the uncertainty was func- replaced with other rare earth ions. Only the three largest contributions to the ground state eigenket are listed. Sm 3+ and Ce 3+ both have a uniquely defined ground state constrained by symmetry despite variation in CEF parameters, but the rest allow for variation in the ground state wavefunction. Of all the ions, Yb 3+ and Dy 3+ have the largest g tensor uncertainty. Note that many ions listed are non-Kramers and are not in general required to have a doublydegenerate ground state, but do because of the pyrochlore lattice symmetry. 200 K (c) Figure 3: Delafossite NaYbSe 2 simulated scattering and fits. (a) shows the crystal structure of NaErSe 2 , the basis for this fit, with a three-fold axis along the c axis, which we set as z.
(b) and (c) show the point-charge model simulated scattering with Yb 3+ as the central ion at 10 K and 200 K, respectively. The three curves show three fits: the bottom (grey) shows the original model and best fit, the middle (green) shows the minimum g zz to within ∆χ 2 < 1, and the top (red) shows the maximum g zz to within ∆χ 2 < 1. tionally infinite.) For most ions, the ground state wavefunction is generally well-constrained by a CEF fit to neutron data.
Two unusual cases here are Sm 3+ and Ce 3+ , where the ground state wavefunction is exactly defined. Because of the D 3 symmetry of the RE 2 Ti 2 O 7 site, one eigenstate doublet is constrained to be | ± 3/2 exactly (this is the second Yb 3+ excited state in table 2). For Sm 3+ and Ce 3+ in the Yb 2 Ti 2 O 7 structure, this ket is the lowest energy state. Thus, even though there is substantial variation in the B m n values, the ground state anisotropy is precisely known. Thus a large uncertainty in the CEF parameters does not necessarily lead to a large uncertainty in the magnetic ground state.

Delafossites
To test whether these large error bars extend beyond pyrochlores, we now consider uncertainties in delafossite structures using the same method. The delafossite AReB 2 structure also has D 3 symmetry for the magnetic site, and for this series we based the point charge model on the NaErSe 2 chemical structure [19]. Thus there is the same number of fitted parameters as in the pyrochlores. The results of the fits are listed in Table 4. The simulated data and best fits for NaYbSe 2 are shown in Fig. 3.
Despite the different number of ligands and different environment, the results for delafossites are similar to pyrochlores: most ions have well-constrained uncertainties in the ground state anisotropy except for Yb 3+ . There are however some exceptions: Ho 3+ fares poorly in the delafossite g zz uncertainty, while Dy 3+ fared worse in the pyrochlore g zz uncertainty.
The uncertainty in the NaYbSe 2 fitted CEF Hamiltonian becomes more interesting when we plot χ 2 red vs the fitted values in Fig. 4. Here there are two local minima with almost identical χ 2 red . The "real" solution has χ 2 red = 0.9682 and g zz = 1.316, g xx = 3.086 (easy plane anisotropy). The alternate solution has χ 2 red = 0.9702 and g zz = 2.626, g xx = 2.600 (isotropic).  This double-solution problem was encountered experimentally in NaYbO 2 [20]; to distinguish these two solutions would be impossible with neutron scattering data alone. In a case such as this, it is important to (i) fully map out the χ 2 contour to identify competing solutions, and (ii) collect additional data or information to identify the correct anisotropy [21,22]. As a side note, Tables 3 and 4 show candidates for strong quantum effects in pyrochlores and delafossites. The g xx values are directly related to J ± expectation values, which are rough measures of quantum tunneling between the ground states. For pyrochlores, Yb 3+ dominates because of its large weight on | ± 1/2 . For the delafossites meanwhile, there are many promising candidates, most notably Nd 3+ , Ce 3+ , and Yb 3+ with their large | ± 1/2 weights. If these point charge calculations are at least approximately close to the real material Hamiltonians, these results give direction on where to find strongly quantum delafossite materials.

Bixbyite Yb 2 O 3
In D 3 symmetry, the Yb 3+ ions appear to have the largest CEF uncertainties. This problem will in principle get worse as the number of crystal field parameters increases in lower symmetry structures. As an example, we considered the first Yb 3+ site in Bixbyite Yb 2 O 3 . This material has two symmetry inequivalent Yb 3+ sites, but we consider the first site which has C 3 symmetry: a three-fold axis about [111] but no mirror planes. This symmetry allows for nine CEF parameters: 3 6 , and B 6 6 . To estimate uncertainty, we follow the same procedure outlined above. The simulated data and best fits are shown in Fig. 5, the best fit CEF parameters are in Table 5, and the calculated g tensor is g xx = 1.4 +1.5 −0.9 , g zz = 7.0 +0.4  200 K (c) The three curves show three fits: the bottom (grey) shows the original model and best fit, the middle (green) shows the minimum g zz fit, and the top (red) shows the maximum g zz fit.  axis-opposite of the pyrochlore and delafossite-but as expected, the uncertainties are even larger. Indeed, the uncertanties of the CEF parameters in Table 5 are so large that most of the CEF parameters are zero to within uncertainty-hardly useful for any detailed modeling. This is not surprising given the number of independent parameters: in this case, we are fitting 12 parameters (nine CEF parameters, a scale factor, and two Lorentzian widths) to 10 observed quantities (three transition energies, five transition intensities, and two Lorentizan widths), which technically is an underdetermined problem. That we are able to perform a fit at all is evidence that the fitted parameters are not totally independent, allowing for finite uncertainties. Thus, a Yb 3+ CEF model with nine independent parameters definitely needs more information than just neutron scattering peaks to constrain a CEF fit.

Discussion and Conclusions
The main conclusion of this numerical study is that uncertainties in CEF fits can be very large even though the fits may look visually good. This is important because these uncertanties should be propagated through other calculations that use the g tensor (such as field-polarized spin wave calculations).
Furthermore, it should be noted that fits to real experimental data will probably have larger uncertainty than the idealized fits we perform here. In real experiments, there are background contributions from phonons and the sample environment, the resolution function is often not precisely known, CEF-phonon coupling affects measured intensities, and peak shapes may be asymmetric. All of these will worsen the agreement between the model and the data. Thus the true uncertainties may be much larger than those we estimate here.
This problem is particularly severe for Yb 3+ compounds as they only have three excited levels. This is unfortunate, as Yb 3+ receives much attention as an effective J = 1/2 host for quantum magnetism. In such cases it is necessary to include additional experimental information, like electron spin resonance as was done for NaYbS 2 [23], nonlinear susceptibility and high-field torsion magnetometry as was done for CsYbSe 2 [24], or saturation magnetization as was done for YbMgGaO 4 [25].
This being said, a secondary conclusion of this study is that not all calculated quantities are affected equally by CEF parameter uncertainties. The clearest examples of this are Sm 2 Ti 2 O 7 and Ce 2 Ti 2 O 7 , where the ground state is precisely known despite uncertainty in the CEF model. This is also true of Yb 3+ : although the pyrochlore and delafossite fits show large uncertainty for g zz , the g xx is more constrained. Likewise, the Yb 3+ ground state eigenkets have relatively modest uncertainties. Therefore even if the overall anisotropy might be in question, the fitted CEF model might still give accurate and useful information about the ground state wavefunction.
In summary, we have shown by simulating and fitting to artificial CEF neutron scattering data sets that CEF fits can have very large uncertainties. In three-fold symmetric environ-ments, Yb 3+ consistently has the largest uncertainties, highlighting the need for additional constraints when fitting its CEF levels. However, the uncertainties in calculated quantities are highly dependent upon the details of the model-some quantities are well-constrained despite uncertainty in the fitted parameters. In all cases, it is important to explore the full χ 2 contour of a CEF model so that uncertainty can be known.