Electron Trap Depths in Cubic Lutetium Oxide Doped with Pr and Ti, Zr or Hf—From Ab Initio Multiconfigurational Calculations

We propose a universal approach to model intervalence charge transfer (IVCT) and metal-to-metal charge transfer (MMCT) transitions between ions in solids. The approach relies on already well-known and reliable ab initio RASSCF/CASPT2/RASSI-SO calculations for a series of emission center coordination geometries (restricted active space self-consistent field, complete active space second-order perturbation theory, and restricted active space state interaction with spin-orbit coupling). Embedding with ab initio model potentials (AIMPs) is used to represent the crystal lattice. We propose a way to construct the geometries via interpolation of the coordinates obtained using solid-state density functional theory (DFT) calculations for the structures where the activator metal is at specific oxidation (charge) states of interest. The approach thus takes the best of two worlds: the precision of the embedded cluster calculations (including localized excited states) and the geometries from DFT, where the effects of ionic radii mismatch (and eventual nearby defects) can be modeled explicitly. The method is applied to the Pr activator and Ti, Zr, Hf codopants in cubic Lu2O3, in which the said ions are used to obtain energy storage and thermoluminescence properties. Electron trap charging and discharging mechanisms (not involving a conduction band) are discussed in the context of the IVCT and MMCT role in them. Trap depths and trap quenching pathways are analyzed.


■ INTRODUCTION
Long-term energy storage and thermoluminescence in phosphors based on cubic Lu 2 O 3 :Pr can be achieved via the introduction of d metal codopants such as Hf, 1−4 Nb, 5 and Ti. 5,6 Under ultraviolet (UV) light irradiation, the Pr 3+ dopant is excited to its 4f 2 and 4f 1 5d 1 excited states. Some of those states overlap with the conduction band of the matrix material (Ia3̅ Lu 2 O 3 ). 7,8 A Pr 3+ → Pr 4+ photoionization is thus possible upon Pr 3+ f−d or f−f excitation. The formation of Pr 4+ was suggested by Wiatrowska and Zych 3 from the respective decrease of Pr 3+ f−d UV absorption in the X-ray-irradiated Lu 2 O 3 :Pr,Hf (i.e., upon high-energy excitation, some part of Pr 3+ were converted to Pr 4+ ). The d orbitals of the codopant (Hf) were assumed to trap the resulting electron, 2 although this was not the only assumed electron trapping mechanism. The same work showed that even without the codopant, Lu 2 O 3 :Pr can exhibit thermoluminescence�albeit very weak. X-ray absorption near-edge structure (XANES) X-ray absorption experiments 9 show the presence of both Pr 3+ and Pr 4+ in Lu 2 O 3 :Pr,Hf and Lu 2 O 3 :Pr, indicating that Pr 4+ is stable in these materials.
The following is a commonly accepted thermoluminescence/ energy storage mechanism in Lu 2 O 3 phosphors. Before the irradiation (in its state without the trapped electron), the codopant is assumed to be in its higher oxidation state (e.g., Hf 4+ ). After the excitation and subsequent Ln 3+ → Ln 4+ ionization (Ln = Pr, Tb), the electron released at the former 3+ Ln site is assumed to travel via a conduction band until it reaches a trap site. The trap sites are characterized by a separated energy level below a conduction band and are typically attributed to the codopant. After the excitation, the codopant is reduced (to Hf 3+ , for example) by the electron received from the conduction band. The reduced codopant remains in that oxidation state until it receives a portion of energy high enough to promote the electron back to the conduction band.
Depending on the codopant, the electron traps in Lu 2 O 3 :Pr can have substantially different depths. 5,6 For the shallow traps corresponding to the thermoluminescence glow curve peaks below 200°C, 2 visible afterglow is observed after UV irradiation of Lu 2 O 3 :Pr,Hf, for example. With the increase in the depth of the traps, permanent energy storage can be achieved. In the latter case, once the material is irradiated, it can forever store a portion of the excitation energy. Similarly, Lu 2 O 3 :Tb can be codoped with transition metal cations. Given the appropriate synthesis conditions, a Lu 2 O 3 :Tb phosphor with thermoluminescence properties can be made. From the experimental data collected for Lu 2 O 3 :Tb,M (M being Ti, 10,11 Zr, 11 Hf, 10−14 V, 15 Nb, 10 and Ta 16,17 ) and Lu 2 O 3 :Pr,M, 1−6 it occurs that the codopant is more crucial in defining the electron trap character. For instance, the most intense thermoluminescence glow peak is located at roughly the same temperatures 5 11 where Zr and Hf codopants resulted in practically identical thermoluminescence in Lu 2 O 3 :Tb,M. Between Pr and Tb, the Pr dopant was selected as the subject of this paper as it is significantly simpler to treat computationally due to the small number of valence electrons�namely, two for Pr 3+ and one for Pr 4+ .
Despite the vast experimental data, a particular mechanism of energy storage and trap depopulation in Lu 2 O 3 :Ln,M is not completely clarified. Albeit likely, it is not directly confirmed that all of the mentioned codopants act as electron traps (i.e., if the dopant cations get locally reduced to form metastable trap states). Moreover, recent studies for Lu 2 O 3 :Tb,M phosphors indicate that some of the depopulation processes likely involve the trapped electron excitation to the conduction band, 18 while the others do not�a direct transfer of the electron to the Tb 4+ / Tb 3+ recombination center occurs. 19 In this research, we have tackled the charging−discharging (electron trapping and release) mechanism using advanced post-Hartree−Fock correlated electron calculations (multiconfigurational and coupled cluster). Several practical issues were resolved, giving a rather straightforward and universal approach to estimating the trap depth for a metal-to-metal charge transfer (MMCT) mechanism of trap population and depopulation. In other words, we have considered a direct electron exchange between the dopant and the codopant without a conduction band involved. Ti, Zr, and Hf codopants were selected: they feature only two oxidation states each (+3 and +4) and correspond to a broad range of traps (Ti: deep; Zr, Hf: moderate and shallow). Additionally, as the [Ln 3+ , M 4+ ] and [Ln 4+ , M 3+ ] ionic pairs are symmetric in terms of local uncompensated charge (+1 at one of the dopants, zero at the other, both in Lu 3+ sites), no additional energy shifts (due to electron−hole attraction) are needed. 20 Identical traps corresponding to Lu 2 O 3 :Tb codoping with either Hf or Zr 11 were another reason to compare Zr and Hf electron trapping properties.
To model the thermoluminescence emission mechanism in Lu 2 O 3 :Pr,M, the appropriate emission levels of Pr 3+ had to be selected. In cubic lanthanide sesquioxides, Pr 3+ does not exhibit green emission from 3 P/ 1 I manifolds (the matter revised by Srivastava et al. 8 and references therein). Instead, red emission from 1 D 2 is observed. This phenomenon is explained by crossrelaxation by Pedroso et al. 9 Srivastava et al. explain it by the overlap of the 3 P and 1 I manifolds with the conduction band. 8 Both mechanisms result in efficient quenching of the 3 P/ 1 I manifold and the consequent emission from 1 D 2 . In the work by Kulesza et al., 4 a thermoluminescence emission spectrum of Lu 2 O 3 :Pr,Hf ceramics is shown in the 450−800 nm range, indicating only red emission peaks at 600−700 nm and no green emission. The Pr 3+ red thermoluminescence emission in the 600−700 nm range was attributed to transitions from the 1 D 2 manifold by Wiatrowska et al. 2 Kulesza et al. show 5 that both optically stimulated Pr 3+ emission (upon optical emptying of the previously charged traps) and regular UV−vis Pr 3+ emission spectra are very similar to those of the thermoluminescence Pr 3+ emission. In other words, no matter the population mechanism, the Pr 3+ emission in Lu 2 O 3 is always red.
Noteworthy, Toncelli et al. 21 attribute this emission to transitions from the 3 P manifold, which is not surprising. Using the levels from the said paper and all of the permutations between them, it is possible to find transitions from both 3 P and 1 D 2 (to the 3 H 3−4 levels) that match the 600−700 nm peaks, while the selection rules predict higher intensities for the spinallowed triplet−triplet transitions from 3 P. A similar conclusion can be made from energy levels obtained in this work�from energy differences alone, both 3 P and 1 D 2 levels can be the initial ones of the emission in question. Quenching processes must therefore define the observed emission. In summary, we are inclined to assume that 1 D 2 is the initial level of interest. However, we have also analyzed the trap depths for the 3 P/ 1 I population.

■ METHODS
The model of energy storage applied in this work is based on the following ideas. In its uncharged (inactivated, unperturbed) state, the thermoluminescence (storage) Lu 2 O 3 :Pr,M (M = Ti, Zr, Hf) phosphor contains (from a very localized, microscale, atom-based point of view) a Pr 3+ ion in an Lu 3+ site and an M 4+ ion in an Lu 3+ site. Upon excitation, Pr 3+ gets ionized to Pr 4+ , and the produced electron is transferred to the M 4+ ion�which is reduced to M 3+ . Provided there is a barrier for the reverse process, the system will stay in this metastable charged (activated, Pr 4+ + M 3+ ) state until the energy sufficient for the electron to overcome the barrier is delivered to the system. In other words, the electron is trapped at the M 3+ site. The model used does not include a conduction band and hence corresponds to a metal-to-metal charge transfer (MMCT) mechanism of the electron trapping/detrapping processes.
Dopant Ion Clusters and Pairs. The part of the Lu 2 O 3 :Pr,M (M = Ti, Zr, Hf) materials that was modeled explicitly was the dopant cation and the oxygens of its immediate surroundings. The rest of the lattice was represented using embedding potentials. Such an approach assumes that the corresponding dopant clusters in the real material are separated by a certain distance�such that any deformations in the geometry of one cluster do not affect the other cluster. As the thermoluminescence materials of interest are characterized by a rather low dopant content, 5 such an assumption is entirely acceptable. Consequently, independent calculations for each cluster geometry were performed. We have also considered the case of the two dopants occupying nearest neighbor sites. That requires a large embedded cluster containing both dopant ions. Such a case is methodologically quite different from what is presented below and is much more complex and extended. It was thus not included in the current paper.
To  6 8− and MO 6 9− separated clusters. Please note that the cluster−cluster coupling energy is not calculated: the energy (diabatic) of a cluster pair is a sum of the separate noninteracting cluster adiabatic energies at a particular reaction coordinate. This means that diabatic electronic energy surfaces will not provide quantitative energy barriers in the case of avoided crossings of the cluster pair electronic surfaces. The surfaces, however, let us follow the chemical character of the cluster pair electronic states along the reaction coordinate and estimate energy barriers from diabatic crossings.
Potential Energy Surfaces. In the assumed mechanism of ion−ion interactions (MMCT in particular), two ions change their oxidation states upon electron transfer. That, in turn, must result in non-negligible changes in bond length between the ions and the surrounding oxygens. When a geometry of a single cluster is changed in correspondence with a vibrational mode (along a linear reaction coordinate), the electronic energy of the cluster changes. This change can be plotted as a function of the reaction coordinate, which forms a potential energy curve. The process of electron transfer involves two separate clusters and thus depends on two independent reaction coordinates. In a Let us denote the diabatic electronic energies of the ion pair as z = E(x, y), where the Pr 4+ /Pr 3+ reaction coordinate is used as x, and the M 3+ /M 4+ reaction coordinate is used as y. In particular, the diabatic potential energy surface of the activated/charged . The surfaces intersect along a line that (usually) has a minimum (at the proximity of the studied system minima). The difference between the minimum of the E I surface of the charged system [Pr 4+ , M 3+ ] and the minimum of the E I /E II surface intersection is taken as the trap depth in this model.
Cluster Geometries and Pseudomodes. In either thermal detrapping or a more general MMCT process, the 3+ and 4+ equilibrium geometries of a particular ion must transform into each other�i.e., atoms of the ion's surroundings have to move in response to the change in the charge state of the said ion. Thus, a coordinate path transforming the XO 6 9− equilibrium geometry into the XO 6 8− equilibrium geometry (X = Pr, Ti, Zr, Hf) must exist. We have proposed making a linear interpolation between the X 3+ /X 4+ equilibrium geometries to exemplify this coordinate path. As stated before, a z = E(x, y) potential energy surface is a function of two independent reaction coordinates x and y: the former corresponding to the Pr 3+ /Pr 4+ coordinate path, and the latter corresponding to the M 3+ /M 4+ coordinate path (M = Ti, Zr, Hf).
To construct the respective coordinate path, the initial and final geometries of the clusters were obtained using density functional theory (DFT, the next subsection). Having those, it is possible to construct a transformation numerically analogous to a vibrational mode, such that it transforms the cluster geometries (of the same atom at different oxidation states) into each other. We have called such a transformation a pseudomode. It is assumed to represent a thermally induced transition and, thus, might be a part of a more extensive (nonlocal, phonon) vibrational mode. The utilized pseudomodes have produced smooth (parabola-shaped) potential energy curves at the RASSCF level of theory.
Given a Cartesian coordinate matrix G 1 (representing, for example, equilibrium PrO 6 8− geometry, Pr 4+ ) and a Cartesian coordinate matrix G 2 (representing, for example, equilibrium PrO 6 9− geometry, Pr 3+ ), the Cartesian displacement matrix of the pseudomode is given by D = G 2 − G 1 . The interpolated geometries G u for a potential energy curve that uses the pseudomode are created using a scalar parameter u (a pseudomode coordinate): G u = uD + G 1 . Setting u = 1 gives G 2 , u = 0 corresponds to G 1 , and 1 > u > 0 corresponds to the intermediate (interpolated) geometries, while the other values of u give the extrapolated geometries. The same principle was used for M 3+ and M 4+ geometries.
Having a pair of Pr and M ions and the 3+ and 4+ equilibrium geometries for each of them (G PrO 6 9 , G PrO 6 8 ; G MO 6 8 , G , MO 6 9− clusters, respectively. The diabatic potential energy surfaces of activated/deactivated systems can now be expressed using parameters u Pr and u M of the respective pseudomodes instead of the arbitrary x and y: The energies in the equations above were obtained for a series of interpolated geometries G Pr and G M , using ab initio calculations described in the following section. Each of the obtained potential energy curves (of the individual clusters) was fitted with a cubic polynomial to get the respective analytical representation of the curves. The intersection of the E I and E II surfaces (defined in eqs 3 and 4, respectively) was found analytically using the cubic roots formula. The details are provided in the Supporting Information (in the Python script that does the analysis and plotting). The intersection of the two surfaces is a nonplanar three-dimensional curve. The minimum of the intersection curve was used to estimate the trap depths and other energetic parameters.
There are several factors motivating the transformation of cluster geometries along pseudomodes. The approach involving intersecting diabatic surfaces has been used before in intervalence 22,23 and metal-to-metal 20,24,25 charge transfer studies. However, those studies featured highly symmetric CaF 2 -type lattices, where activator site geometry optimization was not required: vibrational deformations were addressed via simple proportional scaling of the whole cube surrounding the metal ion, corresponding to a full-symmetric (breathing) vibrational mode of the activator site. The structure of cubic lutetium oxide is more complex: 26 there are two six-coordinated sites of C 3i and C 2 local symmetries, both of which look like distorted cubes with two unoccupied vertex sites.
A metal-to-metal charge transfer process means that every involved ion changes its oxidation state. For the dopant ions, there are two charge states (before and after) and two site symmetries to be considered. For each of those, the equilibrium geometries of the coordination polyhedron can be different in terms of both bond length and bond angles. Unlike CaF 2 -type lattices, the individual before and after geometries do not necessarily transform into each other via scaling of the bond lengths.
For every ion in question and for every site symmetry of the c-Lu 2 O 3 lattice, normal vibrational modes were obtained and analyzed for the embedded clusters. It turned out that for the same ion and site symmetry, none of the 4+ ion cluster modes can transform (morph) the 4+ ion equilibrium geometry into the 3+ ion equilibrium geometry. The same is true for the 3+ ion cluster modes and 4+ ion cluster geometry. In other words, for the same ion and site symmetry, the normal mode vibrational coordinates of the 3+ ion cluster do not intersect with the normal mode vibrational coordinates of the 4+ ion cluster at any point. The geometries produced by the respective deformations do get close to each other in some cases but are never identical. The pseudomode coordinate path solves this issue, making two arbitrary geometries transform into each other by construction. ) were obtained using planewave density functional theory calculation (DFT) on a doped unit cell of Lu 2 O 3 , doing a full cell relaxation (for both 3+ and 4+ dopants). The input geometry from the paper by Zeler et al. 26 was used. GBRV 27 ultrasoft pseudopotentials for Ti, Zr, and Hf and Topasakal and Wentzcovitch 28 projector augmented wave (PAW) potentials for Pr were used. PBEsol 29 exchangecorrelation functional was utilized. The code was the PW module of Quantum Espresso. 30,31 The detailed description of the code parameters and functional selection is given in the Supporting Information.
For the Hf-doped systems, two independent sets of calculations were performed, with the the only difference between the two being the Hf pseudopotential. Hafnium was regarded as a problematic element by Garrity et al.: 27 namely, a universal Hf 0 ultrasoft pseudopotential did not perform well as Hf 4+ in different oxides and, thus, a dedicated Hf 4+ ultrasoft pseudopotential has been introduced. The pseudopotentials differ, among others, by their valence shells: both have occupied 5s and 5p orbitals and empty 6s and 6p orbitals; the Hf 0 pseudopotential was made with two electrons in 5d orbitals, while the Hf 4+ pseudopotential has empty 5d orbitals and additional 5f empty orbitals. In our case, an Hf 3+ ion in an oxide should be more similar to an Hf 4+ ion rather than a neutral Hf.
However, just to be safer, we made two sets of calculations using each of the pseudopotentials. As a result, two slightly different cell geometries were obtained, and the two sets of data are present.
Ab Initio Calculations. The following post-Hartree−Fock correlated electron calculations were performed with the OpenMOLCAS 32,33 software package. The clusters were placed in the embedding that represented the unperturbed lattice of c-Lu 2 O 3 : a 12 Å layer of ab initio model potentials (AIMPs) 34−36 surrounded the MO 6 clusters, followed by a layer of point charges with the external radius of about 77 Å. The point charges were optimized to minimize the values of electric multipoles (of orders 2, 3 and 4) at the cluster atoms using Lattgen code. 37 This kind of embedding was used to optimize the AIMPs in a selfconsistent embedding ion procedure. 38 The ready-to-use embedding files are given as Supporting Information. For each specific geometry of each cluster, the calculations began with relativistic Douglas−Kroll−Hess 39 integrals (order of Hamiltonian 2, order of properties 2, SEWARD code 40,41 ). Selfconsistent field calculations followed (SCF code, 42−46 Hartree− Fock, without the active electrons, i.e., on clusters containing Pr 5+ , Ti 4+ , Zr 4+ , Hf 4+ : PrO 6 7− and MO 6 8− ). The basis set was ANO-RCC triple-ζ. 47−49 The next step depended on the metal cation of the cluster. The potential energy curves for the ground state of the M dopants were calculated using the CCSD(T) 50 calculation on top of a single-root ground state restricted active space (RASSCF 51−53 ) wave function (using single-orbital active space, effectively an SCF calculation). For the M 3+ systems, full T2 and T1 spin adaptation (according to Raghavachari et al. 54 ) was utilized. In the noniterative triples procedure, the denominators were the diagonal Fock matrix elements (closedshell M 4+ ) or the orbital energies (open-shell M 3+ ).
For Pr clusters, state-average RASSCF calculations followed the Hartree−Fock step. The active orbitals (RAS2 orbitals in MOLCAS notation) were molecular orbitals with a predominant Pr 4f, 6s and 5d character. In the case of Pr 3+ /PrO 6 9− , two electrons populated the active space. In the case of Pr 4+ /PrO 6 8− , one electron populated the active space. All possible roots were included (i.e., the number of roots [states] for the state-average calculation was equal to the number of configuration state functions). For Pr 3+ , singlet and triplet states were considered. For Pr 4+ , doublet states were calculated. Next, restricted active space second-order perturbation theory (RASPT2 55−57 ) calculations followed. For Pr 3+ , only the most significant roots (all of the 4f 2 and 4f 1 5d 1 states) were included due to calculation stability issues. What is noteworthy is that the zeroth-order wave function reference weights in the RASPT2 calculations were typically 0.64−0.67. The frozen core orbitals were O 1s; Hf, Pr 1−4s, 2−4p, 3d, 4d; Ti 1−2s, 2p; Zr 1−3s, 2−3p, 3d�as recommended by the basis set authors. 47−49 The IPEA shift was set to zero. The imaginary shift 58 of 0.3 was used for Pr (required to get smooth Pr 3+ potential energy curves at the RASPT2 level of theory). Finally, RASSI-SO 59,60 calculations were used to include the effects of spin−orbit coupling and to couple the states of different irreducible representations and spin multiplicities: for the C 2 sites, representations A and B were coupled, while for the C i site, representations A g and A u were mixed.
Pr 3+ Electronic State Selection. The potential energy surfaces were represented by the following: E I (u Pr , u M ) (eq 3) is calculated as the energy of PrO 6 9− selected RASSI-SO state plus MO 6 8− CCSDT energy; E II (u Pr , u M ) (eq 4) is calculated as the energy of PrO 6 8− RASSI-SO ground state plus MO 6 9− CCSDT ]. For one of the quenching cases, the Pr 3+ ground state was picked. For other processes, the most relevant excited states were selected. One of them was the top of the 1 G manifold, as deduced from the energy gaps between states. These two cases corresponded to the final Pr 3+ states that resulted in no Pr 3+ red emission (i.e., the nonradiative loss of the trapped electron energy). To get the trap depth corresponding to the experimentally observed red emission in the 600−700 nm range (described by, e.g., Wiatrowska and Zych 2 ), the respective excited state had to be selected. As shown in the Introduction, such a selection is not straightforward. We selected the lowest level of the 1 D 2 manifold�it was also the first level above the 1 G manifold. That was the lowest Pr 3+ level from which the red emission is still possible. The level is referred to below as 1 D 2 . For comparison, we have also selected a higher level from which green emission should be possible. For the C 2 site, we picked a level with relatively high spontaneous emission coefficients (calculated with RASSI-SO). By the singlet−triplet character, it is likely a 3 P level. For the C 3i site, however, the coefficients were exactly zero due to inversion symmetry. We have consequently picked the first level above the 1 D 2 manifold, which is mostly singlet and is, thus, likely a 1 I level.
■ RESULTS AND DISCUSSION Calculated Energy Levels. The calculated energy levels (that included spin−orbit coupling, the RASSI-SO states) for the two Pr 3+ sites, the NIST Pr 3+ free-ion energy levels 61 21 ), while the triplet/singlet character is a good match to that of the NIST levels. 61 The calculated energy levels are also in accord with the ones reported by Pascual et al., 62 where the gap between the 1 D and 3 P/ 1 I manifolds is noticeably larger for Pr 3+ in the C 3i site than that of the C 2 Pr 3+ (Figure 1).
Energy Barriers: Error Estimation. It is worth mentioning that the calculated trap depths (and other energetic parameters) should correspond to the actual thermally driven processes in the real materials�from the purely physical-chemical microscopic perspective. We do not know for a fact that the pseudomode coordinate is the optimal path to the barrier. It is, however, a necessary approximation and very low-cost at that (compared to a hypothetical search for the lowest-energy path).
Certain unavoidable errors always emerge from the method's intrinsic limitations. In our particular case, the energy surfaces are diabatic (i.e., there is no cluster−cluster interaction explicitly modeled). Basically, all the crossings between diabatic energy surfaces of different electronic states of a cluster−cluster pair are allowed. We omit a description of the avoided crossings and conical intersections, which would be present in the adiabatic representation of the cluster−cluster pair. Therefore, our results are upper boundary estimates of the energy barrier heights.
On the other hand, even if the calculated trap depths were error-free, their values do not necessarily translate to the macroscopic thermoluminescence kinetics equations. 63 In the latter, the trap depth has a certain fit parameter flavor to it: it depends on the selected kinetics order, to say the least. Accordingly, trap depths from these calculations should be handled with care in the context of their comparison to the experimental trap depths.
Metal-to-Metal Charge Transfer and Localized Trapping. Each pair of the intersecting surfaces in question has its own unique look and properties. However, they are all similar in principle and feature two paraboloid-like surfaces and an intersection line. For clarity, we show only one of them in Figure 2, where the features are the most evident. Any curious reader is invited to use the supplementary data files and the attached Python script to view the respective interactive 3D plots.
The calculated trap depths, optical activation energy gaps and other energetic properties for the dopant pairs and local symmetries in question are summarized in Table 1 , u M )) potential energy surfaces. In Figure  2, that would be a transition from the red point to the purple I (the C 3i site) or 3 P (the C 2 site) levels.
Using the same principle, trap quench barriers (QBs) were obtained. The QB is the barrier of a process that would detrap the electron and result in no visible emission from Pr 3+ . In Figure  2, the QB would also be the energy of a transition from the red point to the purple point along the red surface�although this time the blue surface represents Pr 3+ in one of its nonemitting levels. If the quench barrier is low, filling the trap would result in fast thermal relaxation to [Pr 3+ , M 4+ ], and thus, the trap will not contribute to the energy storage and thermoluminescence. In other words, QBs are the barriers for thermally induced  Figure 2). Another kind of quench barrier featured Pr 3+ in the upper level of its 1 G 4 excited manifold. There can be no visible Pr 3+ emission from that level or any of the below levels, meaning that the stored energy is nonradiatively lost.
For the dopant pairs in question, an electron can also be transferred in a vertical (surface-to-surface) transition (straight lines in Figure 2), in which the geometries of the two pairs are the same (i.e., the pseudomode coordinates u Pr and u M must be the same for both surfaces in order for the vertical transition to happen). In time scales of those transitions, atom positions are considered static. The optical trap population gap (OTPG, the energy gap of an MMCT absorption that results in a filled trap) is the energy difference of a vertical transition from the [Pr 3+ , M 4+ ] energy minimum (E II minimum, Pr 3+ in its ground state) to the [Pr 4+ , M 3+ ] system at the same geometry coordinates. In other words, OTPGs are the energies of optical absorptioninduced metal-to-metal charge transfer transitions that result in the filled traps upon optical excitation.
Quite a few of the OPTG values in Table 1 lie in the 250−320 nm range, which is used to charge the traps in Lu 2 O 3 -based thermoluminescence materials. It is therefore possible that MMCT processes contribute to optical absorption: in Lu 2 O 3 :Pr, only one 250−300 nm broad band is present in photoluminescence excitation spectra, 9 while in Lu 2 O 3 :Pr,Ti there are two bands, at about 230−270 nm and about 300−350 nm. 6 Another form of a vertical transition is (non)radiative relaxation, provided that the final state potential energy surface is below the initial state potential energy surface at the configurational coordinates of the initial state minimum. The  ] pair is fully stable (a thermodynamical minimum, not metastable)�the latter should not be considered a trap state. Pr 3+ would be oxidized by M 4+ , provided that both are present in the sites corresponding to the described situation. This is always the case when Pr is in the C 3i site and the M codopant is in the C 2 site (see Table 1 In the case of the Ti codopant, we do not observe the expected experimental trap depth of 1.8−2.5 eV. 5,6,11,19 The calculated depth is either very low (0.033 eV Ti C 2 and 0.91 Ti C 3i ; Pr C 2 ), or the quench barriers are almost 0 eV (both Pr and Ti in the C 2 sites, or both ions in the C 3i sites). The [Ti C 2 , Pr C 3i ] system exhibits unidirectional oxidation and a large trap depth (3.5 eV). Thus, Ti either does not work as an electron trap in Lu 2 O 3 or (which is more likely) the detrapping mechanism must involve a conduction band. With the latter, the effect of quenching processes described in this work can be mitigated via an increase of the dopant−dopant distance�which is in line with the experimental low dopant concentrations (e.g., Kulesza et al. 5 18,11 This result is in line with the diabatic surface intersections providing "upper limit" estimates of the trap depth. Similarly to the Hf codopant, in the case of the Zr codopant, the same pattern is observed in the dependence of the energetic parameters (trap depths, quench barriers, optical and quenching gaps) on the site symmetries. The trap depths are too high (with respect to the experimental value, 1.44 eV 11 ) or the quench barriers are low. An interesting and noteworthy fact is that the calculated depths of the Zr-based traps are noticeably different from the respective values for Hf ( Table 1). The differences for the deep trap depths are in the range of 0.5−1.5 eV for [Zr/Hf C 2 , Pr C 3i ] and [Zr/Hf C 3i , Pr C 2 ] systems. Such differences should be clearly distinguishable from the thermoluminescence glow curves. However, Sojka et al. indicate 11 that Zr and Hf codopants result in identical glow curves in Lu 2 O 3 :Tb,M (M is Zr or Hf). According to the results presented here, the glow curves in the mentioned paper are unlikely to correspond to the [Pr 4+ , M 3+ ] → [ 1 D 2 Pr 3+ , M 4+ ] electron trapping and release mechanism�at least without the participation of a conduction band. If the mechanism was such, the trap depths and glow curves would not have been identical. This conclusion is quite important, and it can be argued that the pseudopotentials used to obtain the geometries might be a source of an error of some sort. We have consequently performed a calculation on the Zr dopant in the Hf geometry. Yet again, distinctly different trap depths were obtained, indicating that the difference originates from the Zr and Hf physical-chemical properties and not solely from the site geometries. The two elements are clearly not identical in their M 3+ −M 4+ transitions. It is hence reasonable to expect that with the conduction band participation, Zr and Hf would also exhibit different energetics.
We would like to point out that the presented calculations correspond to the electron trapping and detrapping that do not involve a conduction band. The conduction band cases can be modeled using the same approach with Lu 3+ /Lu 2+ instead of Pr.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article That problem, however, is more complex (in several aspects) and deserves a dedicated paper. Configuration Diagrams. The data in Table 1 can be visualized using three-dimensional plots ( Figure 2). However, such 3D plots are only readable and clear with two surfaces and one intersection line. Adding more surfaces makes the plot overwhelmingly complex. Additionally, a traditional way to visualize transitions is configurational diagrams (Figure 3), which are two-dimensional. Constructing the 2D diagrams from the 3D plots in question has its challenges. The points of interest in the 3D plots are the surfaces minima (E I and E II minima) and the minimum of the surface intersection line (Figure 2). A transition from the surface minimum (e.g., E I minimum) to the E I /E II surface-surface intersection minimum can take different paths that lie on the surface (e.g., E I ). The energy difference between these two points is the barrier energy, no matter the path. To visualize the transition, we do not need to know the exact path as well. Thus, for each surface, we have chosen a path produced by the intersection of that surface by a vertical plane that contains both the surface minimum and the E I /E II surfacesurface intersection line minimum. This vertical plane is taken as the 2D graph (image) plane, and the respective path is the configurational diagram in the graph.
In the studied ion pairs, for every elementary MMCT, there is always the charged and uncharged state. Those are, for example, [Pr 4+ , M 3+ ] (charged, filled trap, E I , the red surface in Figure 2) and [Pr 3+ , M 4+ ] (uncharged, empty trap, E II , the blue surface in Figure 2). For two surfaces, there is a minimum for each and one intersection minimum (Figure 2). These three points do not have to lie in the same vertical plane�in our results, they do not. Albeit the intersection points are shared for the charged and uncharged states, individual vertical planes and individual paths are required for each surface. Moreover, we have considered several states of Pr 3+ . For a particular pair of Pr and codopant in the sites of particular symmetries, there is a separate surface corresponding to each electronic state of Pr 3+ . Accordingly, for each dopant pair and site symmetry, we have independently constructed the individual transition paths and put them on the same 3D plot. The diagrams in the subplots of Figure 3 are produced as side views of such 3D plots. The energy axis and the Pr average bond length (and pseudomode coordinate) axis are in the screen (2D graph) plane, and the codopant (Ti/Zr/Hf) geometry axis is perpendicular to the screen plane. Basically, for an individual 3D curve in the form of u u E , , Pr data columns of the same length, column E u u , M Pr was plotted against column u Pr . The one-to-one correspondence of the data was provided by the vertical cut of the surface.
The resulting plots are shown in Figure 3. The selected dopant, codopant and site symmetry combinations are shown. For the C 2 Pr and C 3i M, the charged state (trapped electron) curves lie much higher than the Pr 3+ , M 4+ curves. In the C 3i Pr and C 2 M, the trapped electron curve minima lie below the Pr 3+ , M 4+ minima. The Hf cases are qualitatively very similar to the respective Zr cases for the shown curves. From the curves shown  Table 1 are provided. Note that in panel (a), the red curve minimum almost overlaps with another curve, and the minimal QB is close to zero. The calculated RASSI-SO Pr 3+ levels (in the corresponding sites) are shown as gray dashes.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article in Figure 3, it is clear that the Ti traps are much deeper than the Zr traps. At the same time, quench barriers are clearly low, while vertical relaxation processes are possible. This confirms the previously formulated conclusions about the traps being prone to quenching in the Lu 2 O 3 :Pr,M(IV) systems (Table 1 and its  description). Intervalence Charge Transfer: Carrier Migration. The previous section describes electron transfers that are traditionally considered in the context of thermoluminescence in Lu 2 O 3 materials, namely, electron trapping on a codopant (e.g., ref 5) . Another possibility is intervalence charge transfer (IVCT), 64 where an electron is transferred between the ions of the same chemical element that differ in oxidation state. With the site symmetries in Lu 2 O 3 , those would be  Table 2 summarizes the energetics of the M−M IVCTs in question. When the 3+ ion (donor) and the 4+ ion (acceptor) are of the same site symmetry, the electron transfer barriers are quite low, in the 0.1−0.5 eV range. This indicates that the trapped electrons can easily migrate among identical sites at temperatures much lower than room temperature. This property is another phenomenon related to the very low dopant concentrations required to achieve efficient energy storage in Pr-doped lutetium oxide. 3 On the contrary, when the site symmetries are different, there is a preferred direction of the intervalence charge transfer. 2 (M is Zr or Hf) and 6−7 eV for the reverse processes. As the lower of the two barriers is still quite high, the electron transfer between Zr or Hf ions of different site symmetries will not occur at room temperature. Comparing the Hf−Hf barriers to the Zr−Zr barriers in the respective site symmetries, we can observe that the values are similar but not identical. This supports the previous conclusion of Zr and Hf exhibiting different trap depths (at the end of Table  1 discussion).
The same intervalence charge transfer analysis has been applied to Pr. However, excited states were considered in that case, resulting in a more complex map of possible transitions. The analyzed processes concern an electron transfer between a certain Pr ion and another ′Pr ion. The cases where both Pr and ′Pr ions are of the same site symmetry are shown in Figure 4  processes, the transitions involving the same excited state of Pr 3+ on both sides were characterized by rather small thermalization barriers in the range of 0.3−0.5 eV ( The reverse process (i.e., the intervalence charge transfer relaxation) can take numerous pathways originating from the numerous Pr 3+ 4f levels�depending on the intersection positions between the states, as visualized in Figure 4. It is clear from Figure 4 that the [Pr 3+ , Pr 4+ ] potential energy curves corresponding to Pr 3+ in the lowest 1 D 2 level and the levels above intersect with the [Pr 4+ , ′Pr 3+ ] potential energy curves corresponding to Pr 3+ in some of its lower (nonemitting) levels. The barriers are low, and the respective quenching of Pr 3+ emitting levels via IVCT to a Pr 4+ should be efficient. This is one of the possible mechanisms for temperature-dependent crossrelaxation. However, as the excited states of lanthanides are usually long-lived, the quenching processes do not exclude a Pr 3+ 4f−4f emission. As the barriers are different for the Pr ions of C 2 and C 3i local symmetries (Table 3), there might be a difference in concentration quenching of the respective sites. Such an effect has been described by Bolek et al. 17 Table 3 and Table 2 also provide another clue for why the dopant concentrations required for efficient thermoluminescence 3 are usually low. The barriers for intervalence hopping of the electron from site to site are low, meaning that the electron is likely to either avoid recombination   (Figure 4). Note that both Pr 3+ and Pr 4+ are present in Lu 2 O 3 , 9 and, given our results, the 3+ sites should preferably be the C 2 sites. This also agrees with the results presented by Kulesza  process is characterized by a moderate barrier of 1.33 eV, meaning that, at room temperature, some 2 would accordingly provide additional (likely broad) absorption/excitation bands in blue and near-UV�at the wavelengths typically used in the Lu 2 O 3 excitation for thermoluminescence (250−320 nm; 10 our result might overestimate the transition energies). With the presented results, we emphasize that IVCT absorption in doped Lu 2 O 3 is expected to lie at about the same range of energies as the lanthanide f−d absorption. The f−d absorption corresponds to allowed transitions, while the IVCT absorption is expected to have lower intensity; 64 the former is hence more likely to be dominant.

■ CONCLUSIONS
In this paper, we have used DFT-derived geometries of metal clusters to construct a vibrational pseusomode coordination path, providing a vibration-like transformation between the coordination geometries of ions at different oxidation states. The paths were used as configurational coordinates to construct potential energy surfaces of ion cluster pairs and estimate metalto-metal and intervalence charge transfer (MMCT and IVCT) energetics for selected dopants in Lu 2 O 3 , namely, Pr, Ti, Zr and Hf. The presented approach can be used with any kind of crystal and site symmetries. As for the particular values of the trap depths, the calculated MMCT trap depths values do not correspond well to the experimental values from the previous topical publications. This is not surprising, as the experimental glow curves are interpreted with the mechanisms involving electron transfers via a conduction band. Our results do not exclude other MMCT detrapping mechanisms that do not involve the conduction band: more defects and dopant coordination geometries can be involved. Our model does not give a trap depth with respect to the conduction band, but it clearly describes the optical charging and relaxation of the trap states as well as their quenching.
Another important conclusion regards similarities and differences in the trap depths and other barriers for the processes that involve Zr and Hf codopants. According to the previously published experimental data, Zr and Hf codoping corresponds to identical glow curves in Lu 2 O 3 :Ln,M. 11 In the presented results, the Zr ions exhibit a pattern in the barrier values for various processes similar to that of the Hf ions�the particular values, however, are distinctly different for the two. In other words, according to the presented calculations, Zr and Hf cannot result in identical trap parameters. The two elements are clearly not identical in their M 3+ −M 4+ transitions. Consequently, local conversion of Zr/Hf ions from 4+ to 3+ is unlikely to be the electron trapping mechanism with those codopants. This conclusion is also in line with the DFT results for Hf in Lu 2 O 3 . 66 The Pr−Pr intervalence charge transfer might play an important role in the electron trapping and emission properties of Lu 2 O 3 :Pr. In particular, our results indicate that Pr 3+ ions are less stable in C 3i sites than in the C 2 sites. If Pr in one of the sites had to be 4+, that would rather be the C 3i site. An IVCT process of the excited acceptor formation� The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpca.2c07979.
Two Python scripts and a set of files with total energies that can be used to construct the surfaces and the configuration diagrams, including instructions on how to use the scripts (ZIP) DFT calculation details, cluster geometries, embedding AIMP files, and configuration diagrams for all Pr−M pairs (PDF) The Journal of Physical Chemistry A pubs.acs.org/JPCA Article