Hiroki
Uratani
*ab and
Hiroshi
Onishi
cde
aDepartment of Molecular Engineering, Graduate School of Engineering, Kyoto University, Kyoto 615-8510, Japan. E-mail: uratani@moleng.kyoto-u.ac.jp
bPRESTO, Japan Science and Technology Agency, Kawaguchi, Saitama 332-0012, Japan
cDepartment of Chemistry, School of Science, Kobe University, Kobe, Hyogo 657-8501, Japan
dResearch Center for Membrane and Film Technology, Kobe University, Kobe, Hyogo 657-8501, Japan
eDivision of Advanced Molecular Science, Institute for Molecular Science, Okazaki, Aichi 444-8585, Japan
First published on 3rd July 2025
The water splitting reaction mediated by photocatalysts is attracting much interest in the context of solar energy utilization. However, our understanding of the charge carrier dynamics underlying the photocatalytic water splitting reaction is still limited. Here, focusing on the perovskite-type oxide NaTaO3, which is an archetypical heterogeneous photocatalyst for the water splitting reaction, the charge carrier dynamics were investigated by a computational approach, i.e., molecular dynamics simulations based on quantum chemical calculations. In particular, the present study sheds light on the formation process of polaron, which is a charge carrier dressed by lattice distortion in its surroundings induced to stabilize the charge carrier itself. The results suggested that the charge carriers are weakly localized and maintain the nanoscale spatial distribution of the charge density, and the change in O–Ta bond lengths is the primary factor in the polaron stabilization. In addition, the structural deformation and the resulting polaron stabilization was observed to be stronger in positive (hole) polarons than negative (electron) polarons.
Here, a computational approach was employed to track the ultrafast process of polaron formation in the perovskite-type oxide NaTaO3, an archetypical compound that serves as a homogeneous catalyst for the water splitting reaction.15 While NaTaO3 is studied frequently in the presence of dopants in photocatalytic applications,16,17 the present study focuses on the fundamental character of pristine NaTaO3. The calculations were conducted based on the combined method of quantum chemical calculations and molecular dynamics simulations, which enables us to develop the real-time simulations of atomic dynamics and the associated change in electronic structures simultaneously.
The structural model for the present simulations was constructed as follows. Using the geometry-optimized cubic unit cell of NaTaO3 (Pmm) as the starting point, an orthorhombic unit cell was constructed by rotating the unit cell by 45° around an (a or b or c) axis and repeating the cell twice along that axis. The structural model for the formation of polarons was constructed by further repeating the cell by four times for each axis, resulting in the model containing 256 units of NaTaO3 (Fig. 1).
![]() | ||
Fig. 1 Structural model of NaTaO3 for simulating polaron formation (after thermalization). The yellow and red spheres represent Na and O ions, respectively. A sphere located at the center of each octahedron represents a Ta ion. VESTA program18 was used for visualization. |
All the calculations, unless otherwise noted, were conducted with the self-consistent-charge density-functional tight binding (simply denoted as DFTB in the rest of the present paper) method,19 which is an approximate analog of density-functional theory (DFT) based on the tight-binding model with the introduction of empirical parameters to replace the computationally costly part of the standard DFT calculations. Since the accuracy of the DFTB results depend on the employed parameters, these were optimized using the ADPT program20 to reproduce the DFT (GGA-PBE21) level electronic structure and structural dynamics calculated using the VASP program package.22,23 Further details of the parameter optimization and the benchmark results are presented in ESI.† While the DFTB offers an orders-of-magnitude acceleration compared to the standard DFT calculations, the computational cost is still a bottleneck when treating the systems including thousands of atoms, because in principle the computational cost for the DFTB increases proportionally to the cube of the number of atoms. To treat the nanoscale model system within a reasonable computational time, an acceleration technique called divide-and-conquer (DC)24 was employed in combination with the DFTB. In the DC approach, the system is spatially divided into a set of small subsystems, and the quantum chemical calculation is performed for each subsystem with the inclusion of the neighbouring subsystems (“buffer region”) to mitigate the error originating from the spatial fragmentation. The subsystem electronic structures are summed up to reproduce the electronic structure of the entire system. In the present study, each orthorhombic unit cell was assigned as a subsystem, and the radius of the buffer region was set to 8 Å. The DC-DFTB calculations were performed under the periodic boundary condition using the DCDFTBMD program package.25,26
To simulate the polaron formation process, the Born–Oppenheimer molecular dynamics (BOMD) calculations were performed at the DC-DFTB level. In BOMD simulations, the atomic dynamics is propagated by integrating Newton's equation of motion step-by-step with a sufficiently small time interval Δt, where the force acting on atoms is obtained according to the electronic structure based on the DFTB calculations at each time step.27 Importantly, the BOMD simulations can capture the change in electronic structure, e.g., the fluctuation of electron density, and the dependence of atomic dynamics on that, both of which are key physics of polaron formation. In fact, the DC-DFTB-based BOMD simulation has been employed to study the polaron formation process in other solid-state systems and successfully provided detailed information on atomistic dynamics and electronic structure involved in the polaron formation.28 The present study follows the approach established in ref. 28 with necessary adjustment for NaTaO3. The simulations were conducted in the protocol as indicated below. First, to thermalize the system, the BOMD simulations were performed for 5 ps with rescaling the atomic velocity every 400 fs to set the temperature to 298.15 K. From the last 4 ps of the obtained trajectory, 81 snapshots were extracted every 50 fs. Using the each snapshot as the initial geometry and velocity, 81 runs of BOMD simulations were performed to simulate the formation of the hole polaron, and another 81 runs were performed for the electron polaron as well. The simulations of polaron formation were conducted by restarting the BOMD under the NVE ensemble starting from each snapshot with the total number of electrons in the system modified such that the total charge is +1 for hole polarons and −1 for electron polarons. For the simulations of polaron formation, the DC-DFTB calculations were performed in the spin-polarized manner.29 All the BOMD simulations were conducted with fixed lattice parameters. The time interval Δt was set to be 1 femtosecond.
Fig. 2 presents the spatial distributions of charge carrier densities for holes (yellow) and electrons (blue) at t = 0 and 500 fs in the simulations of polaron formation. As the charge carriers are represented by removal of an unpaired electron from, or, addition of it to, the neutral and spin-unpolarized NaTaO3 system, the charge carrier densities are represented by the total spin densities of the system. As shown in Fig. 2, the charge carriers are weakly localized; they have spatially inhomogeneous distributions, but are not completely localized to specific atoms or unit cells. Notably, Fig. 2 suggests that the nanoscale model structure employed here (Fig. 1) played a crucial role in observing the weak localization, because such long-range inhomogeneity in charge density cannot be properly described within smaller unit cells. The weakly localized distributions are attributed to the disorder in the structure arising from thermal fluctuations. Fig. 2 also indicates that the electrons delocalized more than the holes, being consistent with the fact that CBM of NaTaO3 has the larger band curvature than the VBM,30 making the effective mass of electrons less than that of holes. Note that the delocality of an electron observed here should be considered as an upper limit, because the DFTB tends to systematically overestimate the delocalization of electrons,31 due to the approximation in its theoretical backbone; the DFTB is grounded on the generalized gradient approximation, which is known to systematically overestimate the delocalization of polaron charge density because of the self-interaction error.9,32–34
![]() | ||
Fig. 2 Hole (yellow) and electron (blue) densities at t = 0 and 500 fs from representative trajectories of the simulations of polaron formation. |
The polaron formation induces spatial localization of charge carrier densities. The degrees of locality, i.e., radii of the hole and electron, were quantified by following the previously proposed approach by one of the authors.28 In this approach, the number of atoms over which the charge density is delocalized, d, is evaluated as where sA is the number of unpaired electrons populated on atom A evaluated with Mulliken population analysis. The index A runs over all atoms in the simulation cell, and sA is normalized as
. The number of atoms d can be translated into the total effective volume of the charge carrier density, V, as V = Vcelld/n where Vcell is the volume corresponding to a single formula of NaTaO3, i.e., 1/256 of the total volume of the simulation cell employed. For the cases of holes and electrons, n = 1 and 3, respectively, where n = 1 and 3 denotes the number of Ta and O atoms within a unit cell, respectively. The meaning of the definition of V can be interpreted as follows. Because of the characters of valence band maximum (VBM) and the conduction band minimum (CBM) of NaTaO3, holes and electrons in NaTaO3 are predominantly populated on O and Ta, respectively. Hence, d/n corresponds to “the number of unit cells” over which the charge carrier delocalizes, and the total “volume” of the charge carrier is obtained by multiplying it by the volume of the unit cell. Finally, the volume V can be translated into the effective radii r by treating the charge carriers as spheres, i.e., r = (3V/4π)1/3. Fig. 3 presents the timetraces of the charge carrier radii, indicating that the radii have significant fluctuations associated with the structural dynamics. The mean radii can be estimated as 0.49 and 1.14 nm for holes and electrons, respectively, quantitatively supporting the observation that the electron is more strongly delocalized than the hole.
In general, polaron formation is characterized by structural deformation induced by the charge carrier in its surroundings. Here, from the BOMD trajectories, the structural deformation was characterized by approach reported in ref. 28. Because the VBM and CBM of the present system are primarily located on the framework composed of B and C sites, i.e., Ta and O, we focused on the change in O–Ta bond lengths. In this approach, the bond length around the charge carrier rpol, which is considered to be involved in the polaron formation, is evaluated as the weighted-averaged bond length by the spin population of each atom. For instance, in the case of holes, where the charge carrier is considered populated on O atoms, where rAO–Ta denotes the bond length between the atom A, which is here assumed to be O atom, and the nearest-neighbor Ta atom. In the case of electrons, the same procedure as above was applied with swapping O and Ta. To clarify the effect of charge carriers on the structural deformation, “control experiments,” i.e., the BOMD simulations for a charge neutral system starting from the same initial geometry and velocity as the simulations of polaron formation, were also conducted. Fig. 4 shows the timetrace of O–Ta distance (rpolO–Ta) around holes (orange) and electrons (blue), and for comparison, the averaged O–Ta distance calculated for the charge-neutral BOMD simulations (black). Each result was averaged over the set of 81 trajectories starting from different initial geometries and velocities, and the error bars indicate 95% confidence interval of the averaged value. As indicated in Fig. 4, in the case of holes, the charge-weighted average was significantly longer than the neutral average at t = 0, and the charge-weighted average further increased for t > 0, finally arriving at a plateau at t = 50 fs. This behavior implies that the hole polaron forms by the following steps: the generated charge carrier is first localized to the region where the Ta–O bonds are long because of the thermal fluctuation, and then, relaxation of the localized charge carrier further increases the bond length. On the other hand, the O–Ta distance around electrons (blue) is only slightly longer than the charge-neutral result. Furthermore, the O–Ta distance around electrons did not increase with time, suggesting that in the present simulations the localization of electrons is dominated by thermal fluctuation and the additional structural deformation induced by the electron is insignificant compared to the thermal fluctuation. One of the reasons for this difference may be the different degree of spatial locality of holes and electrons, i.e., electrons are tend to be delocalized compared to holes (Fig. 3), weakening the local change in O–Ta bond lengths induced by electrons.
Note that in the previous study of polarons in a lead iodide perovskite (CH3NH3PbI3),28 the structural deformation induced by polaron relaxation was observed for both holes and electrons, while the reported radius of electrons (1.1 nm) was close to that in NaTaO3 (1.14 nm) estimated in the present study. The difference in the structural dynamics of electron polarons between CH3NH3PbI3 and NaTaO3 may be attributed to the difference in the nature of Pb–I and Ta–O bonds. Pb–I bonds in CH3NH3PbI3 have a partially covalent nature and the CBM of CH3NH3PbI3 is composed of antibonding orbitals between Pb and I,35 suggesting that the change in the Pb–I bond length strongly influences the spatial distribution and energy levels of the CBM orbital. On the other hand, Ta–O bonds in NaTaO3 are considered predominantly ionic, suggesting that such effects of structural deformation are weaker than in the case of CH3NH3PbI3.
The polaron formation is a relaxation process that induces the stabilization of a charge carrier owing to the structural deformation in its surrounding area. Fig. 5 indicates the timetraces of hole and electron energies. The hole energy (Eh) and the electron energy (Ee) are defined as the differences in electronic energy E of charged and neutral states at the same geometry R, i.e., Eh/e(R) = E(R, N ± 1) − E(R, N) where N is the total number of electrons in the system at the charge-neutral state. As shown in Fig. 5a, holes underwent a rapid and significant stabilization, where the hole energy reached a plateau after t = 50 fs. This is a qualitatively similar behavior to the timetrace of O–Ta distance around holes (Fig. 4), where the structural deformation was completed at t = 50 fs, suggesting that the relaxation of hole polarons is dominated by the elongation of Ta–O bond lengths. The hole polaron stabilization can be quantified to be ∼70 meV as the difference between the hole energy at t = 0 and that after t = 50 fs. On the other hand, the electron polaron energy (Fig. 5b) did not show a significant change with time, suggesting that the stabilization of electron polarons is not substantial. This result is consistent with the fact that the structural deformation associated with electron polaron formation was insignificant (Fig. 4).
![]() | ||
Fig. 5 Timetraces of hole (a) and electron (b) energies averaged over 81 BOMD trajectories. Error bars indicate 95% confidence intervals. |
The previous time-resolved (TR) infrared (IR) spectroscopy by Yamakata and co-workers36 suggests that the excited electrons in NaTaO3 have IR absorption ranging from 4000 to 1200 cm−1. They assigned the absorption to the excited electrons trapped in mid-gap states beneath the conduction bands, which are shallow enough to allow thermal excitation of trapped electrons to the conduction bands. Within the timescale of their measurements, trapping to the mid-gap states are suggested to occur instantaneously. They also found that the hole-originated absorption is absent, suggesting that the holes are also trapped to deep mid-gap states where thermal excitation of holes to valence bands is essentially not possible. Notably, Onishi and co-workers suggested the possibility that these mid-gap states originate from polaron self-trapping.37 While the underestimation of polaron binding energy is inevitable due to the methodological limitations of the DFTB, our results (Fig. 5) suggests that the hole polaron is more strongly self-trapped (i.e., having the higher polaron stabilization energy) than the electron polaron, which is qualitatively consistent with the above-mentioned experimental observation. Note that the polaron binding energies in NaTaO3 have been also calculated at the level of DFT with the hybrid exchange–correlation functional,9 providing binding energy values that are comparable to experimentally observed IR absorption energies. In addition, our simulation results suggested the polaron formation completes within 102 fs, which is consistent with the “instantaneous” formation of trapped electrons in the timescale of the TR–IR measurements, given that they originate from polaron self trapping.
In summary, quantum-chemistry-based MD simulations were employed to characterize the polaron formation process in perovskite-type oxide NaTaO3 from the viewpoints of charge density distribution, structural deformation, and energetics. In terms of charge density distribution, it was found that the charge carriers show weak localization to nanoscale spatial regions, which is attributed to the structural disorder induced by thermal fluctuation. The result also suggested that the localization of electrons is weaker than that of holes. The structural deformation was analyzed focusing on the change in O–Ta bond lengths associated with polaron formation. In the case of holes, the O–Ta bond lengths tended to be longer than average even at t = 0, and increased for t > 0, suggesting a two-step process of hole polaron formation; the hole is first localized to the region where the O–Ta bond lengths are incidentally long due to thermal fluctuations, and the localized hole further elongates the O–Ta bonds in its relaxation process. On the other hand, the O–Ta bond lengths around electrons were only slightly longer than average, and did not increase with time, suggesting that the electron polaron character is dominated by thermal structural fluctuations, and the additional structural deformation induced by polaron relaxation is insignificant. The difference in the structural dynamics of the electron polaron in NaTaO3 and in lead iodide perovskite (CH3NH3PbI3)28 implies that the nature of bonding between B-site and C-site ions critically influences the polaron dynamics. Finally, the observed polaron energetics were consistent with those expected from the above-mentioned discussion on the charge density distribution and the structural deformation. Holes underwent a relaxation of ∼70 meV, with a time trace approximately synchronized with that of the O–Ta bond length around holes. In contrast, the relaxation of electrons was insignificant, reflecting the fact that no substantial structural deformation was observed associated with the electron polaron relaxation process; however some attention should be paid to the overestimation of electron delocalization originating from the theoretical limitations of the DFTB technique.
This study was supported by JSPS KAKENHI Grant No. 22H00344 and 25H00902. H. U. is also grateful to JST PRESTO Grant No. JPMJPR23Q2 and JSPS KAKENHI Grant No. 21H05407 and 23H03976.
Footnote |
† Electronic supplementary information (ESI) available: Detailed procedure and testing results of the DFTB parameter optimization. See DOI: https://doi.org/10.1039/d5cp01859e |
This journal is © the Owner Societies 2025 |