Models of Polaron Transport in Inorganic and Hybrid Organic–Inorganic Titanium Oxides

Polarons are a type of localized excess charge in materials and often form in transition metal oxides. The large effective mass and confined nature of polarons make them of fundamental interest for photochemical and electrochemical reactions. The most studied polaronic system is rutile TiO2 where electron addition results in small polaron formation through the reduction of Ti(IV) d0 to Ti(III) d1 centers. Using this model system, we perform a systematic analysis of the potential energy surface based on semiclassical Marcus theory parametrized from the first-principles potential energy landscape. We show that F-doped TiO2 only binds polaron weakly with effective dielectric screening after the second nearest neighbor. To tailor the polaron transport, we compare TiO2 to two metal–organic frameworks (MOFs): MIL-125 and ACM-1. The choice of MOF ligands and connectivity of the TiO6 octahedra largely vary the shape of the diabatic potential energy surface and the polaron mobility. Our models are applicable to other polaronic materials.


■ INTRODUCTION
Localization of excess charge (electrons or holes) to form small polarons is a common feature of dielectric crystals in general and transition metal oxides in particular. 1−3 There is a competition between the potential energy gained through polarization of the crystal from charge localization and the kinetic energy gained by delocalizing a charge within the extended valence or conduction bands. The localization of charge alters the chemical bonding (effective radii of the elements present) and is associated with a local structural distortion, which accompanies the charge as it moves through the crystal (Figure 1). Small polaron formation, with respect to a delocalized band electron, is enhanced by a large band effective mass and a strong low-frequency dielectric response. 4,5 These features are found in TiO 2 , where the conduction band formed of the Ti(IV) d 0 orbitals is combined with large Born effective charges that give rise to a large dielectric constant. 6,7 The behavior of excess electrons and holes in TiO 2 has been the subject of intense study owing to applications in photoelectrochemistry, photovoltaics, and catalysis. 8−12 Although spatially localized, polarons can influence the optical absorption, reactivity, and electronic conductivity, making knowledge of them essential for materials optimization. 13−15 Pristine crystals are diamagnetic due to the closed shell Ti(IV) d 0 configuration. The one-electron reduction to form spinpolarized Ti(III) d 1 centers allows for selective spectroscopic probes including electron paramagnetic resonance (EPR). For example, an EPR study of an undoped single crystal of rutile TiO 2 by Halliburton and colleagues measured polaron trapping below 15 K, with a characteristic activation energy of just 24 meV. 16 The polaron in TiO 2 surfaces can also be directly observed by scanning tunnelling microscopy, making it an ideal test bed for characterization. 17,18 While TiO 2 has been the subject of various polaron-focused studies, theoretical studies on polarons in metal−organic frameworks (MOFs) are rare, often only looking at the stability of polarons in the framework calculated via cluster models. 19−21 Similar structural motifs to TiO 2 can be found in MOFs. For example, MIL-125 is an electrochromic crystal containing rings of corner-sharing octahedra connected by the 1,4-benzenedicarboxylate (bdc) ligand, 22 while ACM-1 features continuous Ti−O chains with 4,4′,4″,4‴-(pyrene-1,3,6,8-tetrayl) tetrabenzoic acid (H4TBAPy) as a pillaring ligand. 23 To the best of our knowledge, polaron transport studies on such periodic porous framework systems have not yet been performed. However, a number of MOFs based on Ti metal nodes with oxygen linkers have been synthesized and are known to host polaronic states upon electron or hole injection 24 or photoexcitation, 25 opening up the possibility to study polaron hopping from Ti to Ti sites via bridging oxygens comparable to those present in TiO 2 .
The present study concerns the behavior of polarons in TiO 6 motifs within crystalline environments. In particular, we perform an in-depth analysis of the potential energy surfaces in the pristine and the F-doped rutile TiO 2 and the prototype TiO 6 containing MOFs MIL-125 and ACM-1. Analyzing the result of hybrid density functional theory (DFT) with semiclassical electron transport theory, we recover the diabatic potential energy surface, which includes key information on the polaron dynamics. We also highlight the possible avenue to designing materials with optimal polaron mobility for a given application. ■ METHODOLOGY Potential Energy Surface for Polaron Hopping. Nudged elastic band (NEB) calculations were performed using DFT with a plane wave basis-set as implemented in VASP, which was modified by the VTST toolkit. 26−29 The projector augmented-wave method was employed with Ti 3d 3 4s 1 and O 2s 2 2p 4 explicitly treated as valence electrons. 30 The unit cell was expanded to 3 × 3 × 6 supercell, and the Brillouin zone was sampled at the Γ point. The plane wave expansion was cut off at 500 eV, and for the exchange-correlation functional, a screened hybrid exchange-correlation functional, HSE06, was employed. 31 For NEB, seven intermediate states were considered between each polaron localization site.
To efficiently calculate MOFs with a large unit cell (and porous structure), NEB calculations were performed using the numerical atomic-centered basis set in FHI-AIMS 32−36 via its integration into the ASE package. 37 Polarons in MIL-125 were calculated in its conventional unit cell of 240 atoms, while ACM-1 was studied in a 2 × 1 × 1 supercell consisting of 156 atoms. Both systems were sampled at the Γ point, and the same HSE06 functional was used. Additionally, both frameworks were studied with the occupation matrix control method (details in the Supporting Information). 38 Polaron Transport. Small polaron transport within a crystal is commonly described in a quasi-static picture by employing Marcus theory. 39 Within this theoretical framework, the polaron hopping from an initial to a final site is described as the transition between two corresponding diabatic states, which can yield a transfer rate following Landau−Zener theory. For the nonadiabatic approximation to hold, the electronic coupling between the ground and excited adiabatic states must be considerably weaker than the reorganization energy associated with the transfer. 40 The overall electron transfer rate can be described as where H AB is the electronic coupling, k B is the Boltzmann constant, λ is the reorganization energy, T is the temperature, and ΔG 0 is the change in Gibbs free energy between the initial and final states. Obtaining λ depends on having an accurate representation of the diabatic potential energy surface of the final state, in the configuration of the initial state. ΔG 0 depends on the minima at the initial and final states, which are identical between diabatic and adiabatic surfaces. Finally, coupling energy H AB can be obtained by comparing the diabatic and adiabatic potential energy surfaces at the saddle point between initial and final states. These quantities are represented in Figure 2. The diabatic activation energy is often a telling indicator of the efficiency of Marcus transport and is related to the reorganization energy and electronic coupling as where e is the elemental charge, R is the distance of hopping, and n is the number of neighbors. We have set n = 1 throughout this work.
If we could obtain a full diabatic and adiabatic representation of the potential energy surface, we would be able to readily calculate transfer rates for polarons. However, the ground state DFT calculations can only obtain one or the other, depending on the degree of localization of the polaron and the choice of DFT parameters. For instance, Deskins and Dupuis 41 tuned the U parameter in GGA+U calculations of rutile and anatase TiO 2 to enforce localized electronic density up to the saddle point within a linear interpolation. This has the advantage of obtaining a diabatic potential energy surface (PES), which would otherwise necessitate excited state calculations to characterize. Unfortunately, it comes at the cost of using an on-site Coulomb correction of U = 10 eV, where effective physical values for bulk transition metal oxides are between 3.1 and 6.4 eV. 42 Deskins and Dupuis employ additional cluster model calculations based on the formalism of Reference 43 to obtain the missing electronic coupling. Furthermore, a direct comparison cannot be drawn with the analogous adiabatic PES obtained with a lower choice of U, due to the mismatch between levels of theory.
We chose to employ a hybrid exchange-correlation functional, obtaining localized polaron geometries, to produce a 1-D adiabatic PES via either linear interpolation or nudged elastic band. The diabatic potential energy surfaces are approximated as parabolas by sampling the energy values on the corresponding side of the saddle point. The values are then assigned weights according to a Boltzmann distribution, mirroring the increased anharmonicity away from the minimum, and the parabola is fitted by linear regression (details in the Supporting Information). The fitting procedure has been imple- mented in CARRIERCAPTURE, which also calculates Marcus transfer rates and charge carrier mobilities. 44 ■ RESULTS AND DISCUSSION Polaron in a Pristine Crystal. The closest Ti−Ti distances in rutile TiO 2 are in the [001] direction of the crystal. We, therefore, focus on this most efficient electron diffusion channel. To generate an electron polaron, in addition to simply adding an extra electron, we also considered explicit doping by substituting one neighboring oxygen with fluorine (i.e., F O ), as depicted in Figure 3, which we will discuss in the next section.
The geometry of the Ti center and O nearest neighbors is depicted in Figure S1, which agrees with a previous computational study 45 and within 0.2% error from an experimental value. 46 We observe that a small stretch from the bulk geometry of the order of 0.01 Å in all six Ti−O bonds results in the localization of the electron polaron, with a magnetization of 0.75 μ B . These stretches agreed well with previous studies. 41,45 To investigate the polaron hopping rate, we probe its PES by two methods: linear interpolation and NEB. The two sets of data are presented in Figure 4. The energy barrier for the interpolated and the NEB pathways are 80 and 58 meV, and the reorganization energies are 0.33 and 0.23 eV, respectively. The lowering of the reorganization energy from the strained to the relaxed saddle point has a significant impact on the transport properties of the fitted Marcus model, as is reported in Table 1. The hopping rate increases from 9.64 × 10 11 s −1 to 1.59 × 10 12 s −1 . The value obtained from NEB is closest to reality, where the polaron will adopt the least energetic path. Indeed, the diabatic activation energy for the interpolation is 80 meV, compared to 58 meV for NEB and 24 meV in EPR experiments. 16 Experimentally, polaron mobilities have been reported in the range of 0.01−10.0 (cm 2 V −1 s −1 ). 47−49 Although the calculated value by NEB and linear interpolation falls into this range, the majority of mobility measurements are calculated from Hall measurements, where the effect from other carriers can be difficult to separate, so direct comparison requires care. Depending on the condition of the sample and experimental setup, treatment of surfaces, domain, defect scattering, phonon scattering, and finite temperature effects may be required. 47,49,50 Furthermore, at higher temperatures, we expect a sizable contribution from a delocalized state that is only 0.15 eV above the ground state polaron structure, which would require a separate treatment. 51 Deskins and Dupuis calculated a hopping rate of 7.65 × 10 11 s −1 with an activation energy of 90 meV, 41 which is close to our interpolated results. However, the increased energy barrier translates to a 1.15 eV reorganization energy, significantly larger than our 0.33 eV. This is compensated by an electronic coupling of 200 meV, an order of magnitude larger than those found by interpolation and NEB. The large increases in both reorganization energy and electronic coupling cancel out, leading to a modest difference of 1.99 s −1 with our own hopping rate from interpolated geometries. This demonstrates how small changes in the activation energy can be amplified in the hopping rate through the exponential dependence.
Polaron in the F-Substituted Crystal. Next, we calculated the potential energy surface polaron hopping in three Ti site neighboring donors in F-doped rutile TiO 2 ( Figure 5 and Table 1). Here an effective dipole is created between the positively charged donor (F O + ) and the negatively charged polaron (e − ). As such, hopping to and away from the F    Figure 3. b Pristine. c Doped.
Chemistry of Materials pubs.acs.org/cm Article site is no longer equivalent. The reorganization energy of 0.23 eV for the Ti2−Ti3 (escape) hop was the same as in the pristine system, indicating the negligible effect of F beyond this point due to effective dielectric screening by the crystal host. The polaron mobility was noticeably higher for Ti1−Ti2 (0.098 cm 2 V −1 s −1 ) and Ti3−Ti2 hop (0.076 cm 2 V −1 s −1 ). The low activation energy for the former and large coupling for the latter is the origin of the pronounced mobility. Despite the electropositivity of the fluorine donor, the polaron is destabilized at its direct proximity, Ti1, by 20 meV, and the global minimum is located at the second nearest neighboring site, Ti2. This is in contrast to the case of a hydrogen donor where polarons are stable at the closest Tisite. 52 The difference originates from the large size of the fluorine dopant. The ionic radius of F − is slightly larger than that of O 2− , so F − creates a strain field that elongates the bonds around itself with the cost of reducing the volume of the neighboring titanium octahedra ( Figure S1). As the polaron prefers to elongate the surrounding bonds in TiO 2 , the strain field due to the F site competes with the stability of the polaron, and the polaron is pushed away from the F site. Yet, the interaction is not totally repulsive, and Ti2 is preferred over Ti3 due to the Coulombic interaction between the polaron and the F − . If we assume a fully classical picture, ionized fluorine dopant is in the −1 charge state, which is +1 relative to the −2 for oxygen in TiO 2 , and a polaron is −1, which therefore causes an attractive interaction. The actual difference was only 10 meV, so we expect reduced effective charge and a sizable screening effect from the surrounding polarizable crystal. The stabilization due to this effect is comparable to thermal energy at room temperature, so we expect that F only weakly binds the polaron and effectively ionizes at room temperature.
Polarons in Metal−Organic Frameworks. Both MIL-125 and ACM-1 show local chemical environments similar to TiO 2 , where metal sites are bridged by oxygen ( Figure 6). In the former, two distinct hopping processes can be identified, where one is via two bridging oxygens and the other is via a single bridging oxygen atom. Due to the differing Ti site distances, we call them "short" (2.67 Å) and "long" hop (3.67 Å), respectively. ACM-1 on the other hand only has one likely hopping process along a continuous chain of Ti sites bridged by single oxygen atoms (3.43 Å).
Relevant parameters such as hopping length, activation barrier, coupling, and reorganization energy are given in Table  2 for both linear interpolation and NEB. The potential energy surface for short hopping in MIL-125 is plotted in Figure 7. It is clear that linear interpolation fails to describe the hopping of a polaron in the considered MOFs, as it significantly overestimates all parameters compared to NEB. This suggests the non-negligible contribution of a nonlinear relaxation in the linkers/ligands. The results obtained via NEB however are reasonable and can be used for comparison with TiO 2 . Comparing these in Tables 2 and S1 between MIL-125 and TiO 2 , the short hop of MIL-125 had a comparable activation barrier (62 meV) but a smaller coupling of 8 meV. The long hop had comparable coupling energy (21 meV) but a higher activation barrier (91 meV). Both cases resulted in lower mobility. On the other hand for ACM-1, the activation barrier was actually smaller than TiO 2 ; however, the resulting mobility was smaller (0.002 cm 2 V −1 s −1 ) due to weak coupling of 3 meV. These results suggest the importance of having both a small activation barrier and a large coupling.
Based on chemical intuition that the TiO 6 motif is similar across all systems, it is tempting to explain the trend in the activation barrier in terms of the hopping distance. Comparing the two hopping paths in MIL-125, it seems that the hopping barrier follows this trend. However, the activation barrier of the short hop in MIL-125 is higher than that of TiO 2 while having a smaller hopping length, and ACM-1 has the smallest barrier overall despite its relatively large hopping length (3.43 Å). All this implies that the polaron itself is not the same across the systems and that the differences in polarization and localization must be taken into account to explain the resultant hopping behavior.
A stark difference is discerned when we plot the polaron density (spin density) at the PES minima (Figure 8). In TiO 2 , the antinode of the d-orbital wave function is pointing toward  the neighboring Ti site, whereas in MOFs, the d-orbital is extending in the plane perpendicular to the direction of the short hopping. This points toward a different occupation of dorbital levels arising from a change in crystal field splitting between the inorganic and hybrid materials. The symmetry mismatch is likely to have reduced the wave function overlap between the initial and the final site, leading to a decrease in electronic coupling. This is corroborated by our results, which show smaller couplings for the short hop in MIL-125 and in ACM-1, despite activation barriers comparable to those of TiO 2 . As a consequence, hopping rates in the studied MOFs are smaller, of the order of 10 11 s −1 (MIL-125) and 10 9 s −1 (ACM-1) (Table S1). It should be noted however that regardless of the largely d-orbital nature of the polaron at the PES minima, we were unable to reproduce the polaron behavior along the reaction pathway solely by modeling the occupation matrix of the d-orbitals, which suggests some contribution from coupling with the neighboring O 2p orbitals ( Figure S3).
Further insight can be gained by studying the extension in space of a single polaron in the differing structures. We calculated the percentage of polaron density on the main Ti site using the Bader charge analysis in TiO 2 , MIL-125, and ACM-1, and it was 43%, 95%, and 89%, respectively. 53 It is also visible in Figure 8 that the polaron in TiO 2 has "tails" on the neighboring Ti sites in the [001] direction, which have also been reported previously. 54 This points toward polarons in MOF building blocks localizing stronger than in inorganic materials, which is unsurprising due to the suppression of longrange polarization in such low-density nanoporous frameworks. From these results, we suggest that the strength of localization and the orientation of the polaron largely affect the activation barrier. The dimensionality of the TiO 6 chain also plays a role as the strength of electron correlation changes with dimensionality. 55 Indeed, we have previously shown that polarons in 2D TiO 6 layers in Sr 3 Ti 2 O 7 localize differently to TiO 2 . 56 We, therefore, expect the isolated ring geometry in MIL-125 and isolated 1D chain geometry in ACM-1 itself to have the effect of strengthening the localization of polarons.
The above results highlight the diverse behavior of polarons accessible by changing the chemical environment of the TiO 6 chain. MOFs are modular materials, and the ubiquitous choice of the ligand allows for various Ti−Ti distances, TiO 6 orientations, and TiO 6 chain topologies. Although beyond the scope of this study, just like in TiO 2 , 57 intrinsic defects may   Chemistry of Materials pubs.acs.org/cm Article obstruct the polaron transport in MOFs. The analysis in Figure  8 is possibly useful to screen MOFs for further materials discovery as the calculation does not require costly NEB calculations. We have used Bader charge analysis to decompose the polaron charges, but other means of volume decomposition such as Voronoi decomposition or simple integration within the Wigner−Seitz radius are likely to reproduce similar results. While a direct control of the dorbital occupation matrix does not reproduce adiabatic barriers accurately, the values for diabatic transition state energies and reorganization energies appear to be promising, suggesting that it might be a useful tool for low-cost screening approaches. If the lattice relaxation away from the linearly interpolated reaction coordinates is small and does not qualitatively change the polaron hopping behavior, the ∼50 times smaller computational cost of linear interpolation can make it an efficient estimation tool. On the other hand, if one wants to quantitatively compare the initial and the final states, our method to extrapolate the diabatic PES from NEB adiabatic barriers is shown to work well for TiO 6 motifs. Finally, we want to review approximations made in this study and provide a possible avenue for future developments. Treatment of the one-dimensional reaction coordinate is a strong approximation, and when other phonon modes have a comparable electron−phonon coupling strength, they may modify the activation barrier, like in the case of "promoting modes" in nonradiative charge carrier transitions. 50 The diabatic potential energy surface was used throughout this work; however, when the coupling energy is large, adiabatic treatment may suit the problem better. Polaron mobility depends on defect concentration. 47 When defect concentration and species are known, we could perform similar analysis to that of F-doped TiO 2 in this study to obtain the polaron mobility of the whole system. If the polaron is further delocalized to a Frolich polaron, mobility is likely to be explained better by models based on response function. 58

■ CONCLUSIONS
In this work, we have analyzed the formation of transport of Ti(III) polarons in different crystal environments. The diabatic polaron hopping behavior was extracted from an adiabatic firstprinciples potential energy surface. Our method does not restrict any exchange-correlation functional and allows the optimal choice to describe the polaron localization in correlated systems. We showed that F-doping in TiO 2 is different from H-doping and weakly traps a polaron in the second-nearest neighbor site.
Through analysis of MOFs we found that the activation energy could vary largely depending on the connectivity of the TiO 6 chain and the choice of ligands/linkers. The dynamics of polarons in MOFs will be further influenced by the long-range framework topologies, as well as structural defects that can be present in high concentrations. Our work does suggest that a wide range of carrier mobility will be accessible in MOFs through appropriate crystal engineering.
Linear interpolation qualitatively reproduces the nudged elastic band potential energy surface for TiO 2 , suggesting it to be a computationally inexpensive method to screen systems with smaller lattice relaxation. We also show that polaron density could be a descriptor to quantitatively understand the difficulty of hopping. Throughout this work, we focused on the electronic polaron in titanium oxides, but our approach is not restricted to these systems. Hole polarons are known to favor even stronger localization than electron counterparts, which suggests that our method is applicable to them. The validity in other host materials could be assessed by comparing the coupling energy with those reported in this work.
Details on the TiO 6 geometry, potential energy surfaces, Boltzmann weighting, and matrix control method (PDF)