Polaron formation following optical absorption is a key process that defines the photophysical properties of many semiconducting transition metal oxides, which comprise an important class of materials with potential optoelectronic and photocatalytic applications. In this work, we use hematite (α-Fe_{2}O_{3}) as a model transition metal oxide semiconductor to demonstrate the feasibility of direct optical population of band edge polaronic states. We employ first-principles electron–phonon computations within the framework of the density functional theory+*U*+*J* method to reveal the presence of these states within a thermal distribution of phonon displacements and model their evolution with temperature. Our computations reproduce the temperature dependence of the optical dielectric function of hematite with remarkable accuracy and indicate that the band edge optical absorption and second-order resonance Raman spectra arise from polaronic optical transitions involving coupling to longitudinal optical phonons with energies greater than 50 meV. Additionally, we find that the resulting polaron comprises an electron localized to two adjacent Fe atoms with distortions that lie primarily along the coordinates of phonons with energies of 31 and 81 meV.

## INTRODUCTION

Semiconducting transition metal oxides are a promising class of materials for the development of next-generation light-harvesting devices due to their exceptional photochemical stability and sustainable large-scale synthesis from earth-abundant precursors, and have thus been the subject of intense research throughout the past two decades.^{1–4} However, the scope of solar energy conversion applications for transition metal oxides has been limited by their inherent tendency to form polarons, quasiparticles comprised of a charge carrier (electron or hole) bound to a proximal distortion of the host lattice.^{5,6} Polarons arise from strong carrier–phonon interactions and are primarily classified by the spatial extent of their characteristic lattice distortion. A large polaron is one whose carrier wavefunction and lattice polarization extend beyond a single lattice constant. Conversely, a small polaron is contained within the volume of approximately one primitive unit of the crystal structure.^{6} The size of a polaron influences its mobility: Large polarons exhibit band transport properties similar to free carriers, whereas the mobility of small polarons is limited to thermally activated carrier-hopping processes. The theoretical description of polarons dates back to 1933, when Landau first proposed their existence;^{7} however, atomically precise descriptions of polaronic states and their formation in transition metal oxides that fully account for extended solid-state structure and can accurately reproduce experimental spectroscopic data are lacking. Such models are crucial to developing a complete understanding of the electronic structure and photophysics of transition metal oxide semiconductors.

Transition metal oxides are particularly challenging to model within the context of band theory due to their unique electronic structure. In conventional crystalline semiconductors (e.g., II–VI, III–V, and group-14 elemental systems), the electronic structure of the bandgap is dictated by hybridization of valence *s* and *p* orbitals of the atomic constituents of the crystal. Significant orbital mixing within the *sp* subspace results in dispersive bands that accommodate efficient covalent charge transport. In contrast, the electronic structure of transition metal oxides is strongly influenced by additional interactions involving the 3*d* orbitals of transition metal ions, particularly those with an open-shell configuration. As a result, the near-gap bands are generally more ionic in character and, thus, weakly dispersive.^{8,9} This phenomenon has a series of profound consequences: (i) Weakly dispersive bands give rise to a high density of strongly localized electronic states within the lattice, (ii) charge carriers (electrons or holes) occupying these bands experience relatively high effective masses and strong on-site carrier–carrier interactions, and (iii) the energy and mobility of these carriers become inherently coupled to the motion of the surrounding nuclei (i.e., phonons).^{9} Altogether, these consequences comprise two of the most pervasive and fundamental challenges to computational modeling of the electronic structure of transition metal oxides: strong carrier–carrier and carrier–phonon interactions.

Localized electronic states in solids, such as polarons, are notoriously difficult to describe within the confines of approximate density functional theory (DFT). This difficulty is due, in part, to the failure of standard exchange and correlation functionals to describe strong Coulombic interactions between localized electrons. Such self-interaction errors typically lead to an over-delocalization of valence electrons in the ground state predicted by standard DFT. Open-shell transition metal oxide semiconductors are particularly susceptible to these errors due to the strongly correlated nature of their metal *d*-orbital subspace. Fortunately, this problem is now routinely mitigated by the widely used DFT+*U* and DFT+*U*+*J* methods, whereby an on-site Coulombic repulsion (Hubbard) parameter *U* and a site exchange (Hund’s) parameter *J* are added to the total energy functional of standard DFT. Given the additive nature of the corrective functionals, implementation of these methods into existing computational codes is simple; however, the task of choosing appropriate values for the *U* and *J* parameters can prove challenging. Many authors choose to tune the parameters until a particular computed observable agrees with experimental data, while others simply opt to use previously reported values successfully applied in modeling similar systems.^{10} Although prevalent in literature, both routes are poorly justified. Both *U* and *J* should be taken as inherent properties of the system itself and of the particular pseudopotentials and functionals employed in the model. A more rigorous approach is to explicitly calculate the parameters *ab initio* by way of density functional perturbation theory (DFPT) in order to ensure the corrected functional precisely accounts for the magnitude of self-interaction errors present in the specific system of interest.^{11–13}

Inclusion of vibrational degrees of freedom presents further challenges to accurate modeling of polaronic states using traditional computational treatments. Thermal lattice vibrations are described by a stochastic superposition of phonons and represent a departure from Bloch periodicity. The resulting disorder causes states of previously well-defined wavevector to become mixed, often leading to a significant renormalization of the electronic structure. Standard DFT is severely limited by the fundamental approximation that atomic nuclei are completely immobile in the crystal lattice (i.e., fixed at their equilibrium positions). As such, the effects of dynamic thermal fluctuations on the nuclear potential are neglected.^{14} Among the most significant of these neglected phenomena is the contribution of phonon-coupled optical transitions to computed optical spectra. These effects are crucial to the accurate modeling of the temperature dependence of optical functions of materials, such as transition metal oxides, that exhibit strong carrier–phonon coupling and polaron formation in ground and/or photoexcited states.^{14,15}

Here, we use hematite (α-Fe_{2}O_{3}) as a model transition metal oxide to develop a method to compute temperature-dependent electronic structures of strongly correlated semiconductors based on a first-principles atomically precise treatment of carrier–phonon coupling. Among transition metal oxides, hematite has emerged as one of the most well-studied materials due to its potential application as a photoanode for water oxidation.^{16–25} However, bulk charge transport in hematite is constrained by small polaron formation. Although some studies have suggested that both electrons^{26–32} and holes^{31–34} form polarons in hematite following excess charge injection or oxygen-vacancy doping, other theoretical investigations indicate that localization of holes injected into hematite in the absence of excess electrons is unfavorable.^{27,35} Recent reports of transient absorption spectroscopy measurements that probe in the near-infrared (NIR)^{36,37} or XUV^{38} suggest that photogenerated electrons and holes can co-localize in hematite to form self-trapped excitons, also known as exciton polarons. Other transient XUV^{39–42} and pump-push photocurrent (PPPC)^{43} spectroscopy studies have shown that small polarons form alongside thermal relaxation of photoexcited electrons. Our group previously demonstrated that small polaron bands can arise in a pristine hematite lattice without the need for prior charge injection, defect doping, or photoexcitation.^{44} The computational methods developed here reveal several highly localized optical transitions that appear near the band edge of hematite upon inclusion of a thermal distribution of atomic displacements along phonon eigenvectors. By varying the temperature that governs the thermal distribution of phonon displacements in the lattice, we can reproduce accurately the thermal difference spectra of hematite we reported previously.^{44} This work unambiguously confirms that these band tail transitions give rise to the formation of electron small polarons through strong coupling with longitudinal optical phonons.

## THEORETICAL BACKGROUND

### The DFT+*U* and DFT+*U*+*J* framework

The DFT+*U* method is characterized by the addition of a parameterized correction to the total energy functional of standard DFT that is intended to mitigate self-interaction errors arising from strongly correlated electronic states. The resulting DFT+*U* functional can be expressed as

where *E*_{DFT} is the standard DFT total energy functional, $nr$ is the electron density, *E*_{U} is the Hubbard correction, and $nmm\u2032I\sigma $ are the orbital occupation matrix elements for each atom *I* within the Hubbard manifold.^{45–47} Here, *σ* represents the electron spin index, ↑ or ↓. In the original formulation proposed by Anisimov and co-workers, *E*_{U} was not rotationally invariant with respect to the localized Hubbard subspace and, as such, depended heavily on the choice of basis for the Hubbard orbital wavefunctions.^{13} Today, the most widely used DFT+*U* method is the rotationally invariant form, first developed by Liechtenstein and co-workers^{48,49} and further simplified by Cococcioni and de Gironcoli.^{11} In this formulation, *E*_{U} is expressed as

where *U*^{I} is the atomic Coulomb repulsion and **n**^{Iσ} is the occupation matrix of the Hubbard atom *I*. In this form, *E*_{U} effectively acts as a penalty function that favors either full or null occupancy of the Hubbard subspace over fractional occupations, thus driving the system toward a state more akin to that predicted by exact DFT. The elements of **n**^{Iσ} take the general form of a projection of the electronic Kohn–Sham valence wavefunctions ($\psi kv\sigma $), denoted by crystal momentum **k** and band index *v*, onto a choice of localized basis functions,

In Eq. (3), $fkv\sigma $ represents the Fermi–Dirac occupation number. Throughout this work, we employ a basis of atomic iron (Fe^{0}) *d* orbitals, $\phi mI$, where the index *I* now specifically denotes iron atomic sites within the lattice. As such, the projection operator takes the simple form given in the following equation:

In this work, we determine the value of *U*^{I} for the iron 3*d* subspace of hematite within the framework of the linear response approach developed by Cococcioni and de Gironcoli.^{11} The aim is to compute, from first principles, a value for the Hubbard parameter that corrects for the unphysical curvature of the total energy induced by non-integer occupation numbers. This goal is accomplished by way of constrained density functional calculations that provide a means of tracking the total energy of the system as a function of the local occupation of the Hubbard atoms. In practice, these constraints are implemented by applying localized perturbations, *α*^{I}, to the Kohn–Sham potential of individual Hubbard sites, *I*, isolated within a supercell. Thus, the response function given in the following equation affords a quantitative measurement of the total energy curvature with respect to changes in the on-site occupancy of the Hubbard manifold:

Elements of the response matrix are calculated as numerical derivatives by evaluating the occupations (*n*^{Iσ} = Tr[**n**^{Iσ}]) subject to a discrete set of applied potentials centered tightly around 0 meV, typically within ±100 meV, starting from the self-consistent potential of a converged standard DFT ground state. The Hubbard parameter for a given site, *I*, is subsequently calculated as the difference in the inverse of two unique response functions, i.e.,

The first (*χ*_{0}) represents the bare, noninteracting response evaluated after a single iteration of the self-consistent field calculation and the second (*χ*) is the self-consistent, interacting response evaluated following convergence. Subtraction of the two response functions effectively removes the total energy curvature associated with the instantaneous rehybridization induced by the applied potential, leaving only the Hubbard *U* of the interacting system.

At its conception, the methodology described above explicitly neglected higher-order terms in the Coulomb interaction (i.e., those described by Hund’s *J*) in favor of a simplified description of the on-site Coulomb repulsion. The resulting Hubbard *U* was taken to be an effective linear combination of both the true Hubbard *U* and Hund’s *J*,

Although Eq. (7) is ubiquitous in literature, its form has yet to be rigorously justified.^{13} Himmetoglu *et al.* later proposed an extension of the linear response method to the calculation of *J*.^{12} Notably, they found that the correct insulating ground state of cubic CuO could be reproduced only with the inclusion of Hund’s *J* as the original DFT+*U* formulation neglected stabilizing magnetic interactions, thus predicting a metallic ground state. These interactions are reintroduced to the functional *E*_{U} according to the following expression:

The spin index *σ*_{min} specifically denotes the minority spin channel. The method for calculating *J*^{I} closely follows that for the calculation of *U*^{I}. Here, the response matrix elements take the form

which defines the response of local magnetization ($mI=nI\u2191\u2212nI\u2193$) to an applied magnetic perturbation of magnitude *β*_{J}. Hund’s *J* is then calculated according to the following equation:

Importantly, the functionality of this technique is limited to the calculation of *J* from a nonmagnetic ground state in order to preserve linearity of the response matrices. This limitation is a direct consequence of the fact that the total energy of a system is not variational with respect to its magnetization as ground state magnetization will typically tend toward saturation.^{13}

We employ this extended linear response method to independently calculate the Hubbard *U* and Hund’s *J* parameters for the Fe 3*d* subspace of hematite. Starting from a fully relaxed standard DFT ground state, local potential shifts (*α*^{I} and *β*^{I}) ranging from −60 to 60 meV are sequentially applied in 20-meV increments to a single Fe center isolated in a 2 × 2 × 2 rhombohedral supercell of hematite. In order to obtain an accurate noninteracting response function, the accuracy threshold of the first iterative Hamiltonian diagonalization of each perturbation was fixed to its final value at convergence of the preceding DFT ground-state calculation. Due to the symmetry-equivalence of the iron centers in hematite, the response function can be computed from perturbing a single iron atom. In order to account for the background response of the lattice, we perform an additional set of perturbations in which equivalent potential shifts are applied to every atom in the supercell *except* the isolated iron. Therefore, each response function takes the form of a 2 × 2 matrix.

To impose self-consistency on the calculated values of *U* and *J*, the response functions are iteratively recalculated from a DFT+*U* and DFT+*J* ground state, respectively. In this scheme, we apply *U*_{in} (*J*_{in}) and recalculate an updated value, *U*_{out} (*J*_{out}). This sequence is repeated until the output value is within 1% of the input value. We note that an iterative recalculation of *J* from a DFT+*J* ground state without the inclusion of the Hubbard *U* parameter may not be generally applicable to all materials. Many systems are driven toward a metallic ground state as the value of the applied *J* parameter increases (i.e., the bandgap narrows with increasing *J*). Consequently, iterating on a DFT+*J* ground state is rendered impossible if the bandgap of the system closes. In such cases, one may choose to first compute a self-consistent *U* and then iteratively compute *J* from a fixed DFT+*U* ground state. In the case of hematite, however, increasing both *U* and *J* widens the bandgap and preserves an insulating ground state.

*Ab initio* electron–phonon coupling via vibrational averages

The earliest electron–phonon coupling computations were based on perturbative methods aimed toward the direct calculation of electron–phonon matrix elements.^{50–52} Here, we use the non-perturbative method of vibrational averages, which involves computing electronic observables within a supercell subjected to a set of atomic displacements that approximate an average thermal configuration, thus circumventing explicit calculation of the interaction matrix. This approach was first developed by Williams^{53} and Lax^{54} and has recently been extended to the computation of temperature-dependent band structures and optical absorption spectra by Zacharias *et al.*^{55–57} and Monserrat *et al.*^{58–60} Notably, perturbative methods of computing electron–phonon matrix elements are limited to include only terms up to a given order. In the vibrational averages approach, this limit is surpassed by virtue of performing full electronic structure calculations on a given configuration. Thus, higher-order terms are implicitly captured in all computed observables.^{61}

Following the adiabatic Born–Oppenheimer approximation, the temperature-dependent expectation value of an electronic observable can be defined as

where $\Phi iu$ represents nuclear states of energy *E*_{i} described by a nuclear configuration expressed in terms of normal mode coordinates, **u** = {*u*_{mq}}. Subscripts *m* and **q** represent the phonon band index and wavevector, respectively. $Z$ is the canonical partition function and *k*_{B} is the Boltzmann constant. The central goal of the method of vibrational averages is to construct a nuclear configuration (or set of configurations), $u\u0304T$, that produces a value of the observable equivalent to the vibrational average at a given temperature, *T*, i.e.,

Within the quadratic approximation of the harmonic oscillator, Eq. (12) can be converted into a set of 3(*N* − 1) equations that determine $u\u0304T$, where *N* is the total number of atomic coordinates,

In Eq. (13), *S*_{mq} takes the form of a sign matrix with elements ±1 and $nB\Omega mq,T$ represents the Bose–Einstein occupation number of a phonon with frequency Ω_{mq} at temperature *T*.^{58–60}

Real-space representations of vibrational displacements are restricted to phonons with real-valued eigenvectors (i.e., those with **q** = 0). Therefore, the sampling of finite-wavevector phonons requires the construction of a supercell commensurate with the desired **q**-grid. Herein, we omit the **q** subscript as all phonons that are considered are effectively “folded” into the center of the first Brillouin zone of the chosen supercell. Equation (13) then yields a set of $23N\u22121$ atomic configurations, where *N* now denotes the number of atoms in the supercell whose nuclei have been displaced from their equilibrium positions according to the following equation:

where *M*_{κ} is the mass of atomic species *κ* and *e*_{κm} are the phonon eigenvectors obtained from diagonalization of the dynamical matrices.^{61–64}

In practice, the methodology described above demands careful consideration of the balance between accuracy and efficiency. Dense sampling of the phonon dispersion curve requires the use of large supercells; however, the configuration space defined by Eq. (14) scales dramatically with system size and is staggeringly vast for even the smallest supercells. In recent years, a number of sophisticated approaches toward overcoming this challenge have been reported. Zacharias *et al.* demonstrated an efficient Monte Carlo integration scheme,^{55} a one-shot approach to generating optimized sign matrices,^{56} and a rigorous algorithm to construct a single temperature-dependent atomic configuration.^{57} Monserrat proposed a method of so-called thermal lines, whereby the temperature dependence of electronic observables is computed along a small set of lines in the configuration space (i.e., by keeping the sign matrices constant as temperature is varied).^{58} In this work, we adopt an approach that combines Monte Carlo integration with the method of thermal lines in order to generate a small subset of configurations that closely approximate the thermal average.

## COMPUTATIONAL METHODS

### Electronic structure computations

First-principles DFT and DFT+*U*+*J* computations were conducted utilizing the pseudopotential plane wave method implemented in the *Quantum ESPRESSO* package.^{65–67} Throughout this work, we employed the exchange–correlation functional proposed by Perdew, Burke, and Ernzerhof and revised for solids (PBEsol)^{68,69} and Optimized Norm-Conserving Vanderbilt (ONCV)^{70,71} pseudopotentials. Converged values of the plane wave cutoff energy and *k*-point grid density were determined for a 10-atom rhombohedral primitive cell of hematite with nuclei clamped at their equilibrium positions. A high plane wave energy cutoff of 1100 eV was chosen to ensure total energy and total interatomic forces were converged to within 10 meV and 0.10 meV/Å per atom, respectively. Self-consistent field (SCF) optimization of the Kohn–Sham wavefunctions was converged on a 4 × 4 × 4 Monkhorst–Pack^{72} *k*-point grid, sampling 28 symmetry-weighted *k*-points in the first Brillouin zone of the primitive cell. For all SCF calculations involving supercells, the *k*-point grid density was reduced commensurate with the dimensions of the supercell. Iterative atomic relaxations were performed utilizing the Broyden–Fletcher–Goldfarb–Shanno (BFGS)^{73–75} algorithm until the total interatomic forces were less than 0.03 meV/Å. Lattice parameters were fixed to an average of experimental values reported for hematite.^{76,77} The antiferromagnetic sublattice of hematite (↑↓↓↑ along the principle axis) was constructed by explicitly treating spin-up and spin-down iron centers as unique atomic species with a fixed antiparallel spin configuration.

### Vibrational structure computations

Full phonon dispersion curves were calculated with density functional perturbation theory (DFPT) as implemented in the *PHonon* code distributed with the *Quantum ESPRESSO* package.^{65–67} Starting from a fully relaxed DFT+*U*+*J* ground state, dynamical matrices were evaluated on a uniform 3 × 3 × 3 grid of *q*-points. Finer sampling was then achieved via Fourier interpolation on a denser **q**-grid.

## RESULTS AND DISCUSSION

### Electronic and vibrational properties of the static DFT+*U*+*J* ground state of hematite

We begin by using the linear response method for determining *U* and *J* to compute electronic and vibrational structures and spectra of the static ground state of hematite, in which all atoms are clamped in their geometrically relaxed positions. Linear response calculations performed within a 2 × 2 × 2 (80-atom) rhombohedral supercell of hematite produced a Hubbard *U* parameter of 3.120 ± 0.006 eV and a Hund’s *J* parameter of 1.535 ± 0.008 eV. Self-consistency was achieved within five iterations, following which both *U*_{out} and *J*_{out} were consistently within 0.1% of their respective preceding values (see Fig. S1 in the supplementary material). We confirmed that these values were converged with respect to the volume of the supercell by performing an identical set of linear response calculations within a 3 × 2 × 2 (120-atom) supercell. The output parameters were within 2.5% of those reported above. Therefore, we conclude that a 2 × 2 × 2 supercell is sufficient for hosting single-atom perturbations and adequately suppresses spurious interactions of the localized potential with its periodic image.

The electronic band structure and optical dielectric spectrum computed from the static DFT+*U*+*J* ground state of hematite are shown in Figs. 1(a) and 1(b). We chose to apply a rigid shift of +0.5 eV to all conduction band eigenvalues in order to bring the computed dielectric function into agreement with the measured spectrum. Herein, this shift is applied to all electronic band diagrams, electronic density of states plots, and computed dielectric functions. In order to fulfill the *f*-sum rule governing total oscillator strength, all computed optical spectra are subsequently renormalized by a factor of $1\u22120.5eV\u210f\omega $.^{55,78} With the applied shift, the computed bandgap is ∼2.4 eV and is characterized by a direct transition midway between the zone-center and the **X**-momentum critical point; however, given the weak dispersion of the conduction band edge, several indirect gaps of similar magnitude exist throughout the Brillouin zone. Consequently, there is little value in definitively classifying the static ground state of hematite as either a direct or indirect semiconductor.

Projection of the DFT+*U*+*J* Kohn–Sham wavefunctions onto an atomic orbital basis [Fig. 1(a), right panel] indicates that the uppermost valence bands (>−4 eV) are predominantly of O 2*p* character with minor Fe 3*d* hybridization. These mildly dispersive bands become maximally hybridized in a narrow region (∼−5 eV) before giving way to the more ionic Fe 3*d* orbitals that comprise the less dispersive bands at the bottom of the valence continuum (<−5 eV). Unoccupied Fe 3*d* orbitals give rise to the minimally dispersive conduction bands and clearly exhibit crystal-field splitting characteristic of high-spin octahedral Fe^{3+} ions. The t_{2g} and e_{g} conduction edges are separated by ∼1 eV, in agreement with reported experimental values.^{79,80} In accordance with ligand-field theory, the higher-energy e_{g} bands exhibit significantly higher O 2*p* hybridization than their more ionic t_{2g} counterparts. Higher-energy (>4 eV) conduction bands are comprised of highly dispersive Fe 4*s* orbitals and contribute negligibly to the density of states, particularly when compared to the sharply peaked density of the Fe 3*d* conduction bands. As a result, the visible and near-UV absorption spectrum of hematite is dominated by optical transitions from O 2*p* valence bands to Fe 3*d* conduction bands.

The single-particle optical dielectric spectrum of the static DFT+*U*+*J* ground state of hematite was obtained by evaluating the momentum matrix elements coupling each pair of valence and conduction band wavefunctions according to the following equation:^{65–67}

where *V* represents the volume of the primitive cell, *N*_{k} is the total number of **k**-points sampled, and $p\u0302$ is the momentum operator. Equation (15) does not account for the finite lifetime of the optically excited state; therefore, we chose to apply a homogeneous spectral Gaussian function with a linewidth of 0.3 eV to the computed spectrum in order to approximate lifetime broadening. Herein, this broadening function is applied to all computed spectra. We have previously reported the optical dielectric spectrum of hematite obtained via Fresnel analysis of transmission and reflectance spectra measured for a series of polycrystalline thin films of hematite.^{44} As shown in Fig. 1(b), all notable features of the measured dielectric spectrum above 2.4 eV are reproduced in the computed spectrum, with only minor differences in intensity and energetic positions. Importantly, the absorption onset near 2 eV is absent from the computed dielectric function. In the next section, we demonstrate that optical transitions in this region are recovered following the inclusion of electron–phonon coupling.

The DFT+*U*+*J* phonon dispersion curve of hematite is shown in Fig. 1(c). As indicated by the projected vibrational density of states, lower-energy phonon modes (<30 meV) correspond primarily to the displacement of Fe atoms, while higher-energy modes (>50 meV) correspond to displacements of O atoms. The majority of the optical phonon branches exhibit minimal dispersion and produce strong peaks in the density of states. Most notably, the two highest-energy branches are nearly dispersionless, giving rise to an anomalously high density of states at ∼81 meV. Additionally, two phonon bandgaps (i.e., regions of near-zero vibrational density) appear in the phonon dispersion, indicated by the horizonal arrows in Fig. 1(c). The first is a narrow gap that follows the sharp drop in the density of states just below 50 meV. A second gap appears just below 57 meV and is effectively much wider, as the region between 57 and 62 meV is crossed only by a single, highly dispersive band near zone center. Gaps in the phonon density of states are known to inhibit the decay of high-energy optical phonons into lower-energy optical and acoustic phonons.^{81–83} Pathways for phonon relaxation are subject to energy and momentum conservation; therefore, allowed routes for the decay of vibrational population above a gap may be severely limited, potentially leading to a buildup of long-lived optical phonons. This phenomenon can profoundly impact the thermalization of hot carriers. If not permitted to decay, high-energy optical phonons emitted during thermal relaxation of an excited carrier can, in turn, prolong the lifetime of the carriers if reabsorbed. As such, we predict that phonon population near these gaps will exhibit strong coupling to photoexcited carriers.

The shaded regions of the dispersion curve indicate branches of predominantly longitudinal character: Red and blue regions correspond to longitudinal optical (LO) and longitudinal acoustic (LA) branches, respectively. Details of this characterization can be found in the supplementary material. Interestingly, the LA branch exhibits numerous avoided crossings with lower-energy (20–30 meV) LO branches, indicating a high degree of LA/LO phonon hybridization away from the zone center. Acoustic/optical hybridization introduces additional carrier–phonon coupling to the acoustic branches and, similar to the occurrence of a phonon bandgap, removes viable pathways for the decay of optical phonons into acoustic phonons.^{84}

A comparison of the previously measured Raman spectrum of a polycrystalline thin film of hematite under band-edge excitation (2.21 eV) with the DFT+*U*+*J* vibrational density of states is shown in Fig. 1(d).^{44} Vertical lines indicate the positions of all computed **Γ**-point Raman-active phonons. The computed highest-energy mode is positioned at ∼81 meV; therefore, we attribute all higher-energy (>90 meV) features of the measured spectrum to second-order, multi-phonon scattering processes.^{44} We compare the first-order region of the Raman spectrum (<90 meV) to the zone-center density of Raman-active LO and transverse optical (TO) modes, obtained by sampling a dense **q**-grid tightly centered around the **Γ**-point ($q\u22642.5\xd7104cm\u22121$). Each of the bands observed in the measured spectrum is accurately reproduced in the computed density spectrum. Notably, the measured bands at 31 and 51 meV correspond to regions of predominately LO-phonon density in the computed **Γ**-point spectrum.

Although first-order Raman scattering is restricted to phonons of near-zero wavevector, second-order scattering processes can access the entire phonon density of states while still conserving the wavevector. In the simplest cases, these processes involve the simultaneous scattering of two phonons of equal and opposite wavevector from the same branch. In Fig. 1(d), the second-order region of the measured Raman spectrum is overlaid with the 2-LO and 2-TO density of states. Both density spectra reproduce the second-order Raman features; however, the 2-LO spectrum exhibits the same intensity distribution as that of the measured spectrum. Additionally, the lowest-energy band of the second-order spectrum at ∼100 meV coincides with twice the computed phonon bandgap at 50 meV, suggesting that two-phonon Raman scattering involves only modes above this gap. Overall, the computed phonon density of states accurately accounts for all of the features observed in the first- and second-order regions of the Raman spectrum. In the following sections, we demonstrate that the resonantly enhanced second-order Raman spectrum corresponds to thermally activated optical transitions into polaronic states near the band edge of hematite.

### Thermal difference spectra of hematite computed via *ab initio* electron–phonon coupling

With the self-consistent values of *U* and *J* and the vibrational structure of hematite in hand, we are now positioned to model the temperature-dependent electronic structure of hematite using the method of thermal lines combined with Monte Carlo averaging. The experimental data used to benchmark these calculations are a series of temperature-dependent thermal difference spectra reported in our previous work.^{44} Here, we define a thermal difference spectrum, Δ*ϵ*_{i}(*ℏω*, *T*), to be the difference between the optical dielectric spectrum measured at a particular temperature and that measured at room temperature, i.e.,

Stochastic Monte Carlo averaging of the computed thermal difference spectrum (TDS) of a 2 × 2 × 2 supercell of hematite at 80 K converged rapidly with respect to the number of thermal lines sampled. As shown in Figs. 2(a) and 2(b), the result obtained by averaging spectra computed across a random distribution of 25 thermal lines is nearly identical to that obtained by sampling only six lines. In fact, we demonstrate that even a *single* thermal line can sufficiently reproduce the Monte Carlo average [Fig. 2(c)]. A 2 × 2 × 2 supercell of hematite permits real-space representation of 237 phonon modes within the primitive Brillouin zone, specifically those at the **Z**, **F**, **L**, **L**_{1}, and **Γ** critical points. In Fig. 2(d), we demonstrate that stochastic sampling of the configuration space of a larger, 3 × 2 × 2, supercell (357 phonon modes) produces a similar average TDS at 80 K. We therefore conclude that a 2 × 2 × 2 supercell is sufficiently large for computing temperature-dependent optical spectra of hematite. Finally, we note that all the computed TDSs shown in Fig. 2 accurately reproduce the previously reported measured spectrum at 80 K^{44} and recover the band edge optical transitions absent from the computed static dielectric spectrum (Fig. 2 insets).

We previously reported TDS of a hematite thin film measured from 30 to 573 K.^{44} Here, we compute temperature-dependent TDS averaged across five thermal lines of a 2 × 2 × 2 supercell of hematite, each closely approximating the Monte Carlo average. Several representative comparisons of the computed and measured spectra are shown in Figs. 3(a)–3(d). At low temperatures (<294 K), the positions and relative intensities of the features in the measured TDS are accurately reproduced in the computed spectra [Figs. 3(a) and 3(b)]. We note that, due to the UV absorption edge of the cryostat windows, these measurements were restricted to photon energies below 3.5 eV. High-temperature (>294 K) spectra were collected without the use of a cryostat and, therefore, extend to 6.4 eV. In Fig. 3(c), we compare the measured and computed TDS at 473 K, illustrating that the computed spectrum reproduces the near-UV features of the measured TDS, albeit with a significant deviation in the band edge intensity.

We characterized the temperature dependence of the computed and measured spectra by integrating the absolute value of Δ*ϵ*_{i} over the interval 1.0–3.5 eV [Fig. 3(d)]. We define the total intensity to be the negative absolute value for spectra collected at temperatures below 294 K and the positive absolute value for spectra collected at temperatures above 294 K. The measured and computed temperature-dependent integrated intensities display remarkable agreement, with minor deviations arising at temperatures above 400 K. Both sets of data were fit to a Bose–Einstein distribution function according to the following equation:

Here, *A* is a constant that corresponds to the value of ∫Δ*ϵ*_{i} at 0 K and Δ*n*_{B} represents the change in the Bose–Einstein occupation number of a phonon of energy *ℏ*Ω relative to its occupation at room temperature (294 K). The temperature dependence of the computed TDS fits the thermal distribution of a 50.3 ± 0.1 meV phonon, while the measured spectra correspond to that of a 49.2 ± 1.6 meV phonon. These values suggest that the temperature-dependent growth of optical transitions near the bandgap of hematite coincides with population of phonons above a 50-meV threshold.

### Correspondence between thermal difference and resonance Raman spectra of hematite

In this section, we use the method of thermal lines to show that including displacements along only the eigenvectors corresponding to the phonon modes that are most strongly coupled to the near band-edge excitation of hematite, as evidenced by resonance Raman excitation profiles, can reproduce all of the features observed in the thermal difference spectra. Our previously reported resonance Raman spectra of a hematite thin film collected with excitation photon energies spanning the visible range are shown in Fig. 4(a).^{44} Relative scattering intensities of the observed bands exhibit a strong dependence on the excitation photon energy. Most notably, the second-order spectrum (>90 meV) is observed only for excitations that are resonant with band-edge transitions (2.0–2.5 eV). In Fig. 4(b), we show the excitation profile of the integrated intensity of each second-order Raman band. The bands at 100, 129, and 136 meV share a similar excitation profile and are strongly enhanced at 2.21 and 2.53 eV. Importantly, each of these profiles closely trace the band edge TDS features. The 163-meV band is similarly enhanced at 2.21 eV; however, it exhibits markedly less enhancement at 2.53 eV. This discrepancy is addressed in the following section. These results strongly suggest that the TDS and second-order Raman spectrum of hematite arise from the same phenomenon.

The majority of the first-order modes (<90 meV) in the observed Raman spectrum share a similar excitation profile (see Fig. S2 in supplementary material), with the exception of the three bands shown in Fig. 4(c): at 31, 51, and 81 meV. Both the 31- and 51-meV bands correspond to regions of high LO-phonon density computed at the zone center [see Fig. 1(d)]. Similar to the second-order modes, the 31- and 51-meV bands exhibit notable enhancement at photon excitation energies of 2.21 and 2.53 eV, respectively. The excitation profile of the 81-meV phonon is nearly identical to that of the 100-, 129-, and 136-meV bands and also closely reproduces the relative intensity of the band edge TDS. We note that this band may be, in part, a two-phonon band arising from a combination of the 31- and 51-meV bands. Complete assignment of all conceivable multi-phonon combinations is beyond the scope of this work. Instead, we draw the general conclusion that the second-order Raman spectrum and band-edge features of the TDS arise from optical transitions coupled to phonon population above the 50-meV phonon bandgap.

We further support this conclusion by performing a series of calculations similar to those used above to generate the computed thermal difference spectra, but with the inclusion of only phonons of energy greater than 50 meV where the spatial phases of these phonons have been synchronized such that their displacements are localized to a single iron center. We find that atomic displacements along the coordinates of the LO phonon modes above the 50-meV gap are sufficient to reproduce the features of the thermal difference spectrum [Fig. 4(d)]. This result, combined with the correspondence of the 2-LO density of states with the intensity distribution of the second-order Raman spectrum [see Fig. 1(d)], indicates that the optical transitions near the bandgap exhibit strong coupling to LO phonons with energies greater than 50 meV.

### Thermally activated polaronic optical transitions in hematite

Longitudinal optical (LO) phonons are known to play a critical role in the formation of polarons due to the strong lattice polarizations induced by their propagation.^{85} In this section, we show that these polarizations result in the formation of a series of polaronic band tails that give rise to the features of the TDS. We connect our model of the temperature-dependent electronic structure of hematite to polaron formation by identifying bands in the electronic dispersion curves that appear at finite temperature and evaluating the localization of wavefunctions associated with these bands.

To assign the transitions that appear in the optical dielectric spectrum at finite temperature, we compute electronic dispersion curves at two temperatures (294 and 573 K) in the configuration space of a 2 × 2 × 2 supercell of hematite along a single thermal line: that which best approximates the Monte Carlo average of the TDS at 80 K [see Fig. 2(c)]. The supercell band structures were unfolded to span the first Brillouin zone of the primitive cell according to the method proposed by Zacharias and Giustino.^{57} The resulting dispersion curves, along with their corresponding densities of states, are shown in Figs. 5(a) and 5(b). The bands are represented by a spectral function, *ρ*, that effectively represents a momentum-resolved density of states. Significant thermal broadening appears in the dispersion curve at finite temperature, giving rise to numerous phonon sidebands and illustrating the profound strength of electron–phonon coupling in hematite, particularly within the Fe 3*d* orbitals. Two minimally dispersive band tails appear just below both the t_{2g} and e_{g} conduction bands, with wavefunctions strongly localized to two neighboring Fe^{3+} ions. Additionally, a weakly dispersive tail appears at the valence band edge, moderately localized to the O^{2−} ions surrounding the two Fe^{3+} ions. The features of the computed TDS arise primarily from optical transitions coupling the uppermost valence bands to the localized conduction band tails.

As seen in the temperature-dependent density of states plotted in Fig. 6(a), the t_{2g} and e_{g} band tails become less dispersive and more localized with increasing temperature, and they exhibit a significantly stronger thermal shift than their corresponding continua. The shift of these band tail states to lower energy with increasing temperature effectively decreases the energy gap between the valence band maximum and the conduction band minimum, thereby accounting for the bandgap shrinkage that previous reports have invoked to explain features observed in the thermal difference spectra.^{86,87} Importantly, these conduction band tails are present at 0 K and therefore dictate the zero-point bandgap renormalization. Thus, the features of the TDS do not indicate the *formation* of new optical transitions with increasing temperature, but rather the *evolution* of transitions that are always present. At the zero point, the t_{2g} band tails are separated by 30 meV. This splitting increases monotonically with temperature, with the lower-energy band shifting at a faster rate than its higher-energy counterpart. The e_{g} band tails exhibit a similar behavior, but they are initially split by a larger magnitude (83 meV). As such, both pairs of band tails give rise to similar features in the observed TDS, with the e_{g} absorption band notably wider than the t_{2g} band [Fig. 6(b)]. At all finite temperatures, the valence band tail and uppermost t_{2g} conduction band tail are separated from their respective continuum edges by ∼80 meV. We therefore propose the unique excitation profile of the 163-meV second-order Raman band [Fig. 4(b)] is the result of an additional double-resonance enhancement at near band edge excitation.

Real-space depictions of the charge density ($\u2211k\psi kc2$) of the conduction band tails at 294 K are shown in Fig. 6(c). Both pairs are strongly localized within the volume of a 2 × 2 × 2 supercell, with the e_{g} wavefunctions marginally more disperse due to hybridization with surrounding O 2*p* orbitals. The localization arises from a significant phonon-induced distortion of the octahedral ligand fields around isolated iron centers. Notably, the charge densities of the two t_{2g} band tail states are localized around two adjacent pairs of Fe atoms, indicating that these states may provide a pathway for transport of localized charges via a thermally activated hopping mechanism.

Finally, we demonstrate that population of the t_{2g} conduction band tail directly forms an electron small polaron. An excess electron introduced to a thermally distorted 2 × 2 × 2 supercell at 294 K fully localizes to the lowest-energy t_{2g} band tail. Following geometrical relaxation of the singly charged state, all atoms but those in the vicinity of the localized charge return to their pristine equilibrium positions. The polaron charge density is shown in the left side of Fig. 7(a) and is nearly identical to the charge density of the vacant band tail, with the amplitude more equally distributed across the two Fe centers. Additionally, the charge density is fully contained within the volume of a single primitive cell, indicating the formation of an electron small polaron. Conversely, an identical charge introduced to a pristine 2 × 2 × 2 supercell, with nuclei clamped at their equilibrium positions, delocalizes across the entire lattice [see the right side of Fig. 7(a)]. Therefore, we conclude that a thermal population of phonons is capable of inducing localized lattice polarizations that serve as nucleation sites for small polaron formation.

The displacement vectors associated with the relaxed polaron are shown in the inset of Fig. 7(c). The distortion extends over two octahedral Fe sites and is characterized by the compression of the Fe–Fe distance and the symmetric displacement of oxygen ions from the center of the distortion. We express the contribution of a particular phonon to the polaronic distortion (*C*_{mq}) according to Eq. (18), where Δ*τ*_{p} is the atomic displacement associated with the polaron and *d*_{pmq} are normalized phonon displacement vectors,

Here, the subscript *p* specifically denotes atoms within the polaronic distortion: two iron atoms and ten oxygen atoms. Figure 7(b) contains a plot of the phonon dispersion curve of hematite shaded to indicate the phonon branches that contribute most significantly to the polaronic distortion (*C*_{mq} > 0.5). For each phonon band *m*, Fig. 7(c) plots the sum of *C*_{mq} over all **q**-points as a function of phonon energy. As demonstrated in Figs. 7(b) and 7(c), the strongest contributions to the distortion arise from 81-meV phonons at the **Γ-** and **Z**-momentum critical points. These phonons account for displacements of the oxygen atoms in the distortion. The iron displacements are comprised of multiple low-energy optical phonons, the most prominent of which is the zone-center 31-meV phonon. Additional contributions arise primarily from the hybridized acoustic and optical branches (20–30 meV) at the bottom of the dispersion curve.

## CONCLUSION

We have demonstrated that, when combined, DFT+*U*+*J* and the method of vibrational averages can accurately reproduce the effects of electron–phonon coupling in hematite. Our results indicate that a thermal distribution of phonons gives rise to highly localized band tails at the edges of the t_{2g} and e_{g} conduction bands. We attribute these bands to strong lattice polarizations arising from population of the LO phonons above the 50-meV phonon bandgap. Following an electronic excitation, population of these band tails leads directly to the formation of an electron small polaron. The resulting distortion lies primarily along the coordinates of the 81- and 31-meV optical phonons. Notably, this distortion encompasses two adjacent Fe atoms, which implicates a potential pathway for electron polaron hopping, a crucial mechanism for photoconductivity in this material.

This work represents a first-principles theoretical description of carrier–phonon coupling in hematite that supports the conclusion of our previous work: Thermally activated optical transitions directly populate electron small polaron states in hematite. By computing an atomically precise description of carrier–phonon coupling, we are able to unambiguously determine which phonons contribute to the polaronic distortion. These assignments are confirmed by experimental resonance Raman and thermal difference spectra first reported in our previous work.^{44} We anticipate that, when extended to other transition metal oxides, this computational approach will reveal fundamental insights into mechanisms of polaron formation in optically excited states of these materials.

## SUPPLEMENTARY MATERIAL

Details of the convergence of the *U* and *J* parameters calculated from the linear response method, characterization of the phonon modes as LO, LA, TO, or TA, and additional resonance Raman profiles of first-order modes are provided in the supplementary material.

## ACKNOWLEDGMENTS

Financial support for this work was provided by the National Science Foundation under Grant No. CHE-2044462. The authors are grateful for the computational resources and technical support provided by the University of Rochester Center for Integrated Research Computing.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jacob L. Shelton**: Conceptualization (equal); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Kathryn E. Knowles**: Conceptualization (equal); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Supervision (lead); Visualization (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.