Bound polaron formation in lithium niobate from ab initio molecular dynamics

Polarons influence decisively the performance of lithium niobate for optical applications. In this work, the formation of (defect) bound polarons in lithium niobate is studied by ab initio molecular dynamics. The calculations show a broad scatter of polaron formation times. Rising temperature increases the share of trajectories with long formation times, which leads to an overall increase of the average formation time with temperature. However, even at elevated temperatures, the average formation time does not exceed the value of 100 femtoseconds, i.e., a value close to the time measured for free, i.e., self-trapped polarons. Analyzing individual trajectories, it is found that the time required for the structural relaxation of the polarons depends sensitively on the excitation of the lithium niobate high-frequency phonon modes and their phase relation.


Introduction
The electrooptic, photorefractive, and nonlinear optical properties of lithium niobate ( LiNbO 3 , LN) [1,2] are exploited in various applications such as surface acoustic wave devices, laser frequency doubling, nonlinear optics or optical switches for gigahertz frequencies. Many LN properties, e.g., optical absorption, frequency conversion and photoconductivity, are profoundly affected by polarons, i.e., self-trapped charge carriers that are screened by a phonon cloud [3][4][5][6][7][8]. The excess charge giving rise to polaron quasiparticles may be the result of chemical reduction or optical excitation [5,9]. This, together with the dependence of the polaron density on the availability of possible lattice trapping sites-which may be locally modified by doping [3]opens many possibilities to engineer the optical coefficients in space and time [10] and allows for fascinating applications such as nonvolatile holographic storage [11,12]. Therefore, there is a strong interest in the detailed understanding of the atomic structure of the polarons as well as their electronic properties, optical response, formation and transport [3,5,[13][14][15][16][17][18][19][20][21][22][23][24].
Concerning the notation, the present work follows Schirmer et al. [3]: Free, i.e., self-trapped polarons are the simplest polaron species in LN. They result from extra electrons trapped at regular Nb Nb ions of the LN lattice. Bound polarons are formed by single electrons localized at Nb Li point defects compensating the lithium deficit observed in the congruently melting composition of lithium niobate. Recent ab initio calculations [7,8] established detailed structure models for these electron polarons that account for both the measured optical absorption peaks and the electronparamagnetic-resonance data.
Generally, the relaxation of hot conduction band electrons into polarons is assumed to occur in three steps [24]: (1) Upon scattering processes the electron reaches the bottom of the conduction band, (2) the electron localizes onto a locally distorted site, e.g. at a regular Nb Nb site (free polaron) [3], or on a Nb Li antisite (bound polaron) [3], and (3) the structural relaxation induced by the trapped extra electron gives rise to a level within the band gap, which can be detected, e.g., by optical absorption.
In case of lead halide perovskites, the polaron-induced strain field dynamics has recently been measured [25]. The experimental studies for lithium niobate indicate that polaron formation takes place on a femtosecond time scale. However, the microscopic processes involved in the polaron formation are not really understood. Williams and co-workers [24] observed a strong temperature dependence of the bound polaron formation time and report values of 230 and 110 fs for measurements at 20 K and at room temperature, respectively. Similarly, times of less than 100 fs and 70-110 fs are reported in Refs. [21,22] and Ref. [20], respectively, for the thermalization of photoexcited carriers into polarons. Beyer et al. [23] conclude from their measurements that bound polarons form in less than 400 fs after optical excitation. To explain the ultrafast generation of bound polarons they suggest that non-thermalized electrons in the conduction band already form free polarons that transform into bound polarons upon interaction with the breathing mode of the oxygen octahedra in lithium niobate. On the theoretical side, Monte-Carlo methods were applied to study the polaron dynamics in lithium niobate [14,18]. In these studies, a pronounced influence of the sample temperature on the polaron diffusion was found.
The present study aims at the microscopic understanding of the polaron formation process in lithium niobate. In particular, we focus on the trapping process of bound polarons, i.e., the localization of conduction band bottom electrons and the structural and vibrational relaxation of the lattice surrounding the Nb Li antisite defect in dependence on the crystal temperature. To that end, ab initio molecular dynamics is employed.

Methodology
The molecular dynamics (MD) calculations presented in the following are based on the Quantum ESPRESSO [26,27] implementation of density functional theory (DFT). The electron-ion interaction is modeled with norm-conserving pseudopotentials that treat the O 2s and 2p, Li 2s, and Nb 4s, 4p, 4d, and 5s states as valence electrons. The PBEsol functional [28] within the generalized-gradient approximation (GGA) is used as starting point to describe the electron exchange and correlation effects. The electron orbitals are expanded into plane waves using an energy cutoff of 85 Ry. The Brillouin zone is sampled with a 2 ×2× 2 k-point mesh.
Congruent, i.e., Li deficient lithium niobate is modelled using a 2 ×2× 2 supercell containing 80 atoms including a single Nb Li antisite atom. The supercell used here overestimates the defect density typical for congruent LN by a factor of two. This needs to be considered when comparing the simulations to experiment. Also, the supercell approach restricts the representation of long-wavelength phonons. However, as will be shown below, in particular defect localized vibrational modes are instrumental for the polaron relaxation.
To correctly account for the single excess electron of the bound polaron, spin-polarized calculations are performed. The charge of the excess electron and Nb defect is compensated by a homogeneous background charge. For a proper description of the strongly localized polaronic excess electron, we need to go beyond (semi-)local DFT. In Ref. [7] it has been shown, that PBEsol+U with the self-consistently determined U value of 4.7 eV for the niobium 4d electrons reproduces the measured EPR parameters, which sensitively probe the spatial distribution of the unpaired excess electron of the polaron ground state. In the present work, focussing on MD calculations that involve structures far away from the ground state, the Hubbard U has to be chosen such that it is large enough to enable the electron self-trapping, but still sufficiently low that no artificial level-crossings are induced during the simulations. These conditions are fulfilled for U d = 2.2 eV, the value chosen here. As will be shown in Section 3, the variation of U leads to an error bar smaller than 10 fs in the calculated polaron formation times.
Molecular dynamics using a Berendsen thermostat [29] applied to systems without excess charge is used to generate equilibrated configurations, i.e., initial atomic positions and velocities, that were subsequently employed as starting points for trajectories that model the polaron dynamics at temperatures of 20, 100, 200, 300, 600 and 1200 K. 40 starting configurations were generated for each temperature. The polaron MD calculations were performed in a microcanonical ensemble with an excess charge of one unpaired electron placed in the conduction band bottom.
The mean distance between the excess electron and the Nb Li antisite atom given by r m = ∫ | − d |n( )d 3 r is an obvious choice to quantify the polaron localization. The time evolution of this quantity during a specific MD run is exemplarily shown in Fig. 1. It can be seen that the electron remains delocalized for about 70 fs, then localizes rapidly within about 15 fs, and remains localized, albeit slightly oscillating, for the remaining simulation time of 370 fs. While the time evolution of the excess electron density is an obvious and direct way to describe the electron localization, it is numerically expensive and storage-intensive, if performed for hundreds of trajectories.
The local magnetic moment local at the Nb Li antisite atom provides an alternative, and numerically far more convenient measure to quantify the electron localization. The spatial distribution of the unpaired electron gives rise to magnetic moments at each atomic site that are directly accessible in spin-polarized DFT. A local magnetic moment of 0.34 B at the Nb Li antisite atom is found to characterize the polaron ground state calculated within PBEsol+U with U=2.2 eV.
The time evolution of the Nb Li local magnetic moment and the average Nb Li -O bond length is shown in Fig. 2.
The data refer to the same trajectory as shown in Fig. 1. The magnetic moment is rather small within the first 70 fs. Then it rises to a maximum at around 120 fs. This parallels the electron localization shown in Fig. 1, as well as the increase of the average bond length between the Nb Li antisite atom and the surrounding O atoms. At the same time, a flat polaron band forms at the conduction band minimum, see bottom panels in Fig. 2. The formation of this rather dispersionless state is characteristic for the electron localization, and could indeed serve as well as criterion for electron localization. When the magnetic moment reaches a local maximum after about 120 fs simulation time, the dispersionless polaron band is lowered into the band gap, i.e., a localized electron state has formed. Simultaneously, the O cage around the Nb Li atom has expanded to nearly maximum size. Subsequently, for the remaining simulation time, the oxygen atoms oscillate around their equilibrium positions for the polaron ground state. These atomic oscillations are accompanied by oscillations in the Nb Li local magnetic moment. Here expanded oxygen octahedra correspond to high magnetic moments and vice versa.
It is clear from Figs. 1 and 2 that the Nb Li local magnetic moment may be used as criterion for the excess charge localization. By comparing the local magnetic moment with explicit electronic densities for several MDs at all temperatures we are able to define a threshold for the localization of the excess electron. To differentiate between localized polarons and random charge fluctuations, we define the excess electron as localized, if the Nb Li local magnetic moment is above 0.18 B for at least 85 % of a 15 fs time interval. The same strategy is used to determine the polaron lifetime: Once the local magnetic moment is below 0.18 B for at least 85 % of a 15 fs period, the polaron is considered to be quenched.
Apart from the electron localization time, also the time subsequently required for the relaxation of the lattice in response to the localized electron is of interest. The lattice relaxation gives rise to an energy level within the band gap, which can be detected, e.g., by optical absorption spectroscopy. To relate our findings to the experimental situation, it is thus helpful to consider the dielectric function. Its imaginary part, calculated for MD snapshot configurations as well as for the ground-state polaron structure, is shown in Fig. 3. The calculations are performed within the independent particle approximation (IPA) employing the Yambo [30] package. The calculations for the static geometries are performed using a Hubbard U values of 5.2 eV and 4.7 eV for the antisite Nb and the other Nb atoms, respectively. These values are known to reproduce the measured polaron EPR data [7].
Generally, the polaron absorption peak shifts to lower energies with higher polaron band position in the LN fundamental gap. The polaron band position relative to the bulk conduction band minimum can thus be used to characterize the completion of the polaron formation. The polaron optical absorption peak corresponding to snapshot e) in Fig. 3 agrees well with the corresponding absorption of the ground-state polaron calculated to be at around 1.8 eV, which is close to the experimental values 1.7-1.9 eV, see Ref. [7]. The corresponding polaron band position of 0.78 eV below the conduction band minimum also Obviously, the average Nb Li -O bond length, the Nb Li local magnetic moment, and the Kohn-Sham eigenvalues of the polaron state are three strongly correlated measures to characterize the polaron formation, see also Fig. 3. Since the local magnetic moment is directly related to the electron localization and the position of the polaron band is directly related to the measured optical absorption, we use these two criteria to quantify the electron localization time and the polaron formation time, respectively.

Results and discussion
The calculated electron localization times and polaron formation times in dependence on the simulation temperature are shown in Fig. 4. A clear trend is seen: The higher the temperature, the longer the localization and formation times. The electron localization times are slightly more affected. The strongest temperature dependence, however, is seen in the polaron lifetime, which reduces strongly with rising temperature. Here the data have to be interpreted with care, however: Many polarons, in particular at low temperature, are not quenched during the simulation time of 0.2 ps. The data shown in Fig. 4 thus represent the lower limit of the lifetimes. In any event, the fact that for temperatures above 1000 K the polaron lifetimes become smaller than their formation times indicates an increasing polaron instability at elevated temperatures (see also Fig. 5b)).
The electron localization times are essentially evenly distributed with a peak at short localization times below 20 fs. This peak attenuates with rising temperature, leading to an inversion of the ratio between fast and slow localizations at around room temperature, see Fig. 5a). Obviously, the pronounced lattice deformation at higher temperatures obstructs the electron localization. A MD starting from the relaxed ground-state geometry of congruent LN results in spontaneous electron localization.
With increasing temperature, the difference between polaron formation and electron localization times, i.e., the time required by the lattice to screen the electron-called relaxation time in the following-decreases slightly as shown in Fig. 4. Furthermore, the relaxation time scatters in a broad range, between about 10 and 40 fs.  Only the MD runs that lead to polaron formation are considered in the statistics shown in Fig. 5. Polaron formation occurs essentially always in the temperature range between 100 and 300 K. In this case, there is, on one hand, enough thermal motion to allow for polaron formation. On the other hand, the thermal energy is not sufficiently high to quench the polaron immediately after formation, cf. Fig. 4. At low and very high temperatures, however, there are occasionally trajectories that do not result in polarons. In case of temperatures above 600 K, this is a consequence of the thermal motion. As shown in Fig. 4 and 5, polarons in the temperature range between 600 and 1200 K are comparatively shortlived. At 20 K, in contrast, there is little thermal motion and systems with unfavorable starting configurations may not form a polaron state within the simulation time. Here, however, a word of caution is in order. Our MD calculations neglect zero-point motion, which is known to strongly affect the LN properties at low temperatures [31,32].
A word of caution is also in order concerning the choice of the Hubbard U. As already mentioned in Sect. 2, an U value of 2.2 eV is used in the MD calculations. To determine the influence of U on the calculated localization and relaxation times, we performed additional calculations with U=4.7 eV and without any U. It is found that using U=4.7 eV generally speeds up the electron localization by 3-5 fs. However, this U value frequently leads to numerical instabilities in the dynamics calculations, i.e., charge sloshing. If no Hubbard U is used, the localization is slowed down by 5-8 fs. This agrees with the expectation that larger U values enforce localization. In any event, we expect the error bar induced by the usage of the PBEsol+U methodology to be smaller than 10 fs.
Another error bar arises from the finite number of simulations and a possible bias of the starting configurations. To probe the the statistical significance of our results, we performed 50 additional MDs at 200 K simulation temperature. The simulations started from randomly determined configurations, with the displacements normalized according to the Boltzmann distribution. This procedure leads to electron localization and polaron formation times that are on average 4 and 8 fs longer, respectively, than the values obtained by starting from equilibrated configurations. From this, we conclude that the statistical uncertainty of the present calculations is smaller than 10 fs.
With the localization times determined, it suggests itself to relate these to local structural motifs that favor or hinder electron localization. For this reason, correlations between the local environment of the antisite Nb, i.e., the geometry of the immediate oxygen, lithium, and niobium neighbors, and the localization times where searched for. Unfortunately, no simple relation was found, indicating that the electron localization is a complex phenomenon that involves more than the local environment.
The structural relaxation subsequent to electron localization requires between 10 and 40 fs. To identify the phonon modes that are relevant for the polaron screening, the ion movements during the polaron structural relaxation are projected on the LN bulk phonon modes. As discussed earlier, the oxygen cage around the Nb Li antisite atom expands upon electron localization, see bond length evolution in Fig. 2. Correspondingly, we find strong contributions of the O related A 2 TO 1 , A 1 TO 1 , A 1 TO 4 , and in particular A 2 TO 3 modes (cf. Ref. [33] for notation) to the lattice relaxation.
Interestingly, while the presence of the localized excess charge clearly modifies the structure evolution during the MD run and its expansion into bulk phonon modes, it does not completely alter the structural dynamics of the system. This holds for all temperatures and is exemplarily shown in Fig. 6 for 600 K simulations, where we compare four example MD runs with and without excess electron starting from identical starting structures and initial velocities.
The trajectories shown in Fig. 6 lhs correspond to the MD dynamics with excess charge, the rhs trajectories refer to the corresponding dynamics for the system without polaron charge. The presence of the polaron charge results in overall larger O-Nb Li distances, but otherwise very similar structural dynamics (see in particular the two uppermost MD runs (a-d) in Fig. 6). For example, the projection of the structure dynamics on the highlighted modes A 2 TO 1 and A 2 TO 3 show a very similar line shape. In line with the increase of the bond length, the frequency of the oxygen vibrations with respect to the Nb Li atom is lowered upon charge trapping, from values around 20.5 THz to about 17.5 THz. Apart from this, it can be concluded that the presence of the localized excess charge slightly modifies rather than drastically alters the temperature driven lattice dynamics. This holds in particular for the low-frequency modes and naturally explains the broad distribution of relaxation times. This can be nicely seen from the two bottom MD runs (e-h) in Fig. 6. The electron localization times for these two MD runs is less than 5 fs, allowing to observe the electron-induced modifications of the relaxation from early on. For the trajectory e) the first maximum of the O displacement (and subsequently the relaxation time) occurs at around 50 fs, whereas for the trajectory g) it occurs already at 20 fs. The different values result from the different phase relations of the excess-electron induced relaxation with the A 2 TO 3 mode. They are out of phase for trajectory e), in contrast to trajectory g). In the former case, the (low frequency) A 2 TO 1 mode is instrumental for the polaron relaxation.
If the temperature regime of 100-300 K is considered, the bound polaron formation in congruent LN is calculated here to be completed within 50 fs and 80 fs. This at the lower end of the experimental observations reporting values between 70 fs and 400 fs [20][21][22][23][24]. In additional to numerical uncertainties at least two reasons may contribute to the difference between experiment and theory: On one hand, the present calculations do not account for the electronic relaxation processes that occur before the excess electron reaches the bottom of the conduction band. On the other hand, the comparatively high Nb Li antisite density in the present simulations may play a role: Replacing one Li with one Nb in an 80 atom unit cell results in a Li/Nb ratio of 0.88, smaller than the ratio of 0.94 typically observed for congruent LN. A formation time of 1500 fs was measured for near stoichiometric LN with a Li/Nb ratio of 0.98 [21], suggesting that the bound-polaron formation time increases with decreasing density of Nb Li trapping sites.
Commonly this is attributed to polaron hopping: The excess electron first forms a free polaron, which is then transported via hopping processes to the antisite. The distance to the nearest antisite, and thus the formation time, is determined by the Li/Nb ratio [21]. A direct formation at the antisite happens only if the electron is excited near the antisite [18]. The distance at which the electron is trapped by the antisite is referred to as trapping radius and has been estimated to be of the order of 8.5 Å for congruent LN [34]. Given that the lattice parameter in the present calculations amounts to 11 Å, it is obvious that the present study cannot account for the formation of free polarons and their subsequent trapping as bound polarons. Preliminary test calculations with larger 160 atom unit cells allowing for a Li/Nb ratio of 0.94 indicate indeed a tendency to longer polaron formation times. The bound polaron formation times calculated here should thus be considered as lower estimates for congruent LN.

Conclusions
Molecular dynamics on the basis of density functional theory was used to investigate the formation of bound polarons in congruent lithium niobate in dependence of the temperature. Thereby we focussed on the localization of excess conduction band electrons and their subsequent transformation to a polaron. A considerable scatter of the polaron formation times is found. This can be traced to the fact that oxygen-related phonon modes are instrumental for the polaron formation and that the polaron structural relaxation depends on the phase relation between the respective modes. The electron localization time, and to a smaller extent, the polaron formation time increase with temperature. At room temperature, bound polarons will form within about 60 fs on average. Longer times can be expected to result from electronic scattering in the conduction band as well as hopping processes required to transport the excess charge in the vicinity of a Nb Li antisite defect.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will