# CrysFieldExplorer: a software for rapid optimization of crystal field
Hamiltonian

Qianli Ma  Neutron Scattering Division, Oak Ridge National Laboratory, Oak
Ridge, Tennessee, USA,37831  Xiaojian Bai  Neutron Scattering Division, Oak
Ridge National Laboratory, Oak Ridge, Tennessee, USA,37831  Department of
Physics and Astronomy, Louisiana State University , Baton Rouge, Louisiana,
USA,70803  Erxi Feng  Neutron Scattering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee, USA,37831  Guannan Zhang  Computer Science
and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee,
USA,37831  Huibo Cao  Neutron Scattering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee, USA,37831

###### Abstract

We present a new lite python-based program, CrysFieldExplorer, for fast
optimizing crystal electric field (CEF) parameters to fit experimental data.
The main novelty of CrysFieldExplorer is the development of a unique loss
function, referred to as the Spectrum-Characteristic Loss
($L_{\text{Spectrum}}$), which is defined based on the characteristic
polynomial of the Hamiltonian matrix. Particle Swarm Optimization and
Covariance matrix adaptation evolution strategy are used to find the minimum
of the total loss function. We demonstrate that CrysFieldExplorer can performs
direct fitting of CEF parameters to any experimental data such as neutron
spectrum, susceptibility, magnetizations etc. CrysFieldExplorer can handle a
large amount of none-zero CEF parameters and reveal multiple local and global
minimum solutions. Detailed crystal field theory, description of the loss
function, implementation and limit of the program are discussed within context
of two examples.

##  I Introduction

Single-ion magnetic anisotropy is one of the key elements for exotic quantum
states in frustrated systems Gingras  _et al._ (2000); Gardner  _et al._
(2010); Cao  _et al._ (2009a); Bertin  _et al._ (2012); Gaudet  _et al._
(2015); Dun  _et al._ (2021); Paddison  _et al._ (2017); Gaudet  _et al._
(2019). Crystal electric field (CEF) is responsible for single-ion magnetic
anisotropy, and it occurs through Coulomb interaction and Pauli exclusion
between the central cation and surrounding anions that splits the energy
levels of electrons grouped by different orbitals of the central cation. In
the absence of the CEF effect or in a spherically symmetrical field, these
energy levels would otherwise be degenerate. The CEF effect is known to cause
dramatic magnetic anisotropy which collectively produce exotic ground states
in a wide range of quantum materials Bramwell and Gingras (2001); Ross  _et
al._ (2011a); Gaudet  _et al._ (2019); Porée  _et al._ (2022); Keimer and
Moore (2017); Smith  _et al._ (2022) and unconventional high temperature
superconductors Fulde  _et al._ (1970); Mesot and Furrer (1997); Lynn  _et
al._ (1980); Metoki  _et al._ (2004); Goremychkin  _et al._ (2004); Rubio
Temprano  _et al._ (2000); Boothroyd  _et al._ (1992); Pickett (2021).

The calculation of excited energy levels of magnetic ions in a crystalline
environment has been well studied Stevens (1952); Hutchings (1964); Lea  _et
al._ (1962); Walter (1984); Loewenhaupt (1990). Stevens and Hutchings have
illustrated systematically the process of determination of the perturbing
Hamiltonian from the evaluation of the electrostatic potential experienced by
the magnetic ion from the surrounding charges Stevens (1952); Hutchings
(1964). One of the common conventions to express the CEF Hamiltonian is using
Stevens Operators. The crux of the CEF study is to solve the single-ion
Hamiltonian by matrix diagonalization and retrieve eigenvalues and
eigenfunctions from the CEF Hamiltonian to fit the corresponding crystal field
excitations. Inelastic neutron scattering is a suitable experimental technique
because it directly measures transitions between ground states and excited
states of magnetic ions. It has gradually become one of the most popular
experimental methods in conducting crystal field analysis. After obtaining
neutron scattering and other bulk property data, the next step in the CEF
analysis is fitting the CEF Hamiltonian to experimental observables. It has
always been a challenging optimization problem to find a solution set that
best describes the experimentally measured data. Efforts to tackle this
problem include programs such as SPECTRE Boothroyd (2014), McPhase Rotter
(2013), Mantid Arnold  _et al._ (2014), SIMPRE Baldovi  _et al._ (2013) CFca
R.Osborn (2017), FOCUS Fabi (1995) and the latest PyCrystalField Scheie
(2021). These programs use various techniques to calculate the CEF parameters
but they all share the same approach that uses $\chi^{2}$, a type of loss
function, to find the global minimum. One of the disadvantages associated with
choosing the $\chi^{2}$ loss function is it relies on a set of starting
parameters that needs to be adjacent to the true solution. As we will
demonstrate in the next chapter, large energy boundaries in the $\chi^{2}$
type loss function exist that could trap the optimization process on a local
minimum. Typically, the point-charge model and Monte Carlo sampling are used
to generate starting parameters. However, the point-charge model is a
classical approximation by positioning several point charges around a magnetic
ion to represent a distribution of valence electrons. This heavily limits the
accuracy of the point-charge model. Additionally, building a point-charge
model or Monte Carlo simulation can sometimes be none-trivial for
inexperienced researchers. More importantly, when dealing with low local
symmetry such as Ci, C1, the number of none-zero crystal field parameters can
be well into high 20s. The complexity and cost of computation of the multi-
dimensional loss function increases exponentially and can easily overwhelm
currently existing software.

To tackle these challenges, we have designed a special loss function, Spectrum
- Characteristic Loss ($L_{\text{Spectrum}}$), based on the theory of
polynomial characteristics and developed a lite python program called
CrysFieldExplorer. CrysFieldExplorer takes advantage of particle swarm
optimization (PSO) Bonyadi and Michalewicz (2017) and covariance matrix
adaptation evolution strategy (CMA-ES) Hansen (2006) to minimize a combination
of Spectrum - Characteristic Loss and traditional $\chi^{2}$ losses. It is
able to quickly fit neutron spectroscopy data and any other experimental
results such as susceptibility, magnetization, specific heat, neutron
diffraction, etc with CEF model and yields a series of solutions containing
information with local and global minima. It bypasses the step of using a
point-charge model to estimate starting parameters. In the following chapters,
we will show our methods can not only discover the solutions reported in
literature, but also uncover multiple other solutions that could fit all
provided experimental data equally well. This suggests that the data used in
the traditional optimization process may not have enough resolution to
distinguish one from another and require further measurements such as local
magnetic susceptibility with polarized neutron diffraction Cao  _et al._
(2009a). Our findings suggest that the traditional standard of determining one
set of best CEF parameters may be flawed. We, therefore, suggest that future
CEF analysis should include a list of possible solutions and discuss their
relative physical meaning and the reasons for the preferred solution.

##  II Theory

###  II.1 Crystal Electric Field Theory

The CEF potential can be constructed from a point-charge perspective detailed
in Hutching’s work Hutchings (1964). The electrostatic potential due to the
surrounding point-charges can be expressed using tesseral harmonic functions
in Cartesian coordinates. Following this convention, the Crystal Field (CF)
Hamiltonian can be expressed by Stevens operators as Stevens (1952):

$H_{CF}=\sum_{n,m}B^{m}_{n}O^{m}_{n}=\sum_{n,m}[A^{m}_{n}\theta_{n}]O^{m}_{n},$
(1)

where $B^{m}_{n}(|m|\leq n)$ is the CF parameters fitted from experimental
measurements, $O^{m}_{n}$ are the Stevens’ Operators. The CF parameters can
also be expressed in terms of $A^{m}_{n}$ and $\theta_{n}$, where $\theta_{n}$
represents reduced matrix elements and is also called as the Stevens factor.
For rare earth ions, $\theta_{n}$ has been tabulated in Table VI of Ref.
Hutchings (1964).

###  II.2 Neutron Scattering Cross-section

Inelastic neutron scattering spectroscopy is well-suited in studying the
crystal field excitations because neutrons can directly excite electron spins
from one level to another and measure the difference between two energy
levels. The observed intensity and excited energy levels can be used to fit
the crystal field parameter $B^{m}_{n}$s in Eq. 1. The excited energy levels
correspond to the energy difference between the lowest eigenvalues and the
corresponding state of the CEF Hamiltonian. The observed intensity is related
to the partial differential magnetic cross-section expressed as:

$\frac{d^{2}\sigma}{d\Omega
dE^{\prime}}=C\dfrac{k_{f}}{k_{i}}F(|Q|)S(|Q|,\hbar\omega),$ (2)

where $k_{f}$ and $k_{i}$ are the momentum of scattered and incident neutrons,
C is a constant including the Debye-Waller factor, $F(|Q|)$ is the magnetic
form factor of the sample. $S(|Q|,\hbar\omega)$ is the scattering function,
$\hbar\omega$ indicate the energy of neutrons. From $S(|Q|,\hbar\omega)$ the
relative intensities from different states of the CEF exictations can be
calculated. At constant $|Q|$, the $S(|Q|,\hbar\omega)$ is,

$S(Q,\hbar\omega)=\sum_{i,i^{\prime}}\dfrac{(\sum_{\alpha}|\langle
i|J_{\alpha}|i^{\prime}\rangle|^{2})e^{-\beta E_{i}}}{\sum_{j}e^{-\beta
E_{j}}}L(\Delta E+\hbar\omega),$ (3)

$i\rightarrow i^{\prime}$ indicates the transition between state $i$ and
$i^{\prime}$ for the magnetic ion. $\beta=-\dfrac{1}{k_{B}T}$, $k_{B}$ is the
Boltzman constant, T indicates temperature. $L(\Delta E+\hbar\omega)$ is a
Lorentzian function that guarantees the energy conservation when neutrons
induce a transition from state $i$ to $i^{\prime}$, which poses a finite
energy width or lifetime. The $\sum_{\alpha}$ sums all three $x,y,z$ component
of the $J_{\alpha}$ matrix. Although many CEF experiments are measured at low
temperature, high temperature measurements can be important for confirming or
extracting the scattering signal originated from CEF excitations that can be
calculated with Eq. 3.

###  II.3 Magnetization and Susceptibility

Magnetization under an external magnetic field $\mathbf{H}$ can be readily
calculated by combining CEF Hamiltonian and the Zeeman term. The overall
Hamiltonian can be expressed as:

$H=H_{CF}-\mu_{B}g_{J}\mathbb{H}\cdot\mathbb{J}.$ (4)

By diagonalizing Eq.4, we can calculate the eigenstates $E_{n}$ and
eigenfunctions $|i\rangle$ of the CF Hamiltonian in a magnetic field
$\mathbf{H}$. The three components ($\alpha=x,\ y,\ z$) of the magnetization
in Cartesian coordinate system are given as,

$M_{\alpha}(\mathbf{H},T)=g_{j}\sum_{n}e^{-\beta E_{n}}\langle n|J|n\rangle/Z$
(5)

where $Z=\sum_{n}e^{-\beta{E_{n}}}$ is the partition function. The powder
average of the magnetization can be derived from Eq. 5 by calculating the
averaged magnetization on a unit sphere.

The magnetic susceptibility can be calculated by taking derivatives of the
magnetization $M_{\alpha}$,

$\chi_{\alpha\beta}=\dfrac{\partial M_{\alpha}}{\partial H_{\beta}}.$ (6)

The powder averaged magnetic susceptibility can be calculated in the same
fashion by substituting $M_{\alpha}$ with its powder-average version. Detailed
information about the calculation can be found in ref.Dun  _et al._ (2021).

###  II.4 Loss Function

Figure 1: Er3Mg2Sb3O14 comparison of Spectrum - Characteristic Loss and the
mean square root loss functions commonly used in other optimization programs.
Both losses are calculated along the line in a 15-dimensional parameter space
and projected on the B${}^{0}_{2}$ dimension.

The core of an optimization problem is constructing a proper loss function
that is smooth on the parameter space, sensitive to the change of input
parameters and does not produce large energy boundaries around local and
global minima. Currently, a few popular softwares such as SPECTRE Boothroyd
(2014), PyCrystalField Scheie (2021) and Mantid Arnold  _et al._ (2014), etc
are able to fit the CEF Hamiltonian using experimental observables. Different
to the existing software packages, here we introduce a newly designed loss
function, Spectrum - Characteristic Loss. and demonstrate its advantage when
deployed in CEF optimization.

The foundation of the Spectrum - Characteristic Loss is from the theory of
characteristic polynomial - $det(\lambda I-A)=0$, where $\lambda$ is the
eigenvalue of matrix A, $I$ is an identity matrix. The difference between the
lowest and excited eigenvalues are neutron observed energy levels. Using this
feature, we construct the loss function for the neutron measured energy levels
as:

$L_{E}=\log_{10}\left(\sum_{i=1}\dfrac{det\\{(E_{\text{exp}}[i]+E_{\text{cal}}[0])\mathrm{I}-H\\}^{2}}{det\\{E_{\text{exp}}[i]I\\}^{2}}\right),$
(7)

$E_{\text{exp}}$ indicate the observed excitation levels from neutron
scattering, the summation of $i$ starts from 1 which is the first excited
energy level observed by neutron experiment. $E_{\text{cal}}[0]$ is the ground
state eigenvalue from matrix $H$. $E_{\text{cal}}[0]$ may not be 0.

Subsequently, we construct the mean square root deviation for the intensity
and the inverse susceptibility as:

$L_{\text{Intensity}}=\dfrac{\sqrt{\sum({I_{\text{true}}[i]-I_{\text{calc}}}[i])^{2}}}{\sqrt{\sum{(I_{\text{true}}}[i])^{2}}},$
(8)

and

$L_{1/\chi}=\dfrac{\sqrt{\sum({1/\chi_{\text{true}}[i]-1/\chi_{\text{calc}}}[i])^{2}}}{\sqrt{\sum{(1/\chi_{\text{true}}}[i])^{2}}}.$
(9)

To construct the Spectrum - Characteristic loss function, we combine Eq.7, 8
such that $L_{\text{Spectrum}}$=$log_{10}(L_{E})+L_{\text{Intensity}}$. To fit
other bulk measured physical properties, CrysFieldExplorer sums up
$L_{\text{Spectrum}}$ with other mean square root losses computed from the
experimental data and conduct global optimization.

Figure 1 compares the Spectrum - Characteristic Loss with traditional mean
square losses of the excited energy levels and relative intensities for an
example material Er3Mg2Sb3O14 . It represents an optimization conducted on a
15-dimensional parameter space. More physical details about the fitting
results are discussed in a later chapter. Here we show the computed loss along
a straight line passing through the true solution of Er3Mg2Sb3O14 Dun  _et
al._ (2021). The line is chosen to point in a random direction within the
15-dimensional parameter space to simulate an optimization procedure (a real
optimization may not follow a straight line). The loss is then projected onto
the B${}^{0}_{2}$ dimension for visualization. The upper panel shows the
Spectrum - Characteristic Loss function while the lower panel is the standard
$\chi^{2}$ loss function. To compare the loss functions like-to like we also
apply a $log10$ operation on the energy levels such that
$\chi^{2}=log_{10}(\chi^{2}_{E})+\chi^{2}_{I}$. Comparing the upper panel to
the lower panel, one can tell that the Spectrum - Characteristic Loss has a
global structure around the minimum and does not produce sharp energy barrier
at B${}^{0}_{2}\sim$ -0.02. Such an energy barrier can trick the optimization
algorithm assuming it has already reached the global minimum. This can lead to
a wrong solution if the starting parameters fall within the energy well. The
Spectrum - Characteristic Loss also appears to be more sensitive to slight
perturbations in CEF parameters. In the specific direction that the loss is
computed, traditional $\chi^{2}$ loss function is only able to converge at the
true solution which is artificially designed such that the energy levels
exactly match between calculated and observed, causing both
$L_{\text{Spectrum}}$ and $\chi^{2}$ to appear sharp at the true solution. On
the other hand, the Spectrum - Characteristic Loss function can discover at
least three other minima in the parameter space. In reality the real global
minimum is unknown for a measured material, so it is important to have a
complete list including all potential minima for further examination as
different CEF parameters from the list could all reproduce results identical
to experimentally measured properties. As the number of CEF parameters
increases, it requires more constraints to conclusively solve the CEF
Hamiltonian. The multiple solutions CrysFieldExplorer discovered is an
indication that the input observables do not impose enough constraints to
reach the true solution.

In addition, CrysFieldExplorer can add any other experimental data such as
specific heat, anisotropy g-tensors, magnetization, etc in the optimization
process. The loss function for additional data can be added the same way as
mean square root deviation. Users can also adjust weights for different losses
to fine tune the optimization process.

##  III Examples

###  III.1 Benchmark with Rare Earth Pyrochlore: $Yb_{2}Ti_{2}O_{7}$

Geometrical frustration has been of great interest in condensed matter physics
due to the ability to host rich order-disorder states. These exotic states of
matter are results of collective behaviors that arise from frustrated
interactions within the quantum many-body system. Understanding these
phenomena is crucial for developing next generation quantum technologies
Tokura  _et al._ (2017); Ball (2017); Giustino  _et al._ (2021); Cava  _et
al._ (2021); Bassman  _et al._ (2020); Lau  _et al._ (2020); Stanev  _et al._
(2021). Rare earth pyrochlore is the archetype of magnetic frustration in
three dimensions. This system has been reported to host exotic magnetic states
such as quantum spin ice Ross  _et al._ (2011a); Kimura  _et al._ (2013);
Martin  _et al._ (2017), quantum spin liquid Gingras and McClarty (2014); Ross
_et al._ (2014); Smith  _et al._ (2022); Bhardwaj  _et al._ (2022) etc. The
pyrochlore structure generally possesses chemical composition of A2B2O7 (B=
Sn, Ge, Pt, Zr, Ti etc) with a space group $Fd\bar{3}m$ and point group
$D_{3d}$ at A site. We selected this system as a high symmetry end to
benchmark our optimization model. Its local site possesses a 3-fold rotation
axis along the local [111] direction. Due to the symmetry constraints, the
number of none zero CEF parameters is 6. In this chapter, we demonstrate the
use of CrysFieldExplorer by comparing a list of optimized results with
observables in typical measurements such as inelastic neutron scattering,
susceptibility, magnetizations and anisotropy g-tensors.

Figure 2: The corner-sharing tetrahedral structure of Yb2Ti2O7. The blue
spheres indicate Yb3+ ions and the red spheres are the surrounding oxygen
atoms around one Yb3+ site. The z-axis is the three-fold rotation axis along
the $<$111$>$ direction.

One of the representative systems is the Yb2Ti2O7, a promising quantum spin
ice candidate with XY-type single-ion magnetic anisotropy Cao  _et al._
(2009b, a); Ross  _et al._ (2011b); Gaudet  _et al._ (2015). The ground states
and crystal field effect of Yb3+ have been well studied by Gaudet et al. in
Ref.Gaudet  _et al._ (2015). Yb3+ has an electronic configuration of
4$f^{13}$, total angular momentum $\mathbf{\textit{J}}=7/2$. The crystal
electric field effect originates from the surrounding O2- ions of the Yb3+
site. As a result, $2J+1=8$ \- fold degeneracy is expected to be lifted by the
CEF effect from the ground state, resulting in 7 excited levels. Additionally,
Yb3+ is a Kramers’ ion, therefore the 8 - fold degeneracy are all doubly
degenerate into 4 well separated doublets. A total of 3 excitations are
expected to be observed through inelastic neutron scattering spectroscopy.
Oxygen atoms around the local A-site form a distorted cubic structure shown in
Figure 2. The most convenient placement of the local coordination system
placing the local $z$ \- axis along the local three-fold $<$111$>$ rotational
axis perpendicular to the oxygen plane. The resulting Hamiltonian takes the
form of:

$H_{CEF}=B^{0}_{2}\hat{O}^{0}_{2}+B^{0}_{4}\hat{O}^{0}_{4}+B^{3}_{4}\hat{O}^{3}_{4}+B^{0}_{6}\hat{O}^{0}_{6}+B^{3}_{6}\hat{O}^{3}_{6}+B^{6}_{6}\hat{O}^{6}_{6}.$
(10)

Eq. 10 follows the convention of Stevens operators for $\hat{O}^{m}_{n}$,
$B^{m}_{n}$s are CF parameters used to describe the Coulomb potential
generated by the surrounding oxygen atoms. Inelastic neutron scattering
results have been reported in Ref.Gaudet  _et al._ (2015).

B${}^{m}_{n}$ (meV) | A.BertinBertin  _et al._ (2012) | J.GaudetGaudet  _et al._ (2015) | This work 111The parameters indicated by green circles in Figure 1 | This work 222The parameters indicated by black circles in Figure 1
---|---|---|---|---
B${}^{0}_{2}$ | 1.270 | 1.135 | 1.57 | -2.22
B${}^{0}_{4}$ | -0.0372 | -0.0615 | -0.0365 | 0.00241
B${}^{3}_{4}$ | 0.275 | 0.315 | -0.627 | -0.0612
B${}^{0}_{6}$ | 0.00025 | 0.0011 | 0.00184 | -0.00102
B${}^{3}_{6}$ | 0.0023 | -0.037 | -0.0285 | 0.0287
B${}^{6}_{6}$ | 0.0024 | 0.005 | -0.0211 | -0.0368
gx(y) | 4.09(2) | 3.69(0.15) | 3.72 | 2.43
gz | 2.04(3) | 1.92(0.2) | 2.07 | 4.51
Table 1: The calculated and fitted CEF parameters from Ref.Gaudet  _et al._
(2015). The calculated values are obtained from point-charge model while the
fitted values are from fitting inelastic neutron scattering data. As can be
seen the difference between the fitted and calculated B${}_{6}^{0}$ and
B${}_{6}^{3}$ is off by a magnitude.

Table 1 lists previous works on the CEF analysis for Yb2Ti2O7 . Bertin et
al.’s results Bertin  _et al._ (2012) are obtained within the point-charge
approximation and taken as starting parameters for the work by Gaudet etl al
Gaudet  _et al._ (2015). They deployed a typical mean-square-root deviation
minimization algorithm that searched through the six dimensional parameter
space in the vicinity of the starting values from Ref. Bertin  _et al._
(2012). By using a point-charge model to generate starting parameters, one
assumes these starting values are close to the global minimum. However, Table
1 shows that the refined value from point-charge model can be different by an
order of magnitude (i.e, B${}^{0}_{6}$, B${}^{3}_{6}$).

Now we will demonstrate the results of CrysFieldExplorer in the Yb2Ti2O7
system. We adapt the searching algorithm using particle swarm optimization
with our customized Spectrum - Characteristic Loss function discussed in
previous sections. We construct the PSO model with a particle size N=400 and
iteration number of 100 times. The details of the hyper-parameters can be
found in Ref. Bonyadi and Michalewicz (2017). The total loss is constructed as
$L_{total}=L_{\text{Spectrum}}+L{1/\chi}$, where $\chi$ in this case
represents the susceptibility data. These are the same experimental data used
to obtain results in Ref.Gaudet  _et al._ (2015). We present
CrysFieldExplorer’s results in Figure 3 and compare it with those reported in
Ref.Gaudet  _et al._ (2015).

Figure 3 lists all 150 results from CrysFieldExplorer on Yb2Ti2O7 in a log
scale along the y-axis. Each Yb2Ti2O7 CEF parameter is plotted against a
custom defined goodness of fit $\chi^{2}=\dfrac{[Expected\ Value-Observed\
Value]^{2}}{Expected\ Value^{2}}$. Readers should be aware of the definition
of present $\chi^{2}$ is simply the sum of percentage difference between true
solution and calculated solution. It is calculated by comparing the excited
energy levels and relative intensities, neutron spectroscopy, susceptibility
and magnetization using CEF parameters provided in the Ref. Gaudet  _et al._
(2015). This $\chi^{2}$ is merely an indication of the agreement goodness
between calculated and true observables. We chose $\chi^{2}<1$ as a cut off
number in present study as we find results below this number shows good
overall agreement between fit and true solution. The solid grey circles
represent 150 solutions generated using PSO algorithm. The solid red circles
represent 13 solutions with $\chi^{2}<1$. The blue circles mark the parameters
reported from Ref. Gaudet  _et al._ (2015) and is taken as the ”true solution”
in this example. The $\chi^{2}$ for the true solution is given an artifically
0 but shifted by 0.03 so it is visible on the log scale. One of the important
observations from Figure 3 is that there exists three local minima in the
parameter space of Yb2Ti2O7 crystal field. Observing Panel a), b), d), f) from
Figure 3 shows two out of the three local minima are concentrated regions
where CrysFieldExplorer can easily converge. There is a third region that
contains one acceptable solution with $\chi^{2}<1$. The true solution reported
from Ref. Gaudet  _et al._ (2015) marked by blue circle falls in 1 of the
minima. The green and black circles are in two other minima. Interestingly,
while no constraints were set on B${}_{4}^{3}$ and B${}_{6}^{3}$, panels c)
and e) show symmetric behavior w.r.t to 0 as expected by the symmetry of local
crystalline environment.

To analyze the physical meaning of the CEF parameters found by
CrysFieldExplorer, the main panels in Fig. 4 a)-d) show the comparison between
calculated and true results in terms of inelastic neutron scattering data,
inverse susceptibility, magnetization and the ratio of gz/gx(y) for the Yb3+.
Unfortunately the experimentally measured susceptibility and magnetization
data was not given in Ref. Gaudet  _et al._ (2015). We assume the calculated
susceptibility and magnetization using the CEF parameters provided from Ref.
Gaudet  _et al._ (2015) as ”true solution”. Each red lines in a)-c) and the
red data points in d) indicate the calculated properties. The blue lines and
open blue circles in a) - d) represents the true solutions. The black, green
color scheme indicate two other minima results from present study
corresponding to the same green and black data points in Figure 3. We also
list our representative results (green and black data) as CEF parameters in
Table 1. They are marked by a and b respectively. Clearly, Figure 3 and 4
demonstrate our PSO optimizer can successfully converge to the true CEF
parameters region from Ref.Gaudet  _et al._ (2015). A universal starting range
of [-100,100] meV were set for all six B${}^{m}_{n}$s in Eq.1. For general
users, we feel this range can be appropriately determined based on the
magnitudes of excited energy levels. Additionally, as Spectrum -
Characteristic Loss is smooth and can effectively avoid producing steep energy
barriers seen in Fig. 1, a larger range for the CEF parameters should not be
problematic. The average time for the PSO to converge to one solution in the
6-dimensional CEF parameter space was approximately 4.7 minutes on a 12th Gen
Intel(R) Core i7-1265U laptop with single core processing. The maximum
magnitude of the searching range compared to the order of the minimum CEF
parameter B${}^{0}_{6}$ is on the order of 104. CrysFieldExplorer can locate
solutions within this vast range of parameter space within minutes. This is an
example to show the powerful optimization capability with CrysFieldExplorer.

Observing each panel in Figure 4 reveals more details of CrysFieldExplorer’s
results. The main panels of a), b) and c) in Fig. 4 compare the calculated
solutions of $\chi^{2}<1$ with the true solution. These solutions are over-
plotted in red and form a region which makes it easier to identify the
relative positions of the calculated solutions in from the true solutions. The
insets of a), b) c) in Figure 4 compare the true solution with the solutions
a) and b) found and listed in Table 1. Solution a) and b) are indicated as
green and black open circle in panel d). Readers should pay particular
attention to solution a) and b) as they are in two different regions in the
CEF parameter space. In panel a) of Fig.4, the calculated spectrum shows good
agreement with the true spectrum except at $\sim$ 76 meV. We estimate
approximately $\pm$ 3$\%$ variation on the peak intensities from our
calculated solutions compared to the true solution. The excitations at 81.8
meV and 116.2 meV matches well in both peak width and integrated intensity. In
panel b), the inverse of susceptibility shows excellent agreement as well. The
CEF fit and true solution are indistinguishable by eye. Panel c) shows the 10K
magnetization data. The solid grey lines mark an region where the true
solution and calculated solutions with $\chi^{2}<1$ reside in. The inset of
panel a), b), c) shows a similar story, no significant difference between true
solution and solutions a) and b). Panel d) is the most interesting plot as it
shows the ratio of the anisotropy g-tensors $\dfrac{g_{z}}{g_{x(y)}}$. Only in
this panel we can distinguish solution b) from solution a) and the true
solution. Observing panel d), most of the calculated anisotropy g-tensors are
in one dominant region where the ratio of $g_{z}/g_{x(y)}$ is $\sim$ 0.5 while
$g_{x}$ is close to 3.8. Our best solution, represented by the green open
circle, yields a set of g-tensor $g_{x(y)}$ = 3.73 $g_{z}$ = 2.07 that are
within the uncertainty given from the published result of 3.69 $\pm$ 0.15 and
1.92 $\pm$ 0.20 Gaudet  _et al._ (2015).The ratio of $g_{z}/g_{x(y)}$ is of
particular interest. Previous polarized neutron studies from Cao et al.Cao
_et al._ (2009b) reported $g_{x}(y)=4.1$ and $g_{z}=2.25$ resulting in a ratio
of 1.82. Bertin’s Bertin  _et al._ (2012) and Gaudet’sGaudet  _et al._ (2015)
results are listed in Table 1 with $g_{z}/g_{x(y)}$= 2.0 and 1.92. From our
solution a), the ratio can be calculated to be 1.79 closes to Cao et al.’s
result. Now we turn our attention to solution b). Despite having a good
agreement with the experimental values in the bulk property measurements, it
shows a strong ratio of $g_{z}/g_{x(y)}$ $\sim$ 3.5, with $g_{z}=4.51$
significantly larger than solution a) and the true solution. This observation
suggests solution b) represents a local minimum that fits well with limited
constraints. By analyzing the anisotropy g-tensor, we conclude that solution
b) does not agree with the physical model of Yb2Ti2O7 .

Now we turn our focus to solution a) which is in a different region with the
true solution, but we cannot distinguish solution a) from the true solution by
all the given measurements. Table 2 compares the ground state wavefunctions
between solution a) marked by green circle in Figure 3 and that reported in
Ref. Gaudet  _et al._ (2015). Both solutions also produce similar ground state
wavefunctions as reflected by the similarity in the anisotropy g-tensors. The
dominant term in the ground state wave function from both solutions are the
$\pm\dfrac{1}{2}$ terms.

Combining these results with previous discussions suggests solution a) is a
degenerate solution under the existing constraints imposed by the experimental
data reported in the Yb2Ti2O7 system. Only one set of CEF parameters were
reported in each paper from previous studies Cao  _et al._ (2009a); Bertin
_et al._ (2012); Gaudet  _et al._ (2015), but our discovery implies that these
papers may need more data to conclusively determine the CEF parameters. The
ground state of Yb2Ti2O7 has been extensively studied from various
perspectives and our results do not imply whether the previous analyses are
incorrect. However, it is an alarming discovery that suggests future CEF
analysis should be cautious in confirming or denying their findings. In next
chapter, we will show that lack of enough constraints may pose serious
uncertainty in determining the true CEF parameters.

$J_{z}$ | $\dfrac{-7}{2}$ | $\dfrac{-5}{2}$ | $\dfrac{-3}{2}$ | $\dfrac{-1}{2}$ | $\dfrac{1}{2}$ | $\dfrac{3}{2}$ | $\dfrac{5}{2}$ | $\dfrac{7}{2}$
---|---|---|---|---|---|---|---|---
GS eigenfunctionsGaudet  _et al._ (2015) | 0 | 0.0866 | 0 | 0 | $-0.9283$ | 0 | 0 | 0.3616
| $-0.3616$ | 0 | 0 | $-0.9283$ | 0 | 0 | $-0.0866$ | 0
GS eigenfunctions333The parameters indicated by green circles in Figure 1 | 0 | $-0.076$ | 0 | 0 | $-0.924$ | 0 | 0 | $-0.375$
| 0.375 | 0 | 0 | $-0.924$ | 0 | 0 | 0.076 | 0
Table 2: The ground state eigenfunctions from ref.Gaudet  _et al._ (2015) and
present work of solution a) as reported in Table 1. Despite the difference in
the CEF parameters, the ground state eigenfunctions looks similar with
dominant $\pm\dfrac{1}{2}$ terms.

Regardless of which set of CEF parameters is correct, one of
CrysFieldExplorer’s biggest advantages is the capability to directly commence
the fitting procedure starting from Eq. 10 without any a priori knowledge of
the starting parameters. This allows users to skip the step of building a
point-charge model or conduct Monte Carlo simulation, either of which can be
daunting for inexperienced users. CrysFieldExplorer only requires a wide
fitting range to be given. All the parameters of Yb2Ti2O7 were chosen from a
range of [-100, 100] meV which can be estimated from the magnitude of the
leading B${}_{2}^{2}$ term. This range requirement can be easily estimated by
users. Too large or too small beyond this range could result in the excited
energy levels to be unreasonable large or small.

To conclude, in this section we demonstrate that CrysFieldExplorer can quickly
produce a list of CEF solutions for Yb2Ti2O7 including the ones reported in
the literature. Our improved loss function $L_{\text{Spectrum}}$ has a higher
sensitivity to the change of CEF parameters and a better global structure that
allows us to search through a wide range of parameter space. Identifying which
one is the real solution is more complicated. On the one hand, there may not
be enough constraints to conclude which solution is the correct one. On the
other hand, as most of the experimental data are collected from powder sample,
the information about the relationship between local coordinates and crystal
structure is intrinsically missing, thus rendering the the model unable to
reduce further to one result.

Figure 3: The fitted crystal field parameters vs $\chi^{2}$ for Yb2Ti2O7. Each
red point is one converged solution from the PSO with $\chi^{2}<1$. The open
blue circles are the ’true’ solution from Ref. Gaudet  _et al._ (2015). The
open green and black circles are the second and third best solutions found by
CrysFieldExplorer. In panels a), b) d) and f) we observe two separate regions
in the phase space that produce acceptable results. The region where the green
data points reside has not been previously reported in Ref.Gaudet  _et al._
(2015); Bertin  _et al._ (2012). Figure 4: The bulk properties of the
solutions with $\chi^{2}<1$ obtained using PSO. a) The energy vs intensity
spectrum. A fluctuation $\sim$ 5$\%$ of intensity is observed around 76 meV.
b) The inverse susceptibility. The calculated results almost completely
overlap with the ’true’ result. c)The magnetization data. Open blue circle
indicates the true solution, solid red lines indicate the CEF fit. The two
solid grey lines shows an region where all the solutions reside in. d) The
ratio of the anisotropy g-tensors $g_{z}/g_{x(y)}$. The green and black data
in the insets of panels a)-c) represent the solution a) and b) which are also
listed in Table 1.

###  III.2 Analysis on the Crystal Field Parameters of the Low Symmetry Tripod
Kagome Lattice $Er_{3}Mg_{2}Sb_{3}O_{14}$

Now we consider a more challenging case - a tripod Kagome lattice that hosts
various exotic magnetic ground states depending on the constituent rare earth
element Dun  _et al._ (2021). We picked $Er_{3}Mg_{2}Sb_{3}O_{14}$ as the test
case here. The CEF fitting becomes more challenging for a magnetic ion at low
point-symmetry site because the number of none zero crystal field parameters
increases. Fig. 5 shows the local structure of an Er atom in the tripod Kagome
lattice Er3Mg2Sb3O14 with a space group of $R\bar{3}$m. The rare earth ion
Er3+ is surrounded by eight oxygen atoms similar to that in pyrochlore system.
However, the local site symmetry becomes $C_{2h}$ in the tripod Kagome
lattice, much lower than $D_{3d}$ in the pyrochlore lattice. The twofold
rotational $C_{2}$ axis lies in the Kagome plane and is displayed in Fig. 5 as
the green arrow on the center Er atom. As a result, this system requires 15
CEF parameters to fully describe the crystal field Hamiltonian,

$\displaystyle H_{CEF}=$ $\displaystyle
B^{0}_{2}\hat{O}^{0}_{2}+B^{1}_{2}\hat{O}^{1}_{2}+B^{2}_{2}\hat{O}^{2}_{2}+B^{0}_{4}\hat{O}^{0}_{4}+B^{1}_{4}\hat{O}^{1}_{4}+B^{2}_{4}\hat{O}^{2}_{4}$
$\displaystyle+B^{3}_{4}\hat{O}^{3}_{4}+B^{4}_{4}\hat{O}^{4}_{4}+B^{0}_{6}\hat{O}^{0}_{6}+B^{1}_{6}\hat{O}^{1}_{6}+B^{2}_{6}\hat{O}^{2}_{6}+B^{3}_{6}\hat{O}^{3}_{6}$
$\displaystyle+B^{4}_{6}\hat{O}^{4}_{6}+B^{5}_{6}\hat{O}^{5}_{6}+B^{6}_{6}\hat{O}^{6}_{6}.$

According to Hund’s rule, the total, spin and orbit angular momentum for Er3+
are $J$=15/2, $S$=3/2 and $L$=6 respectively. This results in 2$J$+1=16 fold
degeneracy for the ground state of Er3+. Since Er3+ is a Kramer ion, all 16
levels are doublets, 7 excited states are expected to be observed from
inelastic neutron scattering measurement. Dun et al. Dun  _et al._ (2021)
observed 6 out of 7 CEF transitions at 6.4(2), 10.5(3), 21.6(4), 50(1), 65(1)
and 67.5(9) meV. Due to the large amount of none zero CEF parameters, they
were not able to directly fit the CEF transitions with currently existing
software packages such as PyCrystalField. A list of CEF parameters was
obtained from an effective point-charge model which effectively reduced the
number of fitted parameters.

Now we demonstrate how CrysFieldExplorer can help provide insight to this
problem and discuss why directly fitting the CEF parameters with given
experimental data was not feasible. Due to the low site symmetry in
Er3Mg2Sb3O14, the parameter space is magnitudes larger than that in Yb2Ti2O7
and traditional mean-square-root loss function is less sensitive to the 15 CEF
parameters comparing to $L_{\text{Spectrum}}$ as demonstrated previously in
Fig. 1. To handle this case, we need to improve the optimization from our
previous particle swarm optimizer. We continue apply the $L_{\text{Spectrum}}$
and deploy the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) Hansen
(2006), which converges rapidly on complex optimization problems where the
global minimum is extremely sharp in the parameter phase space.

For a local site with low point group symmetry such as monoclinic in the
present case, the choice of principle $xyz$ axis can be somewhat arbitrary. In
the case of tripod Kagome lattice, the $y$-axis is chosen to be parallel to
the two-fold rotational axis $C_{2}$ such that none zero terms of the CEF
parameters are expressed in Eq. 1. The choice of x and z axis is arbitrarily
defined within a perpendicular plane to the $C_{2}$ rotational axis. The x-z
plane is indicated as the green plane in Fig. 5. The dashed and solid arrows
indicate x and z axis are not fixed within the plane perpendicular to y-axis.
Rudowicz et al.Rudowicz (1986, 1985) has demonstrated that within the 15 CEF
parameters, certain rotations exists along $y$-axis that make one of the CEF
parameters 0 . Here we deploy the transformation of rotating the CEF
coordination system w.r.t local y-axis by an angle $\alpha$ such that
$B_{2}^{0}$ is 0 for all the solutions. The rotational angle $\alpha$ is
defined as,

$tan(2\alpha)=\dfrac{B_{2}^{1}}{3\times B_{2}^{0}-B_{2}^{2}}.$ (11)

The rotational transformations of relevant CEF parameters have been tabulated
in Table 2 from Ref.Rudowicz (1985). The CEF parameters reported here have all
been transformed under this standardization rule with each set of parameters
rotated by an angle $\alpha$. A complete list of CEF parameters and the
transformation angle $\alpha$ are plotted in SFig. 1 in the supplemental
material.

We summarize our results in Fig. 6. Different to the Yb2Ti2O7 case,
CrysFieldExplorer can find a large number of solutions across the parameter
space with no obvious segregation. We first look at Fig. 6 c). As Er3Mg2Sb3O14
has 15 (14 independent) CEF parameters, we choose the $B_{2}^{0}$ as an
example and plot the custom defined $\chi^{2}$ loss on a log scale. The solid
grey circles indicate all 150 solutions produced by CrysFieldExplorer in a
time span of approximately 12 hours. We chose 1 as a cutoff number for the
custom defined $\chi^{2}$ and denote the best and worst 10 solutions with
$\chi^{2}<1$ using green and blue open circles. Meanwhile, open blue circles
indicate the true solution. The $\chi^{2}$ for the true solution is
artificially defined as 0 but shifted by a small number of 0.03 in order to be
shown on the log plot.

Figure 6 a) and b) show the neutron data on the energy transfer range from -1
to 33 meV and 25 to 90 meV respectively. Fig. 6 d) is the inverse
susceptibility data. The solid grey lines in these three plots are an overlay
of all the solutions with $\chi^{2}<1$ found by CrysFieldExplorer using CMA-ES
algorithm optimizing $L_{\text{Spectrum}}$. The open blue circle indicates the
solution taken from Ref.Dun  _et al._ (2021) which are treated as true
solution. Fig. 6 a) b) and d) are the three data sets used by
CrysFieldExplorer during the optimization process.

Figure 6 e) compares the calculated magnetizations with the true solution. The
solid red lines and solid red lines mark the top and last 10 solutions with
$\chi^{2}<1$. Following the same color scheme, open blue circles are the true
solution. We drew two solid grey lines to indicate a region where all the
solutions stay within. The determination of the grey lines are guided by eye.

Figure 6 f) shows the ratio of the anistropy g-tensors. The grey squares
indicate the ratio between the largest and smallest while the pink triangle
indicate the largest and 2nd smallest anisotropy g-tensors for all solutions
with $\chi^{2}<1$. The green and red square represent the top and last 10 CEF
fits with $\chi^{2}<1$.

Figure 6 reveals several important features of CrysFieldExplorer when applied
to the Er3Mg2Sb3O14 system. CrysFieldExplorer can perform perfect fits to the
neutron data in panels a) and b) to the degree that is impossible for the
human eyes to distinguish. We performed fitting on both the discrete CEF
energy levels and their relative intensity as well as fitting directly on the
neutron spectrum by applying a resolution limited Lorentzian function to
compute the entire spectrum. There are 73 out of 150 sets of parameters with
$\chi^{2}<1$, all of them produced perfect match to human eyes between the
true solution and CEF fit. Additional constraints were imposed when adding the
inverse susceptibility data to the fit. However, within the $\chi^{2}<1$
range, CrysFieldExplorer was still able to discover multiple solutions that
fit all given experimental observables well.

The powder averaged magnetization data was not considered during the fit
process but used as a benchmark to check the fitting results. Comparing
magnetization data in panel e) of Fig. 6, we start noticing differences
between each solution set. The top 10 solutions show good match to the true
solution at all magnetic fields measured. However, the last 10 CEF fitted
solutions with $\chi^{2}<1$ begin to show small deviations at 4 T magnetic
field. However, we argue our these 10 solutions still fit all the data well
and could still be considered acceptable a few decades ago given the computing
power at that time. This shows the efficiency and accuracy CrysFieldExplorer
exhibits when being deployed on the latest high-performance computers.

Figure 6 f) is particularly interesting because it indicates all of the
solutions produced by CrysFieldExplorer are Ising-like, that is one of the
three $g_{x}$, $g_{y}$ or $g_{z}$ are dominant comparing to the other two. We
do not distinguish $g_{x}$, $g_{y}$ or $g_{z}$ here because there is a certain
degree of freedom in choosing the local frame. From present data, it is not
possible to obtain the relationship between local frame and global frame.
Similar to the Yb2Ti2O7 case, the powder data measured with the mentioned
techniques do not have crystal structure resolution.

Figure 5: The local structure of Er3Mg2Sb3O14 . Green atoms indicate Er,
yellow atoms are Mg and the red atoms are Oxygen. The green arrow on the
center Er atom indicate the $C_{2}$ rotational axis and is chosen as the
y-axis. The green surface is perpendicular to the y-axis. The direction of x
and z axis can be arbitrarily chosen within this plane while preserving the
form of CEF Hamiltonian in Eq. III.2. Figure 6: a), b) Neutron spectroscopy
data plotted within the energy transfer range of [-1, 33] meV and [25,90] meV.
c) Solid grey points are the complete solutions. The red horizontal line marks
$\chi^{2}=1$. Open red circles are the last 10 solutions with $\chi^{2}<1$.
Open green circles are the top 10 solutions. Open blue circle is the true
solution. The y-axis is plotted in log style. d) Inverse susceptibility of the
true solution and solutions with $\chi^{2}<1$. e) Magnetization of the true
solution, top 10 solutions and last 10 solutions with $\chi^{2}<1$. The solid
grey lines are estimates of the upper and lower bounds of the CEF fit with
$\chi^{2}<1$. f) Ratios of the anisotropy g-tensors. Squares indicate
$g_{max}/g_{min}$, triangles indicate $g_{max}/g_{2nd\ min}$. Grey data points
indicate all solutions with $\chi^{2}<1$. Green and red data points show the
top and last 10 CEF fits corresponding to panel c). All the plotted solutions
produce Ising-like anisotropy g-tensors.

To conclude, our results show that with increased free CEF parameters, the
constraints required to narrow down to a single solution also increases
drastically. CrysFieldExplorer was able to produce 73 solutions out of 150
searches that matches neutron data and susceptibility data well. The powder
averaged magnetization data were not used in the fitting because it is
computationally expensive to calculate the powder average. However, single
crystal magnetization data could be added in the fitting without much penalty
on the computational cost to impose more constraints. Despite not being able
narrow down to one or a few solutions, CrysFieldExplorer is able to produce a
list of acceptable solutions. These solutions produce excellent agreement
between the observed and measured data such as inelastic neutron scattering
spectrum, susceptibility, and magnetization measurements. We demonstrate in
the case of Er3Mg2Sb3O14 , neutron spectroscopy and susceptibility do not
impose enough constraints to fully describe the problem, thus making it an
under-determined case when only considering limited experimental observables.
Therefore, we agree with the statement in the original paper that direct
fitting of the CEF parameter is not possible Dun  _et al._ (2021).
Nonetheless, in the Er3Mg2Sb3O14 system, by wrapping $L_{\text{Spectrum}}$
into CMA-ES, we show the strong capability and high efficiency
CrysFieldExplorer exhibits. The approximated time span to generate 150 sets of
solutions is approximately 12 hours. We argue 73 of the 150 solutions produce
acceptable results. Future updates can further improve the efficiency by
leveraging high performance computing on multi-core processors on user-end
computers. CrysFieldExplorer’s results also suggest that when fitting a large
number of CEF parameters, user should exercise extra caution when choosing the
correct solution. In some of the previous research, a solution is believed to
be valid when the neutron scattering data, susceptibility data, magnetization
measurements, anisotropy g-tensors, etc are all consistent. However, utilizing
CrysFieldExplorer’s efficient algorithm, we show that these condition can all
be met with very good agreement between experimental and fitted results in low
site symmetry materials such as Er3Mg2Sb3O14 with 15 (14 independent) none
zero CEF parameters. As the number of CEF parameter increases, more
experimental data should be acquired to fit the CEF Hamiltonian with
experimental model to obtain conclusive results. These can be, but not limited
to, polarized neutron diffraction, specific heat etc. Future researchers
should be cautious to draw conclusions on CEF analysis when considering
limited experimental data.

##  IV Implementation

We aim to provide the most straight forward implementations for
CrysFieldExplorer to the research community. Users first need to specify the
magnetic ions and local site symmetry. This allows CrysFieldExplorer to
generate a list of none-zero CEF parameters. Then the total, orbital and spin
angular momentum $J$, $L$ and $S$ can be automatically determined. The program
will construct the Hamiltonian up to 6th order using the Stevens operators as
well as the matrix representation of $J_{plus}$, $J_{minus}$ and $J_{z}$. Then
users can import the experimentally measured data into the program for
fitting. In the present version of CrysFieldExplorer, it is able to fit
neutron scattering data, susceptibility, single crystal magnetization,
anisotropy g-tensor data all together with different weights. For neutron
scattering data, relative intensities of the observed transitions from the
ground state need to be specified by users from a data reduction software. The
currently required data format is a two-column text file. The experimental
data will be loaded into the optimization program. For high-symmetry cases the
program will select particle swarm optimization and for low-symmetry system a
more efficient CMA-ES will be chosen instead. The loss function will be
predetermined as discussed before. Although CrysFieldExplorer does not require
a set of specific starting parameters, a suitable range for all CEF parameters
still needs to be given. However, this can be estimated from the magnitudes of
the transition energy levels from the inelastic neutron scattering. The
optimization strategy will minimize Spectrum - Characteristic Loss first. Once
Spectrum - Characteristic Loss reaches -10, it will start optimizing magnetic
measurements or any other user defined data. Users are able to prioritize
different data by fine tuning their respective weights. By default,
CrysFieldExplorer will generate 100 converged CEF parameters in .csv file.
These parameters can be checked with either CrysFieldExplorer or a 3rd party
software for consistency. The detailed implementation of CrysFieldExplorer can
be found at https://github.com/KyleQianliMa/CrisFieldExplorer.

##  V Use and Limitations

The program CrysFieldExplorer presented here provides a fast converging,
gradient free method that scans through large CEF parameter phase space to
provide satisfactory results based on user input. It bypasses the need of
determining initial parameters from point-charge model which can be difficult
for inexperienced users. CrysFieldExplorer contains two optimization methods,
PSO and CMA-ES, to minimize a set of specially designed loss functions. The
choice of each optimization method is generally chosen such that for CEF
parameters $<6$, PSO is chosen, otherwise CMA-ES is preferred. Although
CrysFieldExplorer can successfully produce reliable solutions, it cannot
interpret the physical meaning of each solution. Furthermore, if the problem
is under-defined, that is experimental data does not pose enough constraints
to the CEF Hamiltonian, multiple solutions produced by CrysFieldExplorer can
be confusing to determine which solution is the correct solution. Although we
provide two optimization methods, there has not been a standard way to choose
the best optimization methods for different problems. The weight for different
losses can also be tricky to balance. In some test cases it is easier to fit
one experiment than another. Finally, the current version of CrysFieldExplorer
does not provide a user interface. It requires users with basic knowledge of
Python programming language.

##  VI Conclusion

We have developed a lite python-based program CrysFieldExplorer to provide an
efficient procedure of fitting the CEF parameters from experimental data.
CrysFieldExplorer uses two evolutionary algorithms - particle swarm
optimization (PSO) and Covariance Matrix Adaptation Evolution Strategy (CMA-
ES) to minimize a newly designed loss function - Spectrum Characteristic loss
$L_{\text{Spectrum}}$. Comparing $L_{\text{Spectrum}}$ function with
traditional mean-square-root loss, we showed that $L_{\text{Spectrum}}$ is
more sensitive to input parameters, significantly reduces the steepness of
energy barriers on the parameter-loss space and has a better global structure
around local and global minima. This allows CrysFieldExplorer to start the
fitting procedure directly from the CEF Hamiltonian without prior knowledge of
starting parameters. Although a fitting range of the CEF parameter is
required, the range can be generously given such that it matches the
magnitudes of the excited CEF energy levels. To demonstrate the use of
CrysFieldExplorer, we have provided two examples in two different systems to
show that our algorithm can discover multiple solutions that matches
experimental data equally well. Our findings suggest that as the number of to-
be-fitted CEF parameters increases, traditional experimental observables such
as neutron scattering data, susceptibility and magnetization measurement may
not impose enough constraints to fully describe the problem. Furthermore,
neither the powder neutron spectroscopy nor the powder susceptibility data
have crystal structure resolution to define the connection between local ion
and crystal structure frame. This arbitrariness could further produce
confusion when determining the correct CEF fit. Therefore, the present
standard of reporting CEF parameters may not be conclusive as the CEF
parameter phase space is highly degenerate. When conducting CEF analysis with
limited experimental data, further measurements such as polarized neutron
diffraction can provide additional information to help better determine the
correct CEF Hamiltonian.

##  VII Acknowledgments

The research at Oak Ridge National Laboratory (ORNL) was supported by the U.S.
Department of Energy (DOE), Office of Science, Office of Advanced Scientific
Computing Research under the contract ERKJ387, and Office of Basic Energy
Sciences, Early Career Research Program Award KC0402020, under Contract DE-
AC05-00OR22725. A portion of this research used resources at the High Flux
Isotope Reactor and the Spallation Neutron Source, DOE Office of Science User
Facilities operated by the Oak Ridge National Laboratory.

## References

  * Gingras  _et al._ (2000) M. J. P. Gingras, B. C. den Hertog, M. Faucher, J. S. Gardner, S. R. Dunsiger, L. J. Chang, B. D. Gaulin, N. P. Raju, and J. E. Greedan, “Thermodynamic and single-ion properties of ${\mathrm{Tb}}^{3+}$ within the collective paramagnetic-spin liquid state of the frustrated pyrochlore antiferromagnet ${\mathrm{Tb}}_{2}{\mathrm{Ti}}_{2}{\mathrm{O}}_{7}$,” Phys. Rev. B 62, 6496–6511 (2000).
  * Gardner  _et al._ (2010) Jason S. Gardner, Michel J. P. Gingras, and John E. Greedan, “Magnetic pyrochlore oxides,” Rev. Mod. Phys. 82, 53–107 (2010).
  * Cao  _et al._ (2009a) H. Cao, A. Gukasov, I. Mirebeau, P. Bonville, C. Decorse, and G. Dhalenne, “Ising versus $XY$ Anisotropy in Frustrated ${R}_{2}{\mathrm{Ti}}_{2}{\mathbf{O}}_{7}$ Compounds as “Seen” by Polarized Neutrons,” Phys. Rev. Lett. 103, 056402 (2009a).
  * Bertin  _et al._ (2012) A Bertin, Y Chapuis, P Dalmas de Réotier, and A Yaouanc, “Crystal electric field in the R2Ti2O7 pyrochlore compounds,” Journal of Physics: Condensed Matter 24, 256003 (2012).
  * Gaudet  _et al._ (2015) J. Gaudet, D. D. Maharaj, G. Sala, E. Kermarrec, K. A. Ross, H. A. Dabkowska, A. I. Kolesnikov, G. E. Granroth, and B. D. Gaulin, “Neutron spectroscopic study of crystalline electric field excitations in stoichiometric and lightly stuffed ${\mathrm{Yb}}_{2}{\mathrm{Ti}}_{2}{\mathrm{O}}_{7}$,” Phys. Rev. B 92, 134420 (2015).
  * Dun  _et al._ (2021) Zhiling Dun, Xiaojian Bai, Matthew B. Stone, Haidong Zhou, and Martin Mourigal, “Effective point-charge analysis of crystal fields: Application to rare-earth pyrochlores and tripod kagome magnets $\mathrm{R}{}_{3}\mathrm{Mg}{}_{2}\mathrm{Sb}{}_{3}\mathrm{O}{}_{14}$,” Phys. Rev. Research 3, 023012 (2021).
  * Paddison  _et al._ (2017) Joseph AM Paddison, Marcus Daum, Zhiling Dun, Georg Ehlers, Yaohua Liu, Matthew B Stone, Haidong Zhou, and Martin Mourigal, “Continuous excitations of the triangular-lattice quantum spin liquid YbMgGaO4,” Nature Physics 13, 117–122 (2017).
  * Gaudet  _et al._ (2019) J. Gaudet, E. M. Smith, J. Dudemaine, J. Beare, C. R. C. Buhariwalla, N. P. Butch, M. B. Stone, A. I. Kolesnikov, Guangyong Xu, D. R. Yahne, K. A. Ross, C. A. Marjerrison, J. D. Garrett, G. M. Luke, A. D. Bianchi, and B. D. Gaulin, “Quantum Spin Ice Dynamics in the Dipole-Octupole Pyrochlore Magnet ${\mathrm{Ce}}_{2}{\mathrm{Zr}}_{2}{\mathrm{O}}_{7}$,” Phys. Rev. Lett. 122, 187201 (2019).
  * Bramwell and Gingras (2001) Steven T. Bramwell and Michel J. P. Gingras, “Spin ice state in frustrated magnetic pyrochlore materials,” Science 294, 1495–1501 (2001), https://www.science.org/doi/pdf/10.1126/science.1064761 .
  * Ross  _et al._ (2011a) Kate A. Ross, Lucile Savary, Bruce D. Gaulin, and Leon Balents, “Quantum excitations in quantum spin ice,” Phys. Rev. X 1, 021002 (2011a).
  * Porée  _et al._ (2022) Victor Porée, Elsa Lhotel, Sylvain Petit, Aleksandra Krajewska, Pascal Puphal, Adam H. Clark, Vladimir Pomjakushin, Helen C. Walker, Nicolas Gauthier, Dariusz J. Gawryluk, and Romain Sibille, “Crystal-field states and defect levels in candidate quantum spin ice ${\mathrm{Ce}}_{2}{\mathrm{Hf}}_{2}{\mathrm{O}}_{7}$,” Phys. Rev. Mater. 6, 044406 (2022).
  * Keimer and Moore (2017) B Keimer and JE Moore, “The physics of quantum materials,” Nature Physics 13, 1045–1055 (2017).
  * Smith  _et al._ (2022) E. M. Smith, O. Benton, D. R. Yahne, B. Placke, R. Schäfer, J. Gaudet, J. Dudemaine, A. Fitterman, J. Beare, A. R. Wildes, S. Bhattacharya, T. DeLazzer, C. R. C. Buhariwalla, N. P. Butch, R. Movshovich, J. D. Garrett, C. A. Marjerrison, J. P. Clancy, E. Kermarrec, G. M. Luke, A. D. Bianchi, K. A. Ross, and B. D. Gaulin, “Case for a ${\mathrm{U}(1)}_{\pi}$ Quantum Spin Liquid Ground State in the Dipole-Octupole Pyrochlore ${\mathrm{Ce}}_{2}{\mathrm{Zr}}_{2}{\mathrm{O}}_{7}$,” Phys. Rev. X 12, 021015 (2022).
  * Fulde  _et al._ (1970) Peter Fulde, LL Hirst, and Alan Luther, “Superconductors containing impurities with crystal-field split energy levels,” Zeitschrift für Physik A Hadrons and nuclei 230, 155–177 (1970).
  * Mesot and Furrer (1997) Joël Mesot and A Furrer, “The crystal field in rare earth based high-temperature superconductors,” Journal of superconductivity 10, 623–643 (1997).
  * Lynn  _et al._ (1980) J. W. Lynn, D. E. Moncton, L. Passell, and W. Thomlinson, “Magnetic correlations and crystal-field levels in the superconductor (${\mathrm{Ce}}_{0.73}$${\mathrm{Ho}}_{0.27}$) ${\mathrm{Ru}}_{2}$,” Phys. Rev. B 21, 70–78 (1980).
  * Metoki  _et al._ (2004) N Metoki, K Kaneko, TD Matsuda, A Galatanu, T Takeuchi, S Hashimoto, T Ueda, R Settai, Y Ōnuki, and N Bernhoeft, “Magnetic structure and the crystal field excitation in heavy-fermion antiferromagnetic superconductor cept3si,” Journal of Physics: Condensed Matter 16, L207 (2004).
  * Goremychkin  _et al._ (2004) E. A. Goremychkin, R. Osborn, E. D. Bauer, M. B. Maple, N. A. Frederick, W. M. Yuhasz, F. M. Woodward, and J. W. Lynn, “Crystal Field Potential of ${\mathrm{P}\mathrm{r}\mathrm{O}\mathrm{s}}_{4}{\mathrm{S}\mathrm{b}}_{12}$: Consequences for Superconductivity,” Phys. Rev. Lett. 93, 157003 (2004).
  * Rubio Temprano  _et al._ (2000) D. Rubio Temprano, J. Mesot, S. Janssen, K. Conder, A. Furrer, H. Mutka, and K. A. Müller, “Large Isotope Effect on the Pseudogap in the High-Temperature Superconductor ${\mathrm{HoBa}}_{2}{\mathrm{Cu}}_{4}{O}_{8}$,” Phys. Rev. Lett. 84, 1990–1993 (2000).
  * Boothroyd  _et al._ (1992) A. T. Boothroyd, S. M. Doyle, D. MK. Paul, and R. Osborn, “Crystal-field excitations in ${\mathrm{Nd}}_{2}$${\mathrm{CuO}}_{4}$, ${\mathrm{Pr}}_{2}$${\mathrm{CuO}}_{4}$, and related n-type superconductors,” Phys. Rev. B 45, 10075–10086 (1992).
  * Pickett (2021) Warren E Pickett, “The dawn of the nickel age of superconductivity,” Nature Reviews Physics 3, 7–8 (2021).
  * Stevens (1952) K W H Stevens, “Matrix elements and operator equivalents connected with the magnetic properties of rare earth ions,” Proceedings of the Physical Society. Section A 65, 209–215 (1952).
  * Hutchings (1964) M.T. Hutchings, “Point-charge calculations of energy levels of magnetic ions in crystalline electric fields,” (Academic Press, 1964) pp. 227–273.
  * Lea  _et al._ (1962) KR Lea, MJM Leask, and WP Wolf, “The raising of angular momentum degeneracy of f-electron terms by cubic crystal fields,” Journal of Physics and Chemistry of Solids 23, 1381–1405 (1962).
  * Walter (1984) U Walter, “Treating crystal field parameters in lower than cubic symmetries,” journal of physics and chemistry of solids 45, 401–408 (1984).
  * Loewenhaupt (1990) M Loewenhaupt, “Crystal fields in low symmetry systems,” Physica B: Condensed Matter 163, 479–482 (1990).
  * Boothroyd (2014) AT Boothroyd, “Spectre, a program for calculating spectroscopic properties of rare earth ions in crystals,” (2014).
  * Rotter (2013) M. Rotter, “Mcphase,”  (2013).
  * Arnold  _et al._ (2014) O. Arnold, J.C. Bilheux, J.M. Borreguero, A. Buts, S.I. Campbell, L. Chapon, M. Doucet, N. Draper, R. Ferraz Leal, M.A. Gigg, V.E. Lynch, A. Markvardsen, D.J. Mikkelson, R.L. Mikkelson, R. Miller, K. Palmen, P. Parker, G. Passos, T.G. Perring, P.F. Peterson, S. Ren, M.A. Reuter, A.T. Savici, J.W. Taylor, R.J. Taylor, R. Tolchenov, W. Zhou, and J. Zikovsky, “Mantid—Data analysis and visualization package for neutron scattering and $\mu$SR experiments,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 764, 156–166 (2014).
  * Baldovi  _et al._ (2013) Jose J Baldovi, Salvador Cardona-Serra, Juan M Clemente-Juan, Eugenio Coronado, Alejandro Gaita-Ariño, and Andrew Palii, “Simpre: A software package to calculate crystal field parameters, energy levels, and magnetic properties on mononuclear lanthanoid complexes based on charge distributions,” (2013).
  * R.Osborn (2017) R.Osborn, “Cfcal,”  (2017).
  * Fabi (1995) Peter Fabi,  _FOCUS: an interactive crystal electric field parameter fitting package using neutron scattering data_ (Council for the Central Laboratory of the Research Councils, 1995).
  * Scheie (2021) Allen Scheie, “Pycrystalfield: software for calculation, analysis and fitting of crystal electric field hamiltonians,” Journal of Applied Crystallography 54, 356–362 (2021).
  * Bonyadi and Michalewicz (2017) Mohammad Reza Bonyadi and Zbigniew Michalewicz, “Particle Swarm Optimization for Single Objective Continuous Space Problems: A Review,” Evolutionary Computation 25, 1–54 (2017).
  * Hansen (2006) Nikolaus Hansen, “The CMA evolution strategy: a comparing review,” Towards a new evolutionary computation: Advances in the estimation of distribution algorithms , 75–102 (2006).
  * Tokura  _et al._ (2017) Yoshinori Tokura, Masashi Kawasaki, and Naoto Nagaosa, “Emergent functions of quantum materials,” Nature Physics 13, 1056–1068 (2017).
  * Ball (2017) Philip Ball, “Quantum materials: Where many paths meet,” MRS Bulletin 42, 698–705 (2017).
  * Giustino  _et al._ (2021) Feliciano Giustino, Jin Hong Lee, Felix Trier, Manuel Bibes, Stephen M Winter, Roser Valentí, Young-Woo Son, Louis Taillefer, Christoph Heil, Adriana I Figueroa, _et al._ , “The 2021 quantum materials roadmap,” Journal of Physics: Materials 3, 042006 (2021).
  * Cava  _et al._ (2021) Robert Cava, Nathalie De Leon, and Weiwei Xie, “Introduction: quantum materials,” (2021).
  * Bassman  _et al._ (2020) Lindsay Bassman, Sahil Gulania, Connor Powers, Rongpeng Li, Thomas Linker, Kuang Liu, TK Satish Kumar, Rajiv K Kalia, Aiichiro Nakano, and Priya Vashishta, “Domain-specific compilers for dynamic simulations of quantum materials on quantum computers,” Quantum Science and Technology 6, 014007 (2020).
  * Lau  _et al._ (2020) Chun Ning Lau, Fengnian Xia, and Linyou Cao, “Emergent quantum materials,” MRS Bulletin 45, 340–347 (2020).
  * Stanev  _et al._ (2021) Valentin Stanev, Kamal Choudhary, Aaron Gilad Kusne, Johnpierre Paglione, and Ichiro Takeuchi, “Artificial intelligence for search and discovery of quantum materials,” Communications Materials 2, 105 (2021).
  * Kimura  _et al._ (2013) Kenta Kimura, S Nakatsuji, JJ Wen, C Broholm, MB Stone, E Nishibori, and H Sawa, “Quantum fluctuations in spin-ice-like pr2zr2o7,” Nature communications 4, 1934 (2013).
  * Martin  _et al._ (2017) N. Martin, P. Bonville, E. Lhotel, S. Guitteny, A. Wildes, C. Decorse, M. Ciomaga Hatnean, G. Balakrishnan, I. Mirebeau, and S. Petit, “Disorder and quantum spin ice,” Phys. Rev. X 7, 041028 (2017).
  * Gingras and McClarty (2014) Michel JP Gingras and Paul A McClarty, “Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets,” Reports on Progress in Physics 77, 056501 (2014).
  * Ross  _et al._ (2014) K. A. Ross, Y. Qiu, J. R. D. Copley, H. A. Dabkowska, and B. D. Gaulin, “Order by Disorder Spin Wave Gap in the $XY$ Pyrochlore Magnet ${\mathrm{Er}}_{2}{\mathrm{Ti}}_{2}{\mathrm{O}}_{7}$,” Phys. Rev. Lett. 112, 057201 (2014).
  * Bhardwaj  _et al._ (2022) Anish Bhardwaj, Shu Zhang, Han Yan, Roderich Moessner, Andriy H Nevidomskyy, and Hitesh J Changlani, “Sleuthing out exotic quantum spin liquidity in the pyrochlore magnet Ce2Zr2O7,” npj Quantum Materials 7, 51 (2022).
  * Cao  _et al._ (2009b) H B Cao, A Gukasov, I Mirebeau, and P Bonville, “Anisotropic exchange in frustrated pyrochlore Yb2Ti2O7,” Journal of Physics: Condensed Matter 21, 492202 (2009b).
  * Ross  _et al._ (2011b) Kate A. Ross, Lucile Savary, Bruce D. Gaulin, and Leon Balents, “Quantum excitations in quantum spin ice,” Phys. Rev. X 1, 021002 (2011b).
  * Rudowicz (1986) Czesław Rudowicz, “On standardization and algebraic symmetry of the ligand field hamiltonian for rare earth ions at monoclinic symmetry sites,” The Journal of chemical physics 84, 5045–5058 (1986).
  * Rudowicz (1985) Czesław Rudowicz, “Relations between arbitrary symmetry spin-hamiltonian parameters bkq and bkq in various axis systems,” Journal of Magnetic Resonance (1969) 63, 95–106 (1985).