DFT calculation of 229thorium-doped magnesium fluoride for nuclear laser spectroscopy

The 229thorium nucleus has an extremely low-energy isomeric state that could be manipulated with light in the vacuum ultraviolet (VUV) range. Recent measurements based on internal conversion electrons place the isomer energy at 8.28(17) eV (Seiferle B et al 2019), within the transmission window of large-band-gap materials, such as fluoride single crystals. Doping 229Th into VUV-transparent materials realizes a spectroscopy target with a high nuclei density and might form the basis of a solid-state nuclear clock. This paper presents a theoretical study of the optical properties of a thorium-doped MgF2 crystal. Using the Vienna Ab-initio Simulation Package, we perform density functional theory calculations of the electronic and optical properties of Th:MgF2. We determine whether thorium will be accepted as a dopant and identify the charge compensation mechanism and geometry. The simulations indicate, that the band gap of Th-doped MgF2 will be significantly reduced compared to undoped MgF2, below the expected 229Th isomer energy. This result is in striking contrast to a similar study performed for Th-doped CaF2 (Dessovic P et al 2014 J. Phys. Condens. Matter 26 105402).


The quest for the 229 Th isomer
Laser spectroscopy is used for many different applications in nuclear physics [3][4][5][6], but not for direct resonance excitation of nuclear transitions. The reason is that the excitation energies of most known excited nuclear states lie far above the range of modern lasers. The only known exception is a long-lived thorium-299 isomeric state known also as 'thorium isomer', whose energy is extremely close to the nuclear ground state. In 2007, an isomer energy of 7.8(5) eV was obtained from indirect measurements with a NASA x-ray microcalorimeter 4 Author to whom any correspondence should be addressed. [7,8]. Recent measurements of spectra of internal conversion electrons emitted in-flight during the decay of the excited nuclei of neutral 229 Th atoms shifted this value to 8.28 (17) eV [1]. This energy is within the reach of modern optical laser spectroscopy, and the isomer transition could serve as a clock transition in a nuclear frequency standard [9], which might reach an unprecedented uncertainty level of 10 −19 [10]. It would also be a very sensitive probe for testing the stability of fundamental constants [11,12]. Frequency shifts and broadenings of this transition in a solid-state environment might be used in studies of material properties, as is commonly done in Mössbauer spectroscopy [13]. It has also been shown that an ensemble of thorium nuclei doped into a transparent crystal may demonstrate superradiance with non-trivial temporal dependences [14], and may be used for constructing a nuclear gamma-ray laser [15].

The solid state approach and magnesium fluoride as a possible host crystal
The isomer energy may be accurately and unambiguously determined only in direct measurements, where either the energy of nuclear fluorescence photons is measured in direct radiative decay, or the excitation of the isomer using narrowband tuneable vacuum ultraviolet (VUV) radiation is detected. Because the transition is extremely weak, and the range of search is quite broad, one needs to address a high number of thorium nuclei simultaneously. Also, non-radiative decay channels (internal conversion) need to be suppressed, if they are not used for the detection, as proposed in [16]. These conditions can be attained by implantation of a macroscopic amount of thorium ions into some transparent crystal, as has been suggested by Peik and Tamm in 2003 [9]. Luminescence of thorium-doped crystals after synchrotron irradiation has been studied in [17,18], so far without detection of a signal. In the first case, LiSrAlF 6 has been used as a host crystal, with doping concentrations of thorium ions up to 4.1 × 10 17 cm −3 . The study allowed to exclude a thorium isomer transition with a vacuum lifetime between 1 s to 5600 s, at an isomer energy between 7.3 eV to 8.8 eV. In the latter study, a CaF 2 host crystal was doped with thorium with a concentration of 8.8 × 10 15 cm −3 . This study also allowed to exclude an isomer transition with a vacuum lifetime between 0.2 s to 1 s, at a transition energy between 7 eV to 9.5 eV. Both results of [17,18] are consistent with a recent work of Seiferle et al [1], and with recent theoretical estimations of the bare isomer lifetime [19], to be about 10 4 s.
Both experiments [17,18] did not succeed to detect the isomer signal, which is masked by the photo-and radioluminescence background. To enable a successful measurement, one needs to grow a new host crystal with reduced luminescence. This may be attained by improving the crystal growing procedure with a decreased amount of contaminations and defects and/or by choosing a new host crystal, with reduced radioluminescence in the vacuum ultraviolet range. To this end, MgF 2 appears to be a promising candidate. MgF 2 is an insulator with a band gap of about 12.4 eV [20] and it exhibits an exceptionally wide transparency range from a few μm into the VUV region. In comparison with CaF 2 , it shows less Cherenkov radiation because of the smaller refractive index, and its radioluminescence caused by heavier elements lies primarily in the longer wavelength range [21]. MgF 2 is a promising host crystal also for an experiment with 'optical pumping' of the thorium isomeric state via the 29.19 keV nuclear state [22], since it consists of lighter atoms and absorption depends strongly on mass. The penetration depth, i.e. the length over which the intensity of the 29.19 keV radiation drops by a factor of 1/e, is 5.27 mm for MgF 2 , 1.32 mm for CaF 2 and 0.325 mm for LiSrAlF 6 [23], where the density of LiSrAlF 6 was taken as 3.45gcm −3 [24].
MgF 2 crystallizes into a simple tetragonal cell (a = b = c) having a rutile structure, space group 136 [25][26][27], as shown in figure 1. The unit cell contains two formula units. Atomic positions are determined by the c/a ratio and the internal parame- ter u. The cations (Mg 2+ ) occupy the 2a and the anions (F + ) the 4f Wyckoff sites. Each cation is neighbored by six anions forming a distorted octahedron, with two neighbors at distance

Computer modelling
In the present paper, the electronic and structural properties of thorium-doped MgF 2 have been calculated using the Vienna Ab-initio Simulation Package (VASP, version 5.4.1) as devised by Kresse and Furthmüller [29,30]. To determine the equilibrium positions of the atoms in the crystal, VASP calculates the forces acting on the ions from the electronic structure. The movement of the heavier ions is then considered within classical mechanics, which is a consequence of the Born-Oppenheimer approximation [31], where the ionic movement can be decoupled from the electronic part, due to the large mass of the nuclei. For the relaxation procedure we either used a conjugate gradient [32] method or the scheme developed by Wood and Zunger [33]. VASP treats the effects of exchange and correlation within density functional theory (DFT) with different approximations. Most simulations in this paper have been performed using the strongly correlated and appropriately normed (SCAN) functional [34], a variant of a generalized gradient approximation (GGA). This functional has been confirmed to enhance results obtained by other GGAs at a comparable computational effort [35][36][37]. For some cases, we also performed calculations using other functionals and methods, in particular LDA (local density approximation), PBE (Perdew-Burke-Ernzerhof [38]), as well as the Hartree-Fock (HF) theory, the hybrid HSE (Heyd-Scuseria-Ernzerhof [39]) functional and the G 0 W 0approach [40]. A comparison of results obtained by different methods provides insight into the reliability and accuracy of the calculations. In this work, the most relevant physical quantity is the band gap. LDA and PBE are known to underestimate the band gap, while PBE typically outperforms LDA calculations in terms of accuracy, keeping the increase in computational time marginal. SCAN is a relatively new functional and so far, it is assumed to perform slightly better than the PBE. HSE is a hybrid functional where the exact HF exchange is mixed with the PBE functional. It usually predicts a more accurate band gap but is outperformed in required computation time by the previous functionals by several orders of magnitude, due to the demanding HF calculation. The G 0 W 0 approach is concerned with finding the Green's function and includes the calculation of the self energy. Since the first bands above the gap are explicitly calculated, the resulting band gaps are generally more precise but often only serve as a reference, due to the large computational requirements.

Modeling results
We calculate the electronic structure of MgF 2 to determine whether thorium will be accepted as a dopant and which charge compensation mechanism will occur. We consider 28 different variants of charge compensation, and estimate the size of the band gap for these configurations to see whether the resulting crystal stays transparent for the expected thorium isomer energy.
The energy cut-off of the plane waves determines the size of the corresponding basis set and is determined in accordance with E cutoff = 2 k 2 /2m e and is chosen to be 1300 eV, to resolve total energy differences well below 1 meV. The pseudopotentials for Mg, F, and Th are provided by the database within the VASP package.

Undoped crystal: band structure
To determine the best choice for the exchange-correlation functional with respect to chemical accuracy compared to computational time, we analyzed the electronic structure of MgF 2 . If the band structure computed with different functionals remains largely unchanged compared to more sophisticated approaches, we conclude that these functionals are suitable for an accurate simulation. The band structures are presented in figure 2.
For all cases compared, the general shape of the bands was rather similar, with the only difference being the value of the band gap, as shown in figure 2. For the purposes of this study, the best compromise between computational time and accuracy is the SCAN functional, as it improves on the PBE band gap, as shown in table 1, while being orders of magnitudes faster than HSE06 or G 0 W 0 .

Doped crystal: general consideration of charge compensation schemes
We assume a supercell based on the simple tetragonal primitive cell, where thorium is introduced as a dopant. In the experiments [17,18], doping concentrations are significantly lower than what is possible to simulate on modern clusters using DFT. Thus, we determine the doping concentration in the calculation by studying the interaction effects of the dopant sites due to periodic boundary conditions by varying the cell size. We performed calculations for super-cells of sizes up to 5 × 5 × 5 unit cells, corresponding to 750 atoms. The interaction strength was quantified by calculating the spatial shift of the atomic positions in proximity of the impurity in comparison to the largest cell. Displacements of atoms out- side of this sphere can be assumed to be negligible. If this shift is small, the impurities do not interact anymore and the simulation presents a valid approximation to realistic doping concentrations. As the computational demand is greatly affected by increasing the number of atoms, the 5 × 5 × 5 super-cell served as the reference calculation. We determined the average shift of atomic positions within a distance 3.11 Å from the substitutional replacement of magnesium with thorium. We observed a decrease of this average spatial displacement by about an order of magnitude for each super-cell size, i.e. 1.1 × 10 −1 Å, 2.0 × 10 −2 Å and 3.7 × 10 −3 Å for the 2 × 2 × 2, the 3 × 3 × 3 and the 4 × 4 × 4 super-cells respectively. Since this average change in positions decreases rapidly, we choose the 4 × 4 × 4 as a compromise between computational demand and chemical accuracy. This corresponds to a doping concentration of about 0.21 % or 1.9 × 10 20 cm −3 .
For a super-cell this large, we investigated whether it was feasible to sample the first Brillouin-zone with a single point, located at the origin Γ, to decrease the required CPU time. Compared to a 2 × 2 × 2 grid in reciprocal space, generated with the Monkhorst-Pack algorithm [42], the change in energy per atom turned out to be 0.65 meV. The sampling in k-space was therefore chosen to be the Γ-point only.
Doping can essentially take place in one of two ways, as shown in figure 3: substitutional doping, where the dopant replaces an atom in the original lattice site, and interstitial doping, where the dopant occupies a position in between the other atoms.
We calculate energies corresponding to these two doping types using the SCAN functional and we find that the substitutional doping is energetically preferred. Thorium has an oxidation number of +4 [43]. Because a multitude of different possibilities for charge compensation schemes can be  When thorium replaces a magnesium atom with an oxidation number of +2, two electrons will be available to form bonds with some other atoms which will be located in the vicinity of the defect. Such compensations can take place either via a single magnesium vacancy, or via two additional interstitial fluorines, which have an oxidation number of −1 respectively. Practically however, MgF 2 may be contaminated by oxygen with an oxidation number of −2 during the crystal growing process, which gives two more additional charge compensation schemes: single interstitial oxygen somewhere in between the lattice sites, or two substitutional oxygens replacing two fluorine atoms.
To compare different configurations it is required that the systems consist of the same particles in order to take all energy contributions into account. Therefore, the energies of F 2 and O 2 in the gas phase were calculated because they are either ejected into the gas phase or F 2 remains in the bulk in form of MgF 2 .
It is possible to calculate the energies of gases when using periodic boundary conditions by setting the lattice constant to a very high value, such that interactions between the different Table 2. Summary of the total energies of the elements needed for the charge compensation. For the calculation of the F 2 and O 2 molecules in the gas phase, the distance to the next atoms were determined by comparing the energy to a reference calculation with a distance of 20 Å. The magnesium and magnesium fluoride compounds were simulated in the crystal phase. Typically, binding energy is released when the crystal forms, which is also called heat or enthalpy of formation. To validate the results obtained for the individual atoms, the heat of formation for MgF 2 was compared with experimental data. The reference value is 11.69 eV per formula unit [44], whereas our simulation with SCAN results in 12.26 eV.
Th-doped MgF 2 with an assumed charge compensation via a single interstitial oxygen contains 127 Mg, 256 F, 1 Th and 1 O atom, and its characteristic energy is given by where E VASP is the relaxed total energy of the super-cell, E(MgF 2 ) is the energy per formula unit of MgF 2 and E(F 2 ) and E(O 2 ) are the energies of the respective molecules, as given in In turn, the super-cell with charge compensation via two interstitial fluorines contains 127 Mg, 258 F atom and 1 Th with the characteristic energy being calculated with Finally, the super-cell with charge compensation via single magnesium vacancy contains 126 Mg, 256 F atoms and 1 Th atom, with a characteristic energy of Characteristics of different charge-compensation schemes with different configurations of charge-compensating atoms are presented in table 3. The configurations are described by initial (unrelaxed) positions of charge-compensating ions/vacancies and are given in fractional units of the 1 × 1 × 1 unit cell with the thorium atom at the position (0.5, 0.5, 0.5).
Positions are rounded to 0.01 for the sake of brevity. Note that all the calculations were done using the SCAN functional, which underestimates the band gap, as it has been discussed in figure 2.
We find that the energetically most favorable configuration is the magnesium vacancy (configuration #26). Its calculated band gap is about 4.92 eV and is significantly lowered compared to pure MgF 2 due to the presence of thorium based electronic states. However, it is well known that DFT systematically underestimates the band gap [45,46]. An estimation of the magnitude of the true band gap can be obtained by comparing the band gap found experimentally in MgF 2 , which is 12.4 eV [20], with the one obtained for pure MgF 2 in the DFT calculation, which was determined to be 7.62 eV. This results in a simple scaling factor of ≈ 1.63 which can be applied to correct the calculated band gaps [47,48], denoted by Δ + , leading to about 8.01 eV. This is lower than the currently expected isomer energy of 8.28(17) eV [1]. One should point out, however, that such a simple scaling correction may also lead to residual error of about 1 eV, and the isomer energy may lie within the band gap. But even in such a case, laser spectroscopy will be very challenging near the absorption edge in optical transmission, broadened in accordance with the Urbach rule [49].
The second most favorable charge compensation type is +2F interstitial (configuration #22) and its characteristic energy is 0.8072 eV higher, which gives the relative Boltzmann factor exp −E 22 /k B T taken at the crystallization point T = 1536K of about 0.2%.
We also determine the electric field gradient acting on the nucleus. The quadrupole interaction of a nucleus with an external electric field is described by the quadrupole tensor and the electric field gradient [50]. The latter is given by the Hessian (rank two tensor of second derivatives) of the electrostatic  [2] of the ground and isomer states of the 229 Th nucleus in the presence of an electric field gradient calculated for the most favorable configuration #26. Here we use the spectroscopic quadrupole momenta of ground and isomer states equal to Q s = 3.11 eb [53] and Q m s = 1.74 eb [54] respectively. potential V that is external to the nucleus, evaluated at the nucleus position This matrix is symmetric and can therefore be diagonalized with the diagonal values V xx , V yy and V zz , which are conventionally ordered such that |V zz | |V yy | |V xx |. Furthermore, since it has trace zero, only two parameters are independent. Usually these are given by V zz and the asymmetry parameter For a more detailed discussion and notably for the case that the electron wave function has an overlap with the nucleus, see [51]. The electric field gradient can be calculated numerically [52] and the results for the most favorable configuration #26 are as follows: V zz = 75.095 V Å −2 , η = 0.624. Quadrupole shifts corresponding to such an electric field gradient can be calculated using the same method as in [2], their values are presented in figure 4.

Discussion
In previous work, Dessovic et al [2] investigated Th-doped CaF 2 crystals with DFT, which has a cubic crystal structure. CaF 2 is similar to MgF 2 in the sense that calcium is in the same group as magnesium, situated one row below in the periodic table. Therefore, the charge imbalance introduced to the crystal by doping is the same. However for CaF 2 , the preferred charge compensation scheme is not the vacancy of a neighboring cation but the +2F interstitial configuration. For this doping complex, the energy gap remains large enough for VUV excitation of the 229 Th isomer at the expected energy of 8.28 (17) eV. From this we conclude, that doping thorium into seemingly similar materials (different fluoride crystals) does not result in similar properties concerning the band gap. While CaF 2 seems a promising host material for VUV excitation, MgF 2 becomes intransparent at the expected energy of the 229 Th isomer. A similar analysis for LiSrAlF 6 is yet pending.
Recently, Barker et al [43] performed DFT calculations for thorium-doped MgF 2 . In this work, the preferred charge compensating mechanism was also determined to be the Mgvacancy. However a G 0 W 0 calculation showed a band gap of 11.9 eV which would be large enough for the 229 thorium isomer energy. This is in contrast to our result (8.01 eV). We speculate this to be a consequence of their smaller 2 × 2 × 2 unit cell. While this choice of a smaller unit cell allows the authors to perform a computationally demanding G 0 W 0 calculation based on LDA orbitals, it increases the effective doping concentration to 2.08% or 1.9 × 10 21 cm −3 and an interaction between dopants may not be neglected.
We attempt to repeat these calculations and also improve their accuracy by using the SCAN functional. To achieve convergence, SCAN requires a much larger energy cutoff in contrast to LDA, where it was possible to use only 400 eV. In agreement with Barker et al [43], we report that the most favorable configuration is the Mg-vacancy and our calculations also show, that the +2F interstitial charge compensation has almost no energy preference for a linear or bent structure. However, the energy difference between the Mg-vacancy and the +2F configuration is 2.7 eV in our results and thus slightly larger than in the previous study, which reported 2.5 eV. The band gap calculated with SCAN and subsequently scaled as discussed previously is determined to be 9.6 eV, which is indeed increased compared to the larger cell but still smaller than what was published by Barker et al using the G 0 W 0 method. Unfortunately, it is not possible to repeat the G 0 W 0 calculation in this study because of the enormous memory requirements due to the increase in electronic cutoff.
Nevertheless, we show that the doping concentration has an impact on size of the band gap. If the concentration is increased to about 1 in 50, the gap is larger compared to 1 in 750. This increase results in a qualitatively different proposition meaning that computational methods suggest that in the former case, the MgF 2 host crystal stays transparent to the Th nuclear transition energy while in the latter case it may not.

Conclusion
The DFT simulations in this work show that thorium, when doped into a MgF 2 crystal, substitutes a magnesium atom. It is determined that the most probable charge compensation scheme is a vacancy of a magnesium atom, situated closest to the impurity. The second-most favorable configuration is an interstitial position of two additional fluorine atoms, with a relative probability of formation during the crystallization process of about 0.2%. The calculations were done using the VASP package [29,30] and the SCAN [34] meta-GGA functional. The corresponding band gap is calculated to be 8.01 eV. In comparison with the most recent estimation for the transition energy of 8.28(17) eV [1], the band gap is predicted to be too small to permit optical VUV spectroscopy of the 229 Th isomer.
It has to be remarked, that the microscopic mechanism of optical absorption changes qualitatively when introducing Th-doping into the MgF 2 crystal. In the case of undoped MgF 2 , a photon exceeding the band gap energy can be absorbed, promoting an electron from the valence band to the conduction band (band-band transition). This leads to an absorption edge in the optical transmission, which is additionally broadened by excitons near the band gap edge due to electric fields in the crystal. Upon Th-doping, new localized electronic states emerge at the lower and upper band edges, that correspond to electron wave functions, localized at the Mg-vacancy or dopant site, respectively (color center or charge-transfer transition). Again the corresponding optical absorption bands are temperature broadened up to 0.5 eV [55]. Even when considering residual uncertainties in the quantitative predictions of DFT (related to their intrinsic model simplifications), we expect a strong overlap between the emerging absorption bands and the predicted isomer energy at 8.27 (17) eV, rendering nuclear laser spectroscopy very challenging if not impossible in the MgF 2 system. Some approaches involving nuclear excitation by electronic transition [56] and the electron-bridge [57] scheme may help to increase transition rates. However, their approach is sophisticated, especially near the edge of the band gap.