Ab-initio prediction of fast non-equilibrium transport of nascent polarons in SrI$_2$: A key to high-performance scintillation

The excellent light yield proportionality of europium-doped strontium iodide (SrI$_2$:Eu) has resulted in state-of-the-art $\gamma$-ray detectors with remarkably high energy resolution, far exceeding that of most halide compounds. In this class of materials the formation of self-trapped hole polarons is very common. However, polaron formation is usually expected to limit carrier mobilities and has been associated with poor scintillator light-yield proportionality and resolution. Here, using a recently developed first-principles method, we perform an unprecedented study of polaron transport in SrI$_2$, both for equilibrium polarons, as well as nascent polarons immediately following a self-trapping event. As a result, we put forward a rationale for the unexpected high energy resolution of SrI$_2$. We identify nine stable hole polaron configurations, which consist of dimerized iodine pairs with polaron binding energies of up to 0.5\,eV. They are connected by a complex potential energy landscape that comprises 66 unique nearest-neighbor migration paths. \textit{Ab initio} molecular dynamics simulations reveal that a large fraction of polarons is born into configurations that migrate practically barrier free at room temperature. As a result, carriers created during $\gamma$-irradiation can quickly diffuse away reducing the chance for non-linear recombination reactions that are the primary culprits for non-proportionality and resolution reduction. This leads to the conclusion that the flat, albeit complex, landscape for polaron migration in SrI$_2$ is key for understanding its outstanding performance. This insight provides important guidance not only for the future development of high-performance scintillators but also of other materials, for which large polaron mobilities are crucial such as batteries and solid state ionic conductors.


INTRODUCTION
A small, or self-trapped, polaron is a quasi-particle consisting of a localised charge carrier that has strongly polarised its immediate surrounding lattice. Polaron formation and transport are crucial for understanding important quantum-mechanical processes occurring in scintillations, 1 as well as batteries, 2 and solid ion conductors. The polaron can be centred on a specific atom, e.g., hole and electron polarons in TiO 2 or SrTiO 3,3,4 or atomic bonds, e.g., in the case of V K -centres in alkali halides. 5 The latter consists of a hole localised on a dimerised nearest-neighbour halide-ion pair. For small polarons such as these, the adiabatic potential energy landscape comprises potential wells that are usually o1 eV deep (corresponding to the polaron-binding energy) and separated by energy barriers that are generally lower than the binding energy. As a result, polaron migration typically proceeds by a hopping mechanism without delocalisation.
As standard density functional theory (DFT) fails to stabilise polaron states, 6,7 mainly due to the large self-interaction error of the localised state, several alternatives have been suggested. These are either based on introducing a localised potential in the spirit of DFT+U 4,8,9 or, in the case of hybrid functionals, by incorporating a certain fraction of possibly screened exchange interaction. 3,6,7,10 Although the dependence on external parameters in both approaches can to some extent be removed by invoking Koopman's theorem, both approaches currently suffer from severe drawbacks. Approaches based on local potentials are commonly restricted to polarons localised on atomic sites and thus exclude V K -centres in alkali halides. Hybrid functionals do not suffer from this shortcoming but are associated with a much larger computational cost. Nevertheless, migration barriers in highly symmetric systems such as sodium iodide have been studied through brute force hybrid functional calculations. 7 Beyond these very simple structures, where the local polaron geometry and migration pathways have been known for decades, the computational cost associated with sampling the adiabatic potential energy landscape of extended (periodic) systems has severely impacted theoretical studies.
In particular, these difficulties have been a major obstacle in understanding the polaronic properties of state-of-the-art scintillator materials such as SrI 2 , 11-15 LaBr 3 , 16,17 or Cs 2 LiYCl 6 . 10,18,19 These materials exhibit a much larger chemical, structural and electronic complexity than 'classic' alkali halide materials, possibly leading to unusual effects on polaron formation and migration.
Here we capitalise on a recently devised parameter-free energy functional, the polaron self-interaction correction (pSIC) method, that is suitable for studying polaronic properties and is computationally much more efficient than hybrid XC functionals (up to a factor of 500 for the cases in Sadigh et al. 7 ). As a result, we are able to construct a comprehensive map of the adiabatic potential energy landscape for hole polarons in SrI 2 . We show that the energy landscape features pathways that give rise to an unusually fast polaron migration compared with materials with similar polaron-binding energies. This finding represents a critical step towards our understanding of the excellent scintillator performance of this material. More generally, it demonstrates a powerful approach for the exploration of polaron migration in complex materials.

RESULTS
Stable polaron configurations V K -centres in simple alkali halides form when the removal of an electron from the top of the valence band, which is dominated by halogen-p states, creates an open-shell system. This triggers the formation of a covalent 5pσ bond between neighbouring halogen ions and implies the ejection of an unoccupied antibonding pσ* state into the gap.
Although SrI 2 has a much more complex crystal structure than, e.g., NaI, it is reasonable to expect a similar mechanism also in SrI 2 , as the top of the valence band consists of iodine 5p states. 15 Relaxations starting from randomly distorted structures as well as subsequent molecular dynamics (MD) simulations using pSIC indeed only find polaronic states of this kind. A random search is, however, not suitable to identify all possible dimer pairs in a system as complex as SrI 2 and hence a more systematic approach has to be employed. To this end, we considered all symmetrically distinct first nearest-neighbour iodine pairs between two I sites A and B in the ideal SrI 2 structure.
In the ideal SrI 2 structure (Figure 1a), there are 12 symmetrically distinct I-I pairs that are separated by 3.9-5.0 Å, whereas the second iodine coordination shell is 45.7 Å apart. Initial seed structures for polaron configurations were constructed from these pairs by displacing the two I ions towards each other to obtain a separation of 3 Å with all other atoms fixed. Nine of the possible first nearest-neighbour I-I pairs remain stable after pSIC relaxation with the Perdew-Becke-Ernzerhof (PBE) 20 functional, all resulting in a dimer bond lengths of~3.3 Å (Figure 1b), in accordance with the recent results for V K centres in NaI and LiI. 7 The three remaining I-I pairs are unstable with respect to relaxation into one of the nine stable ones. Subsequent relaxations of the nine stable pSIC-relaxed structures with hybrid functional calculations using the PBE0 parametrization 21 do not appreciably affect the geometry.
The formation energies obtained at the pSIC level range from − 0.34 eV (configuration 1) to − 0.17 eV (configuration 9) ( Figure 2). The most stable dimers (1 and 2) originate from configurations with an initial I-I separation of 4.2-4.4 Å. This is somewhat contrary to the naive expectation that the dimer with the shortest ideal bond length would be energetically favourable. The formation energies obtained using the pSIC method compare very favourably with the results from the computationally more demanding hybrid PBE0 calculations ( Figure 2).
To further elucidate this relation, we considered the (quasiparticle) energy levels for PBE and PBE0 in both the neutral and charged state for the most energetically favourable V K -centre configuration ( Figure 3). In the charge neutral case for both PBE and PBE0, five occupied levels, all related to 5p states of the dimerised iodine pair, are ejected into the bandgap, with the topmost occupied level being the antibonding 5pσ* state. PBE and PBE0 closely agree with respect to the positions of the ejected levels relative to the valence bandedge, and also the overlap for the five ejected bands is very close to unity. The charged state, however, highlights the failure of PBE to describe polaronic states. Although the wave function overlap of the now unoccupied antibonding 5pσ* state is still close to one, the energetic position relative to the valence bandedge differs significantly between PBE and PBE0. Similarly, the bonding 5pσ    Figure 4c) with a barrier height of 0.12 eV. Although already this result suggests that polaron diffusion in SrI 2 is much faster than in, e.g., sodium iodide where the lowest migration barrier is 0.19 eV, 7 the availability of additional pathways further separates these materials. Specifically, in thermal equilibrium at 300 K, the population of type 3 polarons as given by the Boltzmann factor is~5%. Polarons of this type are accessible from configuration 1 via 2 with an effective barrier of 0.14 eV (Figure 4b). They are furthermore connected with each other via a path parallel to the c axis with a barrier of only 0.03 eV (Figure 4d), which enables practically barrier-free diffusion at room temperature. Interestingly, the barrier for migration via the 3 → 3′ network is smaller (0.03 eV) than for the 'recombination' reaction 3 → 2 (0.08 eV). These adiabatic potential energy landscape features imply that already at room temperature polarons in SrI 2 can diffuse extremely fast with one-dimensional characteristics. This situation is reminiscent of the kick-out mechanism leading to fast diffusion of gold in silicon. Here the solubility of Au on an interstitial site is lower than on a regular lattice site but with the reversed relation between their mobilities. 22 Although the foregoing considerations apply under equilibrium conditions, in the following section, we will demonstrate that under irradiation conditions one can expect an even larger polaron mobility due to an effective inversion of the equilibrium population.   It is apparent that self-trapping occurs almost immediately, on the sub-picosecond timescale (see Figure 1 in the Supplementary Material). The fluctuations in the nearest-neighbour iodine-dimer distances are thus sufficiently large to readily relax into a polaronic configuration in the presence of a hole charge and the selftrapping process is effectively barrier free.
By averaging over many MD trajectories, it is possible to determine the distribution of 'as-born' (nascent) carriers over different polaron types (Figure 5a). Compared with the equilibrium distribution, which for simplicity can be approximated by assuming ΔG f ≈ ΔE f , it is striking that the probability for polarons of types 1 (59% in equilibrium versus 14% nascent) and 3 (5% vs 47%) is practically inverted. At room temperature, polarons of type 3 are gradually converted to type 2 over the timescale of 5-15 ps, whereas the number of type 1 polarons remains low (Figure 5b). This is compatible with the lack of a direct pathway between configurations 3 and 1 (Figure 4b).
The population inversion has important consequences for polaron migration. The barrier for 3 → 3′ jumps is only 30 meV (Figure 4b) and thus significantly smaller than for 1 → 1′ events (120 meV). (Note that there is no direct pathway between polarons of type 2). Under equilibrium conditions, type 1 polarons outnumber polarons of type 3, which outweighs the lower barriers available to the latter. The population inversion in the nascent distribution, however, implies that migration via the 3 → 3′ pathway is strongly enhanced. This picture is confirmed by explicit MD simulations, which show extremely rapid and practically athermal diffusion of type 1 polarons ( Figure 6) with tens of hopping events occurring over just 15 ps at room temperature, whereas the polarons of type 1 remain stationary under the same conditions.

DISCUSSION
On the basis of the results presented in the previous section, we will in the following argue that both the inversion of the equilibrium distribution and the rapid migration of type    To this end, we first briefly review the scintillation process.
The incident radiation initially produces fast electrons, either directly by fundamental light-matter interactions or indirectly via electron-electron scattering, that traverse the crystal and deposit their energy by exciting electron-hole pairs and plasmons. [23][24][25] This process occurs on the femtosecond timescale and gives rise to a distribution of low-energy carriers in a nanometre-sized region surrounding the electrons tracks. The excitation density increases with decreasing kinetic energy of the fast electron. 26,27 Nevertheless, at this point, the total number of low-energy electron-hole pairs is approximately proportional to the initial photon energy. In an ideal scintillator without losses, all pairs would transfer their energy to the activator ions, resulting in a photon count that is proportional to the incident energy. In real materials, the low-energy excitations and carriers are, however, subject to annihilation via non-radiative quenching mechanisms such as Auger recombination of free carriers or exciton-exciton annihilation. 28 As the rates of these processes exhibit a non-linear dependence on electron-hole or exciton density, 28 the nonuniformity of the excitation density along the track will give rise to a total light yield that is no longer proportional to the energy of the incident photon energy. It is also apparent that the nonradiative recombination rates are maximal during the initial stage, when the excitation density is largest. In particular, it has been conjectured 29 that the first 10 picoseconds are the most decisive for the proportionality of the light-yield response. Thus, reduction of the excitation density would naturally lead to reduced quenching and increased proportionality. This can be accomplished by means of temporary carrier capture on shallow traps as discussed by Williams et al. 29,30 in the context of Tl-doped CsI. The same effect was obtained by co-doping Ce and Sr in lanthanum bromide. [31][32][33] Extending this idea to deep traps, such as self-trapped polarons, has so far not been considered useful. In common halides materials, such as NaI and CsI, self-trapped (hole) polarons feature both large binding and migration energies, a coupling that is expected based on the structure of the potential energy surface. Polarons therefore exhibit both long detrapping times and slow migration. However, as demonstrated in this work, self-trapped polarons in SrI 2 (i) are characterised by trapping times below 1 ps, (ii) have large binding energies comparable to other halides, (iii) are predominantly 'born' as type 3 polarons and (iv) in the latter configuration exhibit essentially barrier-free migration at room temperature.
The existence of multiple low-energy configurations (which is possible due to the low symmetry of the system) in combination with a comparably soft energy landscape (compare the phonon dispersion) gives rise to a nascent polaron population that strongly favours highly mobile polaron configurations (3). We speculate that these polarons can quickly migrate away from the initial track region while remaining localised (limiting Auger recombination), causing the track to broaden more quickly. This in turn reduces the density of electron-hole pairs during the crucial first picoseconds, and thus lowers the rates of the non-linear non-radiative recombination processes that constitute the root cause of the non-linear response. This mechanism provides a rationale for the superior proportionality of SrI 2 compared with materials such as NaI and CsI.
To contrast the polaron dynamics in SrI 2 and NaI, a material with rather poor energy resolution, we have also performed MD simulations of polaron formation and transport in this material at 300 K using 216-atom supercells. These calculations confirm that even in NaI, self-trapping occurs at sub-picosecond timescale at room temperature. However, we never observe a hopping event in any of our NaI simulations, regardless of initial conditions. It is thus evident that polaron transport is much slower in NaI as compared with SrI 2 .
The present work provides the basis for future large-scale simulations based on, e.g., rate-equation models 28,30 or kinetic Monte Carlo methods, 34 which enable quantitative simulations of light yield based on a set of material-specific input parameters. In this context, we close by commenting on the importance of using suitable methods to obtain the relevant rates: assuming identical pre-factors, one would expect the conversion between 2 and 1 (ΔE m = 0.10 eV) to take approximately four to five times longer exp½ΔE 2-1 m =ΔE 3-2 m À Á than from 3 to 2 (ΔE m = 0.08 eV). To match the MD data (Figure 5b) one, however, requires a ratio that is~20-40 times larger. This deviation emphasises the role of vibrational degrees of freedom, which affect the attempt frequencies and even more so the formation free energies. It thereby highlights the importance of direct MD simulations, as opposed to, e.g., solely harmonic transition state theory, for understanding polaron dynamics.

MATERIALS AND METHODS
The failure of DFT to stabilise polarons in wide gap insulators is related to the self-interaction present in common semi-local exchange-correlation functionals, 6 which manifests itself in the convexity of the total energy as a function of the fractional hole charge. 4,9 The recently developed pSIC method 7 builds on this insight and provides a parameter-free approach for studying polaronic properties that is both accurate and computationally efficient. The pSIC total energy functional for a system with an added electron/hole is given by where E DFT [0; R] signifies the total energy of the ionic configuration R, in the reference electronic charge state, which for applications in this report is chosen as the neutral charge state. Furthermore, μ 7 DFT is the electron chemical potential, i.e., the right/left derivative of E DFT with respect to electron number at the reference electronic configuration. The main approximation in Equation (1) stems from the omission of the largely R-independent many-body scissors corrections to the DFT band energies at the valence band maximum and the conduction band minimum. The pSIC predicted migration barrier for hole polarons in NaI was in good agreement with experiment, 7 encouraging the present comprehensive study of SrI 2 .
SrI 2 with its space-group Pbca (no. 61 in the International Tables of Crystallography 35 ) contains 24 atoms in the primitive unit cell with one Sr and two I sub-lattices, each occupying an 8c Wyckoff site with three degrees of freedom (Figure 1a). We used experimental values for the cell metric: a = 15.22 Å, b = 8.22 Å and c = 7.90 Å. 36 All internal coordinates are fully relaxed in our DFT calculations using the PBE parametrization of the generalised gradient approximation, projector augmented-wave potentials, 37 a cutoff energy of 229 eV and no symmetry constraints as implemented in the Vienna ab initio simulation package. 38 Monoclinic (a, 0, 0) × (0, b, − c) × (0, b, 2c) supercells containing 72 atoms were constructed to study polaron formation and migration together with a 2 × 3 × 2 k-point sampling. For the larger cells used in the size-convergence study (up to 960 atoms), only the Γ-point was included (Supplementary Information). Atomic positions were optimised until the residual forces were o0.01 meV/Å.
In the pSIC method, polaron formation energies are calculated according to where R ideal and R pol denote the ideal and a polaron configuration, respectively. As this method relies on electron chemical potentials computed in the neutral reference state, it is not subject to imagecharge corrections. The pSIC results were benchmarked against PBE0 hybrid calculations, from which V K -centre formation energies are obtained as energy differences between charged and neutral supercells 4 Fast non-equilibrium transport of nascent polarons in SrI 2 F Zhou et al according to Image-charge corrections are applied to hybrid calculations with the modified Makov-Payne correction of Lany and Zunger. 39,40 To this end, we computed both the ionic and electronic contributions to the dielectric constant. In the latter case, we approximately included local field effects for the hybrid functionals by adding the difference between DFT dielectric constants with and without local field effects. In this manner, we obtained the diagonal dielectric tensor ϵ = diag (6.42, 6.13, 7.09) resulting in a monopole-monopole correction of 0.207 eV.
To study the polaron migration between dimers A-B → A-C, a systematic exploration of symmetrically distinct iodine triplets A-B-C was performed, where I A -I B and I A -I C are stable dimers. Here two pairs (triplets) are considered equivalent if there exists a space-group operation that maps A-B to A′-B′ (A-B-C to A′-B′-C′). The search for polaron migration paths was performed using the Compressive Sensing Lattice Dynamics code. 41 The corresponding minimum energy path for each A-B → A-C combination was optimised with the nudged elastic band method 42,43 using three intermediate images.
MD simulations were performed at 300 K using (a, 0, 0) × (0, 2b, 0) × (0, 0, 2c) (96 atom) supercells and an integration time step of 4 fs. The temperature of the system was controlled via a Nosé-Hoover thermostat with the Nosé mass chosen such that the temperature fluctuates with a period of~0.1 ps. We generated 36 separate MD trajectories of the system with a hole charge, the initial ionic coordinates and velocities for which were taken from an equilibration run carried out for the neutral system. The pSIC functional was used to calculate the ionic forces in the presence of the charge excitation. The self-trapping process deposits a large amount of energy in a small volume element equivalent to the polaron-binding energy. For the polarons in SrI 2 this amounts to~0.36 eV in the 96-atom supercell that has been used for MD studies. This heat has to be dissipated before equilibrium can be attained. The Nosé-Hoover thermostat used here equilibrates the temperature of the system within 0.1 ps, which is faster than in the real system. As a result, our simulations tend to underestimate the extent of the non-equilibrium character observed in the MD studies reported in this section. MD duration was limited to t20 ps to focus on polaron dynamics most relevant to the scintillator response.