Physics-informed machine learning for inorganic scintillator discovery

Scintillators—a class of materials that emit light under excitation by ionizing radiation—are important for a range of applications including medical imaging, oil exploration, astrophysics, nuclear materials detection for nonproliferation and homeland security, and detection of x-rays, γ-rays, and neutrons in high energy physics experiments.^1–3 A scintillator can absorb radiation energy via creation of electron-hole pairs, which subsequently recombine to emit photons. To facilitate the scintillation process, dopants—also referred to as activators—are frequently used in a host material. Extrinsic lanthanide dopants, such as Ce³⁺, often serve as luminescent centers in fast and bright scintillators,^4–10 owing to their allowed and relatively fast 5d–4f transition. The overall scintillation process can broadly be divided into four distinct processes, namely, energy conversion, thermalization, transfer to luminescent centers, and light emission.^1,2,11 In the energy conversion process, a large part of the ionizing radiation is absorbed by using the scintillator to form “hot” electrons and holes. These charge carriers then thermalize through inelastic processes and subsequently lead to the formation of excitonic states and localize at luminescent centers. The relaxation of excited luminescent centers then leads to the emission of scintillation light through radiative recombination. Therefore, as a strict requirement for a working Ce-activated scintillator, the corresponding 4f (i.e., the ground state with the electronic configuration of [Xe]4fⁿ 5d⁰) and 5d₁ (i.e., the lowest excited state with the electronic configuration of [Xe]4fⁿ⁻¹ 5d¹) levels induced by using the activator must lie within the bandgap of the host material.^12–14 Furthermore, there are a number of additional application-specific desirable properties in terms of stopping power, light yield per unit energy of incident radiation, energy resolution, speed, and cost that a good scintillator material should satisfy. In principle, these requirements can be directly translated to constraints on fundamental material properties such as electronic structure and defect energetics, as well as transport, optical, chemical, and structural properties.^2,15 Furthermore, the valence and conduction band dispersions control the mobility of the thermalized charge carriers. Therefore, these features of the electronic structure are closely related to the light-yield and the time-scale of the scintillation process.^16–20 

While there are a range of materials properties that dictate the performance of a scintillator, the 4f and 5d-shell electron binding energies of doped lanthanides in wide bandgap compounds (relative to the binding energies of electrons at the top of the valence band and the bottom of the conduction band, respectively) are perhaps the most important material characteristics that dictate the suitability of any host compound as a viable scintillator. For instance, scintillation in Ce-doped materials corresponds to a transition from the excited state where the lowest Ce 5d level is filled to the ground state where a single 4f level is filled. The position of the ground state 4f¹ configuration with respect to the valence band maximum (VBM) in a Ce³⁺ doped host determines how quickly and efficiently a hole generated during the energy conversion process will localize at the luminescent center. In Ce-doped oxides, the Ce 4f level usually lies low in the bandgap close to the VBM, as a result of which some of these systems can lead to a very fast and efficient scintillator (for example Lu₂SiO₅:Ce or LSO:Ce). On the other hand, Ce-doped fluorides often exhibit a poor light yield because of the Ce 4f levels that lie about more than 4 eV above the VBM, leading to a very poor hole capture probability at the activator site.²¹ Similarly, too large an energy difference in the conduction band minimum (CBM) and the excited state level would also decrease the probability of Ce³⁺ to capture an electron since more phonon modes would need to be available to dissipate the large energy and therefore reducing the probability of electron capture. Furthermore, if the Ce excited state, with the electronic configuration of [Xe]4f⁰ 5d¹, lies below but too close to the conduction band minimum (CBM), then thermal excitation from the 5d state into the conduction band can, in principle, reduce or even completely quench luminescence, thereby leading to a poor scintillator. Finally, the position of energetically favorable defect levels relative to the 4f and 5d transition levels in the bandgap also dictates the nature of possible trapping mechanisms in the host that can quench or reduce the transfer of energy from the incident radiation to the activator site. These notions are schematically captured in Fig. 1(a) (see, for instance, Ref. 1 for a more detailed discussion on scintillation mechanisms).

Owing to the importance of the 4f–5d transition levels in lanthanide doped solids, much effort has recently been devoted in the literature to understanding the systematics and guiding principles that govern the relative chemical trends exhibited by these levels going from one host chemistry to the other. A significant amount of this understanding has resulted due to the pioneering work of Dorenbos and collaborators,^{13,14,21–27} who have systematically compiled and analyzed the available experimental data on these transition levels in a diverse range of chemistries and crystal structures in a series of papers spanning over a decade. This body of work has also led to physics-based models for predicting properties of the 4f–5d transitions, which can be used to estimate the positioning of these states in the bandgap given limited experimental input.^28,29 However, these models rely on one or more experimentally measured parameters as input to pin-down the lanthanide activator’s 4f and 5d₁ energy levels relative to the band edges of the host compound.

Here, further building on the past work, we develop new machine learning (ML) models that—once trained on available experimental data—can directly predict the key parameters required for making reliable predictions of the Ce-activator’s energy levels for any given host chemistry. In combination with the physics-based empirical models, the developed ML models can reliably predict the chemical trends for the host-dependent electron binding energies within the 4f and 5d₁ energy levels of the lanthanide activator to estimate the relative positions of the activator’s energy levels relative to the valence and conduction band edges of any given host chemistry. As a result, the ML-models enable high-throughput chemical space exploration to search for potential novel inorganic scintillator chemistries, in a much more efficient manner than the current state-of-the art quantum mechanical computations or experimental measurements allow for. After quantifying the prediction performance of the models on unseen data, we demonstrate that the present approach can not only capture systematic and physically meaningful chemical trends across diverse host chemistries but is also equally effective in screening potential candidate compounds in a high-throughput manner while exploring vast chemical spaces. However, before we describe our findings in greater detail, a body of work relevant to the present study has been reviewed in further detail in Sec. II.

The energies of the 4f–5d transitions and the positions of the 4f and 5d transition levels (taken with respect to the VBM and CBM of the host or taken with respect to the vacuum level) of divalent and trivalent lanthanide dopants in compounds depend strongly on the type of lanthanide, its valence state, and the nature of the chemical environment in the compound.^30–33 Despite this large variability, there are also significant systematic trends in these energies, when moving across the lanthanide series as well as when comparing different compounds doped with a common lanthanide ion. In the past decade, these chemical trends have been analyzed in great detail by Dorenbos and others.^22–26 Some of these studies have also resulted in empirical models, which can be used to estimate the relative (i.e., referred to the host’s valence and conduction band edges) or the absolute (i.e., with respect to the vacuum level) binding energies of electrons in the 4f and 5d transition levels.^29,34 In this section, we briefly review some of this past work, of the informatics models developed in the present work.

Going through the lanthanide series in the periodic table of the elements, Fig. 1(b) shows the electronic configurations (in their charge neutral states) and commonly exhibited oxidation states of lanthanide impurity ions in compounds. It can be seen that while the most commonly exhibited oxidation state for most lanthanides is +3, all of them can exhibit both the +2 and +3 oxidation states. In a recent study, experimental data on 4f-electron binding energies for a number of lanthanides in the trivalent [Xe]4fⁿ and the divalent [Xe]4fⁿ⁺¹ configurations in three different chemical environments (namely, the free lanthanide ions, in the pure lanthanide metals, and the lanthanide ions dissolved in aqueous solutions) were collected and analyzed to establish the absolute 4f-electron binding energies in all divalent and trivalent lanthanides [cf. Fig. 2(a)].²⁸ It is seen that going from the vacuum to an aqueous solution to a metallic environment the shape of the 4f-electron binding curve for both divalent and trivalent lanthanides does not change and the effect of chemical environment is captured by almost rigidly shifting (except a slight contraction tilt, to be discussed below) the two curves (with Q = 2+ and Q = 3+) with respect to each other and moving them closer to the vacuum level. As a result of this analysis, Dorenbos put forward a chemical shift model in which the action of the chemical environment on the 4f-electron binding energy was modeled through the energy difference U(n, A) defined as

where E_4f(n + 1, 2+, A) and E_4f(n, 3+, A) are plotted in Fig. 2(a) for the divalent Q = 2+ and trivalent Q = 3+ lanthanide ions in a host “A” and n varies from 0 to 14 for La³⁺ to Lu³⁺, counting the number of valence electrons in the trivalent ions. U(n, A) is also referred to as the 4f–4f Coulomb repulsion energy, as it expresses the Coulomb repulsion experienced by an electron when it is added to the 4f shell of a lanthanide, in going from a trivalent state to a divalent state that already contains n electrons. For a free Eu ion, U(6, vacuum) amounts to 18.05 eV, as indicated in the left panel of Fig. 2(a).

By choosing Eu as the lanthanide of reference and taking U(6, A)—the 4f electron binding energy difference between Eu²⁺ and Eu³⁺—as a parameter (which can be determined experimentally), Dorenbos has developed an empirical model that provides an estimate for the 4f-electron binding energies of all divalent and all trivalent lanthanide ions relative to the vacuum energy in a given environment A.²⁸ Within this empirical model, the 4f-electron binding energies of Eu in a compound A are given by the following equations:

If the 4f-electron binding energy curves shown in Fig. 2(a) shifted rigidly going from one chemical environment to the other, then the above equations would be sufficient to estimate the 4f-electron binding energy of any divalent or trivalent lanthanide ion relative to the vacuum energy in a given environment A. However, in addition to rigidly shifting, the binding energy curves also change their shape [cf. Fig. 2(b)], which has been described as a contraction-tilt effect where the relative binding energy of a lanthanide other than Eu is modified with a correction term that is proportional to the difference ΔR(m) between its ionic radius and that of Eu. Once the contraction-tilt has been accounted for, the vacuum referred 4f-electron binding energy of any divalent or trivalent lanthanide ion can be estimated as

Here energies are taken in eV and ΔR(m) is in the units of pm. Q represents the charge state of the doped lanthanide ion and m is the number of valence electrons [m = n for a trivalent ion and m = n + 1 for a divalent ion, as defined in Fig. 2(a)]. f typically varies between 0.6 and 0.8 for different fits to experimental data and is often taken as 0.7. The term E(Eu^Q, A) is referred to as the chemical shift.

Owing to a highly localized spatial extent of the radial distribution function of the 4f wave function (as compared to that of the 5d wave function) in lanthanide ions, the relative 4f binding energies are almost entirely controlled by the electron-electron interactions within the 4f-shell. Furthermore, these interactions are very characteristic of a given lanthanide free ion, and the overall chemical environment as well as the local coordination environment around the doped lanthanide ion (i.e., the crystal field) has little influence on the relative ordering of these energies. In fact, this is the reason that the above described framework can be used to compute the 4f-electron binding energies of all divalent and all trivalent lanthanide ions relative to the vacuum energy in any given environment A, just with the knowledge of the 4f binding energies of the free lanthanide ions and the environment-dependent U parameter.

The situation for the 5d binding energies is more complex. Owing to the extended spatial extent and the diffused nature of radial distribution function of the 5d wave function, local environment and chemical bonding have a much more pronounced effect on the 5d binding energies of the lanthanide ions as compared to that of the 4f binding energies. In fact, the crystal field interaction with 5d-electrons is about 50 times stronger than with 4f-electrons and therefore 5d binding energies of lanthanide dopants depend much more strongly on the nature of the host compound.²⁹ 

To further illustrate this, Fig. 2(c) depicts the 4f-shell and 5d-shell electron binding energy shifts for a Ce³⁺ ion going from a vacuum (i.e., an isolated free Ce³⁺ ion in a vacuum) to a host chemical environment. It is clearly evident from the figure that going from the vacuum to the host environment, the 5d levels not only shift with respect to the absolute vacuum level, but also split significantly due to the strong crystal field in the host environment. The difference between the lowest 5d energy levels going from a vacuum to the host environment, when the lowest 4f energy levels in the vacuum [i.e., 4f(Ce³⁺, Free)] and in the host [i.e., 4f(Ce³⁺, A)] are aligned, is defined as the crystal field depression or the spectroscopic red shift [denoted as, D(Ce³⁺, A) in Fig. 2(c)]. Similar to the environment dependent parameter U, the red shift D(Ce³⁺, A) is another key parameter, knowledge of which allows us to estimate the 4f–5d transition energies $Efd1(Ce3+,A)$⁠) for the lanthanide ion (in this case Ce³⁺) doped in a host environment A as follows:

where $Efd1$(Ce³⁺, free), $Efd1$(Ce³⁺, A), and D(Ce³⁺, A) are defined in Fig. 2(c) and the 4f–5d transition energy for a free Ce³⁺ ion in a vacuum, i.e., $Efd1$(Ce³⁺, free) is 6.12 eV. Once the binding energy for the Ce³⁺ 4f level and the transition energy $Efd1$(Ce³⁺, A) in a given host are known, the 5d₁ binding energy can simply be estimated by adding the two terms.

The aim of this work is to use available experimental data on the parameters U(6, A)¹² and D(Ce³⁺, A),¹³ in combination with readily available materials descriptors, to develop validated machine learning (ML) models that allow us to accurately estimate these two parameters for any given host chemistry. Once tested and validated, such ML models can be used with physics-based empirical models (discussed above) and density functional theory (DFT) computations to provide quick estimates of the lowest 4f and 5d energy levels of the lanthanide impurities with respect to the valence and conduction band edges of the host, in any given host environment. Such a physics-informed ML scheme can be used in a first line of screening to quickly screen promising novel candidate materials from vast chemical spaces.

The past experimental and theoretical work has revealed that while there are systematic variations in the 4f- and 5d-shell electron binding energies of a Ce³⁺ ion in a given chemical environment, there can also be significant differences in the energies with respect to varying chemistries and crystal structures. The central goal of the adopted machine learning approach is to establish a mapping between easily accessible attributes and the target property [in our case the parameters U(6, A) and D(Ce³⁺, A), taken separately] for a set of compositions that allows us to learn these chemical trends across compound chemistries in order to make rapid estimates of the 4f and 5d electron binding energies for a Ce³⁺ activator in a given host environment. To this end, we start with a set of elemental features, determined largely on the basis of chemical intuition, available past domain knowledge³⁷ and ease of accessibility. This initial set includes compositionally averaged features such as electronegativity,³⁸ first ionization potential,³⁸ electron affinity,³⁹ and atomic polarization⁴⁰ for the elemental species forming the compound. More specifically, a compositional averaged feature is defined as $∑nifi∑ni$⁠, where n_i denotes the number of elemental species of type i and the corresponding elemental feature f_i in the compound formula and the sum runs over i. We note that while first ionization potential and electron affinity naturally provide relevant energy scales for electronic excitations, electronegativity has been used in the past to model relative chemical trends in the valence and conduction band edges across a range of chemistries.^41,42 The rational for including atomic polarization was based on the fact that polarization effects are deemed to be relevant toward the determination of the centroid shift, i.e., the shift in the d band center of a lanthanide ion going from a vacuum to a host environment.⁴³ Furthermore, it has recently been argued that the charge-weighted compositional averaged electronegativity of cations—defined as $∑niziχi∑nizi$⁠, where z_i and χ_i denote, respectively, the nominal charge state and electronegativity of elemental species of type i in the compound formula—can play a role in determining the spectroscopic polarizability which in turn governs the d level positioning of a lanthanide activator in a host compound relative to that in a vacuum.²⁹ Therefore, we explicitly include the charge weighted averaged electronegativities for both the cation and anion species in our feature set. In addition, to account for influences of positively and negatively charged species separately, we also consider composition-averaged features such as electronegativity, first ionization potential, electron affinity, and atomic polarization with the averages restricted only over either cation or anion species of compounds in the training dataset. Finally, empirical radii⁴⁴ and Pettifor’s Mendeleev’s number⁴⁵ were included in the list to account for relative size effects and chemical similarity across different chemical species, respectively.³⁷ 

In addition to the elemental features discussed above, we found it necessary to include the DFT-computed host bandgap as a feature to explicitly account for dielectric screening in a host environment, which directly dictates the strength of the electron-electron Coulombic repulsion. We note that, in principle, either the bandgap or the electronic dielectric constant (as the bandgap is well known to inversely correlate with the electronic dielectric constant) could have been used as a feature for this purpose.⁴⁶ However, we prefer the DFT bandgap as it can be computed at a relatively lower computational cost. Moreover, large datasets of DFT bandgaps computed at a consistent level of theory for a wide range of chemistries are already available through several open-source materials databases.^47–50 In the present work, we resorted to the Materials Project database⁴⁷ for the bandgap feature set enumeration and the automated workflow was handled using the functionality provided within the PyMatGen library.⁵¹ In total, a set of 20 features constituted our initial feature set, further details of which are captured in Table I. This initial feature set was further reduced based on mutual correlations within the feature set as well as correlations of the individual features with each of the properties, as discussed further in Sec. IV.

The above feature set was used as a starting point to train ML models (discussed below in detail) to predict the two targeted properties, namely, U(6, A) and D(Ce³⁺, A). To establish this mapping, two separate datasets that contain experimentally reported values of U(6, A)¹² and D(Ce³⁺, A)¹³ for a set of compounds for which ground state electronic structures and Kohn-Sham bandgaps (computed using the conventional semi-local density functional theory)⁵² are reported in the Materials Project database. The U(6, A) and the D(Ce³⁺, A) property datasets—containing 109 and 174 compounds, respectively—along with the data collected from the Materials Project database, are provided in the supplementary material accompanying the manuscript.

To establish a general and predictive mapping between a subset of the feature set and the targeted properties, here we employ a well established kernel ridge regression (KRR) model^53,54 with a gaussian kernel, as implemented in the Scikit-learn machine learning suite.⁵⁵ The KRR ML algorithm works on the principles of similarity, quantified using a carefully pre-selected set of relevant features. Subsequently, relying on a basic assumption that similar compounds should exhibit similar properties, an interpolative predictive mapping is established while adhering to the best practices of statistical learning. We note that KRR has recently been very successful for a diverse set of materials problems^56,57 and some examples of application of this approach to materials problems include, but are not limited to, automated property predictions of molecular^58–63 and periodic systems,^64–68 development of adaptive force fields,^69,70 radiation resistance and dielectric breakdown strength prediction,^71–73 self-consistent solutions for quantum mechanics,^74–77 and predictions of bandgaps in solids.⁷⁸ 

Within KRR, the ML estimate of a target property ϕ—in our case the U(6, A) or D(Ce³⁺, A)—of a new system j, is estimated by a sum of weighted kernel functions (i.e., Gaussians) over the entire training set, as

where i runs over the systems in the training dataset, and $|Fij|2=||Fi−Fj||22$⁠, the squared Euclidean distance between the feature vectors F_i and F_j. The coefficients α_is are obtained from the training (or learning) process built on minimizing the following regularized expression:

with $ϕiML$ being the ML estimate for the corresponding experimental value $ϕiExp$⁠. The explicit solution to this minimization problem is α = (K + λI)⁻¹ϕ^Exp, where I is the identity matrix and $Kij=exp−12σ2|Fij|2$ are the kernel matrix elements corresponding to pairs of all materials in the training set. The parameters λ and σ are determined in a five-fold cross validation loop using a logarithmically scaled fine grid. Further technical details on the KRR methodology can be found in Refs. 53 and 54.

We test the prediction performance of the developed ML models by making predictions for the Ce³⁺ dopant’s 4f and 5d levels in A²⁺B⁴⁺O₃ perovskite oxides and elpasolite compounds with a generic formula of A₂(B,B′)X₆. Since the developed ML models use DFT computed bandgaps—consistently computed using the Perdew, Burke, and Ernzerhof (PBE)⁵² generalized gradient approximation (GGA) exchange-correlation functional—as an input feature, the bandgap computations for these datasets were performed using the Vienna Ab Initio Simulation Package (VASP).^79,80

For the oxides, we consider six single perovskite chemistries combinatorially formed by taking A = Ca, Sr or Ba and B = Hf or Zr. Furthermore, to better capture the anticipated chemical trends in the predicted energy levels of the Ce³⁺ activator as a function of host chemistry, in addition to the six ABO₃ perovskites, we also include chemistries formed by a 50% admixture of rocksalt ordered A- and/or B-site cation pairs appearing as nearest neighbors in different columns of the periodic table. This leads to a set of 15 perovskite host chemistries employed here to test the predictions of our ML models. The electronic wave functions were expanded in plane waves up to a cutoff energy of 500 eV. The pseudopotentials based on the projector augmented wave method⁸¹ explicitly included following valence electronic configurations for different elemental species; Ca: 3s²3p⁶4s², Sr: 4s²4p⁶5s², Ba: 5s²5p⁶6s², Zr: 4p⁶4d²5s², Hf: 5p⁶5d²6s² and O: 2s²2p⁴. A gamma-centered automatically generated 4 × 4 × 3 Monkhorst-Pack k-point mesh⁸² was used for Brillouin-zone integrations for a $2×2×2$ cubic (or pseudo-cubic, in the case of the orthorhombic Pnma space group) supercell containing 20 atoms. For smaller or larger supercells, the k-point meshes were appropriately scaled to give the same k-point density in the reciprocal space. To obtain a geometry optimized equilibrium structure, atomic positions and the supercell lattice parameters were fully relaxed using the conjugate gradient method until all the Hellmann-Feynman forces and the stress component were less than 0.01 V/Å and 1.0 × 10⁻² GPa, respectively.

The A₂(B,B′)X₆ elpasolite dataset came from a previous study and features halide compounds (X = F, Cl, Br, and I) with the A- and B-sites occupied by alkali metals Li, Na, K, Rb, and Cs and the B′-site occupied by Sc, Y, La, Ce, Nd, Eu, Gd, and Er.⁸³ In this work, we employ a subset of 200 elpasolites for which bandgaps computed using both a conventional semi-local exchange correlation functional and a hybrid functional are available. The technical details of the DFT computations for the elpasolite dataset can be found in Ref. 83.

Any ML model targeted toward learning a specific material property relies on two main ingredients: the learning algorithm itself and a compact and easily accessible numerical representation (in the form of a multi-dimensional feature vector, also frequently referred to as a descriptor) of the materials in the training datasets. While features listed in Table I—selected largely based on available domain knowledge and chemical intuition—certainly constitute a good starting point, it is not entirely unexpected that some of these features are highly correlated with others and therefore the initial dataset may carry some redundant information. Furthermore, the features that exhibit strong correlations with the two targeted properties are naturally more relevant than the others that show only a poor correlation or almost no correlation with the targets. Therefore, to systematically down-select a compact set of most relevant features, we start by carrying out a Pearson correlation analysis for the 20 features within our initial feature set. The Pearson correlation coefficient $P$ is defined as

for two N dimensional column vectors x and y with $x¯$ and $ȳ$ representing the corresponding means. Here scalars x and y represent a pair of features and the subscript i labels a particular compound in a dataset of N total compounds. It is easy to see that the correlation coefficient $P$ always lies between −1 and +1—capturing the degree of negative and positive correlation, respectively.

Figure 3(a) presents results of the Pearson correlation analysis in the form of a correlogram, i.e., a correlation map of pairwise correlation coefficients computed for each unique pair of features appearing in the initial feature set. The lower and upper triangles show the pairwise correlations computed over the U(6, A) and the D(Ce³⁺, A) property datasets, respectively. The different features are identified using the corresponding feature IDs, as defined in Table I. The filled fraction of the circle in each of the pie charts corresponds to the absolute value of the associated correlation coefficient $P$⁠. In addition, for an easier visualization, the blue and red colors are used to indicate positive and negative correlations, respectively; and the lighter the tone used, the less significant the corresponding correlation. From the correlogram, it can be easily seen that a number of feature pairs are positively correlated, while many others exhibit strong negative correlations. As a next step, we use the computed correlations to further down-select a reduced set of promising features, separately for each targeted property, as follows: If any two features are correlated (or anti-correlated) more than a cutoff value of $Pc1$ (i.e., $|P[Fi,Fj]|>Pc1$⁠) then we include only one of the two features in the down-selected set—the one that shows a higher correlation with the targeted property. Furthermore, we impose a condition that a feature must exhibit a correlation of at least $Pc2$ [i.e., $|P[Fi,ϕ]|⩾Pc2$⁠, with ϕ = U(6, A) or D(Ce³⁺, A)]. Here, Fⁱ and F^j represent two N component column vectors with the superscripts i and j indexing the respective feature IDs and N is the total number of compounds in a target property dataset. Note that by gradually decreasing $Pc1$ and/or increasing $Pc2$ one can systematically reduce the size of the initial feature set. Based on an initial trial-and-error analysis, we find that $Pc1$ = 0.85 and $Pc2$ = 0.25 can be regarded as reasonable cutoff choices which allow us to significantly reduce the size of the initial feature set for computational efficiency without sacrificing the predictive performance of the final ML models developed here.

For the specific choices of $Pc1$ = 0.85 and $Pc2$ = 0.25, only 5 features survive out of the original 20 features for each of the two property datasets. Figures 3(b) and 3(c) present the pairwise correlations among the down-selected feature sets for the U(6, A) and the D(Ce³⁺, A) property datasets, respectively. The top bar in each of the two panels also depicts the correlations of the selected features with the target property. The features—identified by their feature IDs as defined in Table I—are arranged (going from left to right) according to their increasing absolute correlation coefficients with the respective target properties. To get a closer look at the feature-property correlations for the features that exhibit a high degree of correlations with the two targeted properties, Figs. 3(d)–3(i) present scatter plots for the top three features plotted separately against the target property for the U(6, A) and the D(Ce³⁺, A) property datasets, respectively. It can be seen from the figures that some of the features show rather remarkable trends with the properties. Interestingly, while the bandgap (E_g) shows up as an important feature toward prediction of both target properties, the average ionization potential (⁠$IPav$⁠) and the average electronegativity (⁠$χav$⁠) standout for their strong linear correlations with the U(6, A) and the D(Ce³⁺, A) values, respectively. Despite these correlations, Fig. 3 also reveals that there is not one single feature that provides a quantitative capability of predicting any of the target properties. As a natural next step, our aim is to harness these newly identified feature-property trends via employing notions of chemical similarity in a ML model to establish a validated and predictive model that allows for rapid predictions of the target properties for chemistries not included in the training datasets.

We use a KRR-based ML model to quantify the correlations between the down-selected features described above and the two target properties. Details of our ML models are discussed in Sec. III B. For each target property, we train a separate ML model starting from the respective feature sets and training data sets. Since we are interested in the simplest possible ML models (for both better generalizability and computational convenience), our goal is to identify smallest possible subsets of features (out of the already down-selected feature sets) that exhibit highest prediction performance, quantified by their ability to accurately predict the respective target properties of the test set compounds (i.e., evaluated on an unseen dataset). Moreover, in addition to accurate estimates of the target properties, it is also desired to quantify the model’s confidence in those predictions (i.e., the error bars associated with the predictions) with respect to variability in the training data. In other words, as an indicator of the robustness of our ML model, we seek a quantitative measure of the extent to which the ML model predictions change if different subsets of the overall training data are used to train the model. Toward this aim, we carryout a comprehensive analysis by building KRR ML models using all possible subsets of the relevant features (i.e., all unique combinations of Ω features with Ω ∈ [1, 5]), chosen from within the respective down-selected feature sets identified above for each of the targeted property. The root mean square errors (RMSEs) for the U(6, A) and D(Ce³⁺, A) predictions on training and test sets for various models are presented in Figs. 4(a) and 4(b), respectively. In order to account for model prediction variability associated with randomly selected training/test splits, Fig. 4 reports the RMSEs averaged over 100 different randomly selected training/test splits for each of the models.

Figure 4 presents the results of the ML model when different dimensionalities (i.e., Ω) of feature sets are considered for both target properties. We first consider the results of using only one feature to train the model. Also note that the results of ML models built on one-dimensional (1D) feature vectors in Figs. 4(a) and 4(b) are presented in an order which is consistent with the order in which these features appear in Figs. 3(b) and 3(c), respectively. In other words, the 1D features are considered in an increasing order of their correlation with the target properties. As a result, it can be clearly seen that the test set prediction performance of the ML models largely follows the trend exhibited by the feature-property correlation coefficients presented in Fig. 3. For the U(6, A) dataset, while the performance of a single feature, namely, $IPav$ (feature ID: 8), is by far better than all others, the predictive performance of E_g (feature ID: 1) is also clearly distinguished. For the D(Ce³⁺, A) dataset, on the other hand, the top two features, i.e., E_g and $χav$ (feature ID: 3), show a rather comparable (though superior than the other three features) predictive power toward the targeted property.

The 2D models that lead to the lowest RMSEs on the test set data have been marked with a “⋆” in Fig. 4(a) [for the U(6, A) dataset] and Fig. 4(b) [for the D(Ce³⁺, A) dataset]. It is not surprising that, in each case, the best performing 2D feature vector is formed by the two top 1D features. It is also interesting to note that, for both cases, going beyond the best performing 2D models does not lead to a statistically significant improvement in the model prediction performance. For instance, while the best binary feature pair (⁠$IPav$⁠, E_g) for the U(6, A) dataset leads to a test set RMSE of 141 meV, the ML models built on the best 4D and 5D (taking all 5 features considered) feature vectors do not result in any improvements, also leading to RMSEs of 140 meV. A similar behavior can also be seen for the D(Ce³⁺, A) dataset. Since, as a general rule, higher model complexity often leads to poor generalizability, in the case of a comparable prediction performance, a simpler model (i.e., built on a lower dimensional feature set) should always be preferred over a more complex one. Therefore, henceforth we focus our attention on the best performing 2D models.

A comparison between the predictive performance of the 2D ML models for the two target properties, i.e., U(6, A) and D(Ce³⁺, A), reveals considerably different test set RMSEs of 0.14 eV and 0.40 eV, respectively. This indicates that while the environment-dependent U parameter can be predicted to a quite high accuracy, errors in the ML predictions for the spectroscopic red shift parameters are relatively high. This is not very surprising given that (i) the 5d binding energies are known to exhibit a significant local environment and chemical bonding dependence (as compared to the 4f binding energies), owing to the extended spatial extent and the diffused nature of radial distribution function of the orbitals and (ii) the features used to train the D(Ce³⁺, A) ML model are completely agnostic toward the local chemical environment around a Ce activator site. While it is nontrivial and computationally prohibitive to explicitly account for an accurate description of the local chemical environment around a Ce dopant for a large set of diverse chemistries that are considered here—as one needs to perform DFT computations with very large supercells to account for dilute Ce concentrations. Furthermore, since the primary aim of this work is to develop a ML model that enables rapid screening of promising chemistries (based on the predictions of Ce 4f and 5d₁ levels) at a minimal computational effort, in an effort to further boost the predictive performance of the spectroscopic red shift parameters here we pursue an alternative route that allows us to implicitly account for the relevant information.

The classification scheme is depicted in Table II and explicitly accounts for the chemical nature and charge states of cations, anions, and other chemical groups forming the host chemistries. As a next step, we augment our best performing 2D feature pair (E_g, $χav$⁠) with each of the six classification-based features (i.e., features f₁ through f₆) presented in Table II separately to train different ML models. The training and test set prediction performances of these models are captured in Fig. 5. Our analysis indicates that only features that lead to any statistically significant improvements, when included with E_g and $χav$⁠, are f₃ and f₄. Interestingly, when both of these features are used together to augment the best performing 2D feature pair, the resulting feature leads to even lower prediction errors. More specifically, the training and test set RMSEs with the (E_g, $χav$⁠, f₃, f₄) features are 0.31 eV and 0.36 eV, respectively as compared to the 0.38 eV and 0.40 eV, respectively, obtained for the 2D feature pair (E_g, $χav$⁠). Therefore, inclusion of the two best performing classification-based features in the ML model leads to a clear 10% improvement on the test set prediction performance. Further augmentation of this feature set by including additional features did not lead to any considerable improvements in the prediction performance. This can also be seen from Fig. 5, where two ML models trained on (E_g, $χav$⁠, f₃, f₄) and (E_g, $χav$⁠, f₁:f₆) feature sets are shown to exhibit comparable performance. Physical significance of f₃ and f₄ can be realized by noting that these classification-based features are intended to explicitly account for complex groups in ionic crystals can share one or more anions with each other or trivalent metal ions present, respectively. While the complex groups are well known to have a significant influence on the value for the crystal field depression,¹³ the trivalent metal ions often serve as a substitutional site for an trivalent activator due to favorable electrostatic interactions with the local environment provided at that site. Finally we note that, unlike for the cases where one targets to include the local coordination information explicitly in the ML model (through DFT computations or by other means), the identified additional features f₃, f₄ are readily available and do not require any additional computational effort.

While results presented in Figs. 4 and 5 capture the average performance and variability (taken over 100 different runs) for our best performing models (marked with a ⋆) for U(6, A) and D(Ce³⁺, A), respectively, in Figs. 6(a) and 6(b) we present representative parity plots directly comparing the experimental values with the ML predictions for the two properties. In each case, 90% of the dataset was used for training (plotted as squares) with the remaining for testing the model performance (plotted as circles). In each case, we used five different ML runs randomly selecting training and test set splits (depicted by different colors). It can be seen from the parity plots that the prediction performance is comparable for the training and test sets, indicating that there is no overfitting (a problem when a ML model performs very well on a training set but exhibits poor performance on a test set). Furthermore, despite their simplicity (given that we have selected the simplest possible ML models with a minimal set of features in each case), the ML models exhibit good predictive power and stability (predictions do not change drastically over different training/test splits). This highlights the robustness of the model.

While both (i) identification of easily accessible features that correlate with U(6, A) and/or D(Ce³⁺, A) and (ii) the ML-based error analysis quantifying the predictive power of these features toward the target properties are useful and novel, the true power of such a ML-based accelerated property prediction approach is its capability of instantly predicting the targets for any arbitrary host chemistry and crystal structure, without needing to pursue a cumbersome approach of computing defect geometries and the electronic structure from first principles, for instance, via resorting to DFT or many-body electronic structure theories. Note that except for the GGA bandgaps, all other features required for the ML models presented here are readily available and require no additional computational cost. Moreover, the GGA bandgaps can be either accessed through a number of open-source materials property databases or computed at a rather moderate computational cost using conventional DFT simulations. Finally, we note that while the present ML framework provides the 4f and 5d₁ levels for Ce³⁺ doping (with respect to the vacuum level) in a given host chemistry, the host VBM and CBM levels still need to be computed at a suitable level of theory such as hybrid functionals or the GW approach—owing to the well-known tendency of bandgap underestimation by local and semi-local exchange correlation functionals within conventional DFT. However, note that since this step only requires a bulk calculation with a primitive cell for each host chemistry, this can also be performed rather efficiently as compared to the case where one needs to explicitly model a dilute Ce doping concentration by considering a much larger supercell.

After quantifying the accuracy of the developed ML prediction models, next we test the models by making predictions for the Ce³⁺ dopant’s 4f and 5d₁ levels in a range of host chemistries which are practically relevant for scintillator applications. In particular, we test the models by making predictions on a limited set of A²⁺B⁴⁺O₃ perovskite compounds and a large dataset of 200 elpasolite compounds with a generic formula of A₂(B,B′)X₆.

For the oxide dataset, we consider six single perovskite and nine double perovskite chemistries combinatorially formed by taking A = Ca, Sr, or Ba and B = Hf or Zr, as described in Sec. III C. Note that from a scintillation point of view, the perovskite host compounds studied here consistently show high stopping power (owing to the high nuclear charge) and high melting temperatures. Furthermore within this set of compounds, owing to the systematic chemical trends across chemistries, the 4f and 5d levels of Ce³⁺ dopants are expected to vary systematically, some of which have also been reported experimentally, providing a strong motivation to include this dataset as a test case in the present study. We also note that from a structural point of view, the compounds comprised of 100% Ba at the A-site and (Ba_0.5,Sr_0.5)HfO₃, are cubic and crystallize in the $Pm3¯m$ space group (No. 221),⁸⁴ while the others display an orthorhombic crystal structure with the Pnma space group (No. 62),⁸⁵ owing to an energetically favored pattern of BO₆ octahedral rotations (i.e., the a⁻b⁺a⁻ tilt pattern configuration according to Glazer’s notation).⁸⁶ In the present study, each of the compounds was therefore modeled using its appropriate DFT ground state geometry and electronic structure.

While the predictions for the 4f and 5d₁ levels can be made by employing the trained ML models to predict U(6, A) and/or D(Ce³⁺, A) first and subsequently using Eqs. (1)–(5), to locate their relative positions with respect to band edges, one also needs to compute the bandgaps as well as the VBM and CBM positions with respect to the vacuum level for each of the host chemistries. Although relative chemical trends in the bandgaps are expected to be correct within conventional DFT, the actual bandgap values are grossly underestimated as compared to the experimentally measured ones. In fact, this is a well-known tendency for the local and semi-local exchange correlation functionals within DFT (including the PBE functional employed here) which has been attributed to the inherent lack of derivative discontinuity⁸⁷ and delocalization error^88,89 within the local or semi-local exchange-correlation functionals.⁹⁰ To obtain better predictions of the bandgaps, we resorted to a more accurate (and also much more computationally expensive) approach utilizing hybrid Heyd-Scuseria-Ernzerhof (HSE) functional calculations^91,92 on top of the relaxed PBE geometry and PBE eigenfunctions.⁵² Within this family of functionals, one has access to the parameters α and ω, representing the fraction of exact exchange and the inverse screening length, respectively, that can be used to tune the details of the electronic structure within a given class of materials. The standard HSE06 functional is obtained with α = 0.25 and ω = 0.207 Å⁻¹. A comparison of the bandgaps computed using the HSE06 functional with the corresponding experimental values reveals that the HSE06 bandgaps—though representing a significant improvement as compared to the PBE bandgaps—are still considerably underestimated. To further improve the agreement between the computed and experimental bandgaps, we systematically vary the fraction of the Hartree-Fock exact exchange α in the DFT hybrid functional (while fixing the inverse screening length ω to that of the standard HSE06 value) to track the corresponding variation in the bandgap. Our analysis reveals that the bulk bandgaps computed with the HSE functional with α = 0.4 result in excellent agreement with the corresponding experimental bandgap values (cf. Table III in the supplementary material).

Once accurate bandgap values are computed, to reliably determine the relative position of band edges in compounds, a reference state that does not change with chemical composition first needs to be identified. While there are several schemes available in the literature,^5,93–96 we use the binding energy of oxygen 1s core electrons to calculate relative positions of the valence and conduction band edges.^97,98 An experimentally reported value of the valence band edge for BaZrO₃ with respect to the vacuum level from Ref. 27 then allows us to pin the zero of energy.

The stacked band scheme of Fig. 7(a) reveals clear trends in the host band edges and the Ce³⁺ activator’s 4f and 5d₁ levels when the composition changes. The results from the HSE06 and the HSE with α = 0.4 functionals show very similar trends across the A- and B-site chemistries and an increased amount of exact exchange in the latter function lead to a further broadening of the bandgaps largely by pushing the conduction band edges further up, toward the vacuum level. Note that there are systematic variations in the relative band edge alignments, which closely follow chemical trends going down group 2 of the periodic table. As we go from Ca to Ba, both the atomic size and the metallic nature increase and the bandgap decreases. Both the lowering of the CBM and the rise in the VBM contribute to the bandgap reduction, although the changes in the CBM position are greater than those in the VBM position. This observation is in line with the fact that while the valence bands are largely constituted by O-2p orbitals, the A-atoms’ 2s and B-atoms’ d orbitals largely contribute to the conduction bands. Unlike the A-site atomic species, Zr and Hf have a very similar atomic radii and chemical nature, as a result only relatively small changes in the band edge alignment with respect to the B-site atomic species are evident. While the trends in the band edges are quite clear, the dopant levels also exhibit systematic trends with respect to the A- and B-site chemistries, albeit at a relatively small energy scale, which become clear in Figs. 7(b) and 7(c) presenting magnified views of the ML-predicted 5d₁ and 4f levels, respectively.

These trends, coupled with the changes in the conduction band edge, have important consequences for the luminescence properties of the oxides. The stacked band scheme implies that the lowest 5d-level of Ce³⁺ in zirconate compounds will be located above the conduction band edge, eradicating any possibility of a 5d–4f emission, which agrees with observations.²⁷ Unlike the zirconates, Ce³⁺ 5d–4f emission is observed in the three AHfO₃ compounds. Note that while for CaHfO₃ the 5d-level is already below the conduction band edge, for SrHfO₃ it is predicted to be just near the CBM, which is allowed to move further below due to lattice relaxation leading to a Stokes shifted emission. In the case of BaHfO₃, the lowest 5d-state is predicted to be buried in the conduction band and as a result the emission should be fully quenched due to auto-ionization of the 5d-electron to the conduction band. This is also in agreement with a previous study,²⁷ but some studies have also reported 5d–4f emission for BaHfO₃ with radioluminescence efficiency strongly dependent on preparation atmosphere and temperature,⁹⁹ the origins of which remain to be understood. Finally, we note that the vacuum referred binding energy (VRBE) of an electron in the lowest energy 5d-level of Ce³⁺ is predicted to be typically −2 eV, consistently for the entire set of chemistries, again in agreement with experimental observations.²⁷ The data plotted in Fig. 7 are provided as a part of the supplementary material for the paper.

As a second example, to demonstrate the efficacy and efficiency of the developed surrogate models, we use a materials dataset of 200 double perovskite halides (or elpasolites) of A₂BB′X₆-type—a class of materials that is highly relevant for scintillation applications.^100,101 This dataset was recently employed in a study that focused on ML of HSE06 (i.e., using α = 0.25) bandgaps using a multi-fidelity learning model and features halide compounds (X = F, Cl, Br, and I) with the A- and B-sites occupied by alkali metals Li, Na, K, Rb, and Cs and the B′-site occupied by Sc, Y, La, Ce, Nd, Eu, Gd, and Er.⁸³ Similar to the previous example of perovskite oxides, we use the halide ion’s deep core 1s state for the relative valence and conduction band edge alignments (separately within each class of halides) and use experimentally reported vacuum levels for LaX₃ (with X = F, Cl, Br, and I) to set the absolute zero of the energy scale across the halide chemistries.²¹ 

Figure 8 presents our results for the elpasolites dataset, grouped according to the four halide chemistries. Within each of the halide groups, compounds are arranged in an increasing order of their conduction band edge position with respect to the vacuum level. The stacked band levels scheme not only provides an overview on how the valence band and conduction band energies change across the halides, but also depicts the Ce³⁺ 4f ground and first 5d excited states within the bandgap or the conduction band. Also note that, owing to the use of bandgaps computed using HSE06 for the elpasolites in Fig. 8, the predictions on scintillator versus non-scintillator compounds made on the basis of relative placements of the Ce³⁺ 4f and 5d₁ levels with respect to the band edges (i.e., the VBM and CBM, respectively) are going to be rather conservative and false negatives can be expected.

Several interesting observations can be made from Fig. 8. First, while the binding energies in the Ce³⁺ 4f-levels appear largely to be constant for compounds with a given halide chemistry, the first 5d excited state energies exhibit significant variations with respect to the nature of the cation species occupying the A, B, and B′ sites. These intuitive predictions are in line with the previously discussed spatially localized and extended nature of the 4f and 5d wave functions, respectively. These physically meaningful trends are naturally learned by the ML models while training on the two experimental datasets. Second, note that the 4f-VBM energy gap is consistently large for fluorides. The position of the ground state 4f energy levels with respect to the VBM determines how quickly and efficiently a hole generated during the energy conversion process will localize at the luminescent center. For the Ce-doped elpasolite fluorides, the 4f levels are predicted to lie around 5 eV or more above the VBM, indicating a lower hole capture probability for these materials leading to poor scintillation performance. We also note that the large 4f-VBM gap predicted here appears to be a general trend across different fluoride chemistries beyond the elpasolite crystal structure (for instance, cf. Fig. 8 in Ref. 21). Finally, from a scintillation performance point of view, the chlorides and bromides are perhaps the most attractive compounds, a number of which are predicted to have a favorable placement of the Ce³⁺ 4f ground and first 5d excited state levels within the bandgap. If fact a survey of the literature reveals that twenty chloride or bromide compounds in the present dataset have been previously reported as Ce³⁺-activated scintillators, viz. Cs₂(Na,Ce)Cl₆,¹⁰⁵ Cs₂(Na,Gd)Cl₆,¹⁰⁶ Cs₂(Li,Gd)Cl₆,¹⁰² Cs₂(Li,La)Cl₆,¹⁰⁰ Cs₂(Li,Y)Cl₆,¹⁰³ Cs₂(Na,Y)Cl₆,¹⁰⁷ Rb₂(Li,Y)Cl₆,¹⁰⁸ Rb₂(Li,Ce)Cl₆,¹⁰⁹ Cs₂(Na,La)Cl₆,¹¹⁰ Cs₂(Na,La)Br₆,¹¹⁰ Rb₂(Li,Gd)Cl₆,¹¹¹ Cs₂(Li,Ce)Br₆,¹¹² Cs₂(Na,Ce)Br₆,¹¹³ Rb₂(Li,Gd)Br₆,¹¹⁴ Rb₂(Li,Ce)Br₆,¹¹⁵ Cs₂(Na,Gd)Br₆,¹⁰⁶ Cs₂(Li,Y)Br₆,¹¹⁶ Rb₂(Li,Y)Br₆,¹¹⁷ Cs₂(Na,Y)Br₆,¹¹⁸ and Cs₂(Li,Gd)Br₆.¹⁰⁴ All of these are also predicted by the current model to be scintillating with the Ce 4f and 5d states within the bandgap, situated suitably close to the VBM and CBM, respectively. These previously known scintillators are highlighted in blue in Fig. 8. The data presented in Fig. 8 can be found in Table IV of the supplementary material accompanying the manuscript.

Perhaps the most important and practical reason for seeking to develop ML-based models for high-throughput predictions of the 4f ground state and first 5d excited state of lanthanide dopants for a given host chemistry is that the experimental and/or first principles computational approaches to predicting/estimating these levels are exceedingly computation-time and resource intensive. The high cost associated with the conventional approaches severely limits our ability to efficiently screen large chemical space in search of potentially new scintillators. Furthermore, while conventional DFT has well-known limitations in accurately describing f-electron physics (and one often needs to resort to more computationally demanding approaches, such as, hybrid functions and beyond-DFT approaches such as the GW method to alleviate this deficiency), for widespread studies of systems involving larger unit cells containing dilute concentration of lanthanide dopants, the quantum mechanical computational route becomes increasingly impractical. Here we have shown that ML-based efficient property prediction models can be particularly promising in this direction and can provide an alternately pathway for rapid candidate materials screening while exploring vast chemical spaces. Indeed, data-enabled attempts have been undertaken in the past targeted toward predicting various properties for diverse scintillator chemistries but have met with marginal success due either to insufficient incorporation of relevant physics into the developed models, or to the use of features that are not easily accessible for new compounds (for example, features that depend on direct experimental measurements).^37,119,120 The work presented here constitutes an important first step in demonstrating that carefully tested and validated ML-models built upon physically relevant yet readily available features and trained directly on the available past experimental data can in fact lead to efficient yet accurate predictions on physically relevant properties that dictate performance in inorganic scintillators. We also note that while here we have considered the locations of the 4f and 5d₁ levels relative to the band edges as key criteria for scintillation, there are also a number of application-specific performance requirements in terms of, for instance, radiation hardness, stopping power, decay time, light output, and emission wavelength that play a crucial role in scintillator design. While all relevant performance metrics must be addressed in successive development stages, the use of ML models, such as those developed here, can be particularly beneficial in the initial round of screening to make go/no-go decisions. Finally, we also note that the prediction performance of the presented ML models can be systematically improved by incorporating both additional training data as it becomes available and additional relevant features, identified via employing recently developed tools^121,122 to explore combinatorially generated vast non-linear feature spaces. Both of these exciting directions are the subject of future work.

See supplementary material for the U(6, A) and D(Ce³⁺, A) property datasets used here to train the developed machine learning models. In addition, the DFT-GGA bandgaps and host densities (collected from the Materials Project database)⁴⁷ for the compounds in the two training datasets are provided and ML predictions for the Ce³⁺ activator’s 4f and 5d₁ levels for the perovskite oxide and elpasolite datasets are also reported.

The authors gratefully acknowledge support from the Laboratory Directed Research and Development (LDRD) program within the Los Alamos National Laboratory (LANL), an affirmative action equal opportunity employer, operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. DOE under Contract No. DE-AC52-06NA25396. Computational support for this work was provided by LANL’s high performance computing clusters.