Exploring far-from-equilibrium ultrafast polarization control in ferroelectric oxides with excited-state neural network quantum molecular dynamics.

Ferroelectric materials exhibit a rich range of complex polar topologies, but their study under far-from-equilibrium optical excitation has been largely unexplored because of the difficulty in modeling the multiple spatiotemporal scales involved quantum-mechanically. To study optical excitation at spatiotemporal scales where these topologies emerge, we have performed multiscale excited-state neural network quantum molecular dynamics simulations that integrate quantum-mechanical description of electronic excitation and billion-atom machine learning molecular dynamics to describe ultrafast polarization control in an archetypal ferroelectric oxide, lead titanate. Far-from-equilibrium quantum simulations reveal a marked photo-induced change in the electronic energy landscape and resulting cross-over from ferroelectric to octahedral tilting topological dynamics within picoseconds. The coupling and frustration of these dynamics, in turn, create topological defects in the form of polar strings. The demonstrated nexus of multiscale quantum simulation and machine learning will boost not only the emerging field of ferroelectric topotronics but also broader optoelectronic applications.

INTRODUCTION

Ultrafast photoexcitation offers an exciting avenue to control materials by markedly modifying their potential energy surface through the creation of electron-hole pairs, which can allow for order of magnitude faster phase changes in comparison to traditional methods and access to hidden nonequilibrium phases (–). One material system of particular interest for nonequilibrium optical control are perovskite oxide–based nanostructures due to their wide range of desirable electronic properties from traditional ferroelectricity to the formation of complex magnetic and polar topological structures such as skyrmions and merons (–). These properties are usually dynamically controlled via traditional phase change–inducing methods such as heating and straining (, , –), which has made perovskite oxides ideal materials for nanostructured devices (, ).

Recent studies in nonequilibrium optical control of perovskite oxides have primarily focused on using weak terahertz pulses to couple with optical phonon modes to directly control polarization patterns (–). Above-bandgap laser pulse excitation with moderate fluence has also been shown to control nanopolarization topologies in SrTiO3 (STO)–PbTiO3 (PTO) superlattices (, ) through introduction of electron-hole screening at the material interface. However, topological polarization control under light-induced far-from-equilibrium conditions has been less explored. This is largely due to the fact that modeling the rich topological changes induced by nonequilibrium electron-ion dynamics under strong above-bandgap excitation requires understanding of multiple spatiotemporal scales, which has remained elusive.

To investigate far-from-equilibrium optical control of complex polar topologies, we have implemented a multiscale excited-state neural network quantum molecular dynamics (XS-NNQMD) framework to study strong optical excitation on the prototypical ferroelectric (FE) perovskite oxide PTO. In our multiscale framework, excited electron–lattice dynamics is studied by quantum molecular dynamics (QMD) simulations that combine Fermi-occupation QMD (FOQMD) to describe the effects of massive electronic excitations on atomic motions () and real-time time-dependent density functional theory (RT-TDDFT) along with Maxwell’s equations to describe light-matter interaction () and provide the initial excitation level for FOQMD simulation and surface hopping–based nonadiabatic QMD (NAQMD) to include nonadiabatic coupling between electrons and phonons (), thereby assessing the relaxation time between the electronic and lattice temperatures until which FOQMD is justified. Ground-state QMD and excited FOQMD simulations were then used to train neural networks to perform large-scale XS-NNQMD simulation to investigate atomistic mechanisms accountable for far-from-equilibrium control of large nanopolarization topologies in PTO. Despite remarkable success of machine learning–based MD simulations realizing atomistic simulations at scale while retaining quantum-mechanical accuracy (), NNQMD incorporating massive electronic excitations has been less explored.

We find that strong photoexcitation lowers the FE energy barrier, which results in a transverse optical (TO) phonon-driven transition to a nonpolar symmetric structure at room temperature. Concomitant to the loss of FE order under far-from-equilibrium excitation, a new order emerges in the form of tilting of the oxygen octahedra. This emergent order plays a key role in a wide variety of polar topological structures. Billion-atom XS-NNQMD simulations demonstrate how different regions in the crystal acquire different tilting orientations, which, in turn, result in frustration of the tilting phase to create topological defects in the form of string-like domain boundaries. Along these frustrated regions, the FE order is preserved and forms polar strings. This type of phase transition is analogous to nonadiabatic Kibble-Zurek phase transitions, where rapid temperature quenching in a uniformly symmetric field introduces string-like domain boundaries in space (, ), but with two key distinctions. First, nonadiabatic transition here is induced through photoexcitation-driven change of energy landscape rather than conventional temperature quench. Second, the dynamics resembles more that of an inverse Kibble-Zurek transition (), where the system is driven from an equilibrium antisymmetric polarized state to a highly nonequilibrium nonpolar symmetric state with topological defects emerging as result of the emergence of a hidden order parameter in the form octahedral tilting. Overall, our multiscale quantum simulation and machine learning framework reveals pattern/defect formation during highly nonequilibrium phase changes (–) applied to perovskite-based FEs but in the less explored regime of light-induced electronic excitation. The uncovered hidden parameter of light-induced octahedral rotations demonstrates the utility of the multiscale framework for exploring ultrafast control of FE materials for applications within the emerging field of FE “topotronics” (–).

RESULTS

To study the effect of strong optical excitation on PTO, we first performed FOQMD simulations, in which massive electronic excitation is represented by electronic occupation numbers that follow Fermi-Dirac distribution with an effective electronic temperature (). Following previous experimental and theoretical studies of nonthermal structural changes under strong photoexcitation (–), 5% of the valence electrons were promoted to the conduction band. We also performed RT-TDDFT simulations () to estimate the corresponding laser fluence (see section SI). We observed rapid thermalization of far-from-equilibrium electronic distribution within 2 fs due to electron-electron interaction, which justifies the use of FOQMD simulation. To further examine the effects of nonadiabatic coupling between electrons and phonons (), we also performed surface-hopping NAQMD simulations. Excited-state electron-phonon dynamics under NAQMD simulation was found similar to those in the FOQMD (section SII). The time scale of phonon-assisted electron-hole recombination is also examined with NAQMD, which was found to be on a much longer time scale than the MD simulations performed here, which is consistent with other studies (). This further justifies the use of the two temperature FOQMD simulation.

To study possible charge transfer under excitation, we performed static ground-state and Fermi-occupation DFT calculations (see Methods for details). Figure 1A shows the ground-state partial density of states (PDOS). Oxygen states are near the valence band maximum (VBM), whereas states near the conduction band minimum (CBM) are composed of Ti states, indicating that photoexcitation will cause charge transfer to from O to Ti.

Fig. 1.

Photo-induced charge transfer.

(A) PDOS. States near the VBM are primarily composed of O, and states near the CBM are primarily composed of Ti, indicating expected O─Ti charge transfer upon excitation. (B) Isosurfaces of hole (blue) and excited electron (red) densities upon excitation. The hole density is primarily localized around O atoms, while the excited electron density is localized around Ti atoms. Spheres with black, silver, and red color represent Pb, Ti, and O atoms, respectively. (C) Instantaneous force vectors as a result of excitation.

Figure 1B illustrates the expected charge transfer by plotting the difference in electron density between the excited and ground states. The hole density (blue) is primarily localized around O atoms, whereas the excited electron density (red) is primarily localized around Ti atoms. Further analysis of the charge transfer using Mulliken analysis is provided in fig. S5. This charge transfer results in an attractive force between O and Ti atoms along the c polarization axis, directly opposing the ground-state polarization, which is illustrated in Fig. 1C. The forces on the other atoms remained nearly zero due to the symmetry of the TiO6 octahedra.

To investigate the lattice dynamics under excitation, we performed FOQMD simulations in the isothermal-isobaric (NPT) ensemble at a temperature of 300 K on a 4 by 4 by 4 PTO supercell and computed the ground- and excited-state phonon dispersions (see Methods). Table 1 lists the average O, Ti, and Pb displacements along the c polarization axis for the first 60 and 500 fs of the NPT-MD simulation. The O atoms are divided into three categories O1, O2, and O3 in accordance to their positions in the five-atom PTO unit cell, which is labeled in Fig. 2A. The displacements during the first 60 fs along with the instantaneous force at time t = 0 illustrated in Fig. 1A indicate an initial activation of A1 TO1–like optical phonons that could potentially reverse or erase the polarization of the PTO crystal, which is diagrammed in Fig. 2A. Average displacements after 500 fs indicate up-conversion into the higher-frequency E TO2 mode, which is also illustrated in Fig. 2A. Collective thermal excitations of A1 and E TO modes have been previously deemed responsible for negative thermal expansion in PTO through c-axis contraction, and their softening has been related to polarization and tetragonality reduction in PTO (, ).

Table 1.

Average displacements of each atom type after 60 and 500 fs along the c polarization axis.

Oxygen types are labeled O1, O2, and O3 in accordance with their symmetry positions in the PTO unit cell, which is pictured in Fig. 2A.

Time	Pb	Ti	O1	O2	O3
60 fs	−0.04 Å	0.09 Å	−0.19 Å	−0.20 Å	−0.14 Å
500 fs	−0.14 Å	−0.26 Å	−0.67 Å	−0.68 Å	−0.52 Å

Fig. 2.

Bulk lattice and polarization dynamics.

(A) Schematic of the initially activated TO1 phonon mode and the up-converted TO2 phonon mode. Oxygen atoms are labeled O1, O2, and O3 with respect to their symmetry position in the PTO unit cell. (B) Snapshot of oxygen rotations with magenta arrow labeling of the rotation orientation for four-unit cells in ab plane 940 fs after excitation. (C) Ground-state and (D) excited-state VDOS. VDOS was computed from Fourier transform, the velocity autocorrelation function of the MD simulation at 300 K. (E) Ground-state and (F) excited-state phonon dispersion. Activated high-frequency optical branches are both softened and hardened to form a phonon bandgap. (G) Excited-state polarization dynamics at 300 K. Average Ti polar displacements along the c axis are initially reversed, followed by a steady-state structure with near-zero average polarization.

Visual inspection of the trajectory also indicates the presence of oxygen rotations. Figure 2B shows a snapshot of four-unit cells in the ab plane of the trajectory after 940 fs with the tilt orientation highlighted with magenta arrows. We also examined the phonon dynamics by comparing the vibrational density states (VDOS) of the ground-state and excited-state simulations, which are plotted in Fig. 2 (C and D, respectively). Figure 2D illustrates further phonon up-conversion beyond the initial TO1 and TO2 phonon activation. The VDOS also appears to form a gap between 10 and 15 THz upon excitation and loss of higher-frequency modes. To further examine this, we plotted the phonon dispersion for the ground state and excited state, which are shown in Fig. 2 (E and F). Hardening and softening of high-frequency optical branches leads to the formation of a large phonon bandgap. The presence of a large phonon gap has previously been attributed to quantum paraelectricity in STO (), further indicating the potential photoinduction of paraelectricity in PTO. Phonon instabilities also occurred at the R and M points (manifested as negative modes), which is consistent with octahedral rotations seen in the MD trajectory (, ). The unstable R and M eigenvectors corresponded to oscillation between in-phase and out-of-phase tilts, respectively, and a gif of the motion is provided in movies S1 and S2. To better understand the formation of octahedral rotations, we performed Mulliken bond overlap analysis of the Ti─O and Pb─O bonded atoms for the excited-state trajectory, which is plotted in fig. S6. In comparison to a ground-state QMD trajectory, weakening of both Ti─O and Pb─O bonds is seen in the excited state; however, the relative Pb─O overlap in comparison to the Ti─O overlap increases, indicating an increase in relative strength of the Pb─O interaction. This is consistent with the notion that tilting in ABO3 perovskites is often thought to result from a necessity of stabilizing short A-O interactions (). Previous experimental study of the ABO3 perovskite EuTiO3 (ETO) has shown a similar photoinduction of octahedral tilting (), and the observed cross-over between polar and octahedral tilting modes under photoexcitation is consistent with static DFT investigations of both PTO and BaTiO3 cubic structures under photoexcitation ().

To examine the polarization dynamics, we plotted the average polar displacements 〈〉 of the Ti atoms along the c polarization axis (Fig. 2G). 〈〉 was computed as the average difference of the Ti atoms from the centroid of the TiO6 octahedra. In a centrosymmetric nonpolar structure, Ti atoms will be located at the center of TiO6 octahedra, while in the polar structure, they are displaced. During the first 200 fs, we observe an immediate reversal of the average polar displacements, followed by relaxation to a quasi-steady state structure with near-zero average polarization along the c axis. For comparison, the polar displacements for ground-state QMD trajectory are shown in fig. S7, for which the crystal remains negatively polarized along the c axis. We also examined the time evolution of the ratio of the lattice constants and the average unit-cell volume, which are shown in fig. S8. We found the lattice to contract and the lattice ratios c/a and a/b to equilibrate at 1 ps, further indicating the formation of a cubic nonpolar phase. A video of the excited-state trajectory of the bulk structure with illustration of the unit-cell polarization vectors is provided in movie S3.

The observed phase change can be best characterized by the order parameter D, which describes the FE ordering of the system in terms of the Ti polar displacement, and ϕtilt, which describes the tilt ordering in terms of the tilt angle. A different ϕtilt can occur about each axis and thus, in general, constitutes three different order parameters; however, certain tilt phases are forbidden by symmetry considerations (). A detailed discussion of the tilting observed is provided in section SVI. We observe a loss of FE ordering, which is manifested in merging of the split peaks in the ground-state Ti─O partial radial distribution function g(r) as shown in Fig. 3A. The loss of FE ordering is countered by emergence of tilt ordering, which can be seen as a splitting of the third peak in the O─O─O bond angle distribution as shown in Fig. 3B. The change in ordering can be described by photoexcitation-driven lowering of the FE potential well barrier (Fig. 3C) concurrently with raising of the tilt potential well height (Fig. 3D). As the coupling of lattice motions is limited by the speed of sound within the material, this sudden onset of tilt ordering by photoexcitation in large nanostructures may result in topological defect formation as different regions acquire different tilt orientations with respect to the original tilt-free lattice.

Fig. 3.

Dual-order parameters of the phase transition.

(A) Ti─O partial radial distribution g(r) for both excited-state and ground-state trajectories. In the ground state, polarized structure Ti─O symmetry is broken, resulting in splitting of the first g(r) peak. The symmetry is restored upon optical excitation as the polarization is lost, resulting in a single peak. (B) O─O─O bond angle distribution for ground and excited states. Third peak splits in the excited state due to symmetry breaking from the optical induction of oxygen rotations. (C) and (D) diagram the phase change in terms of the FE-order parameter D and the tilt order parameter ϕtilt. Optical excitation lowers the FE energy barrier and raises the tilt energy barrier.

To study photocontrol of large nanostructures and potential formation of topological defects, we developed NNQMD models for both ground- and excited-state dynamics based on QMD and FOQMD training data, respectively, following workflow in Fig. 4A. Atomic positions from the FOQMD data are first used to compute highly scalable symmetry functions (see Methods) that are then fed into a neural network, which is trained to the total free energy of the quantum system including entropic contributions. Atomic forces are then computed by taking the derivative of the network and symmetry functions with respect to the atomic coordinate. The training data were generated in the canonical (NVT) ensemble at 300 K. While in the excited state, the NVT ensemble will artificially strain the PTO lattice, experimental boundary conditions such as tensile strain can replicate this condition. There was no major structural change observed differences between NVT and NPT QMD simulations. Details on the NNQMD model training and validation are discussed in Methods and section SVII.

Fig. 4.

Excited-state topological defect dynamics.

(A) XS-NNQMD workflow. (B) Initial configurations of bulk PTO crystal. Local polarization of each TiO6 structure is color-coded by their magnitude along c axis. (C and D) Photo-induced nucleation of TiO6 tilt domains and the formation of polar stripes in the bulk PTO system at 24 ps using XS-NNQMD simulation. The polar stripes are sandwiched by domains with different TiO6 cage tilt orientations. (E) Initial configuration of PTO crystal with circular antipolarization domains. The circular antipolarization domains were first introduced by displacing Ti-atom position and relaxed at low temperature using ground-state NNQMD simulation. Emergent Bloch-type skyrmion at each domain is shown in inset of (E) with the polarization vectors on ab plane before optical excitation. (F and G) Nucleation, propagation, and formation of polar stripes in the system using XS-NNQMD simulation. Unlike in the bulk PTO crystal system, the polar stripes consist of multiple segments of remnant polarization of the original PTO crystal and the circular antipolar domains. (H) Photo-induced domains and domain boundaries color-coded by the tilt orientation TiO6 cage as order parameter. (I) Atomic coordinates around the domain boundaries highlighting TiO6 cage tilt orientations in red. For clarity, only O-atom positions and lines connecting neighboring O atoms are shown.

We applied the validated NNQMD model to study optically induced phase-change dynamics in both a large (749 Å by 749 Å by 12.46 Å) uniformly polarized PTO crystal (Fig. 4, B to D) and a PTO supercell (354 Å by 354 Å by 12.46 Å) engrained with circular nanodomains of opposing polarization (Fig. 4, E to G). The initial condition of opposing polarization inside the circular domains seeded the creation of Bloch-type skyrmions in the ground-state NNQMD model as illustrated in the inlet of Fig. 4E. These types of Bloch skyrmions have been observed inside bulk PTO section of STO/PTO nanostructures (). While the initial condition of opposing circular polarized regions would usually require an interfacial material such as STO, recent experimental work has shown the ability to engrain arbitrary domains in FE crystals (). Thus, induction of Bloch-type skyrmions in a bulk PTO crystal is not out of the range of experimental possibility and provides a nontrivial ground-state topology whose dynamics under optical excitation can be compared to that of the optically excited dynamics in a single-bulk crystal. A brief discussion of the Bloch-type skyrmion formation is provided in section SVIII. Videos of the dynamics for both the bulk crystal and the Bloch-type skyrmion domains under optical excitation are provided in movies S4 and S5, respectively. The dynamics for both systems illustrate a wave-like nucleation of the nonpolar regions as optical excitation lowers the FE energy barrier; however, the polarization is partially preserved where the wave propagations meet in the form of polar strings as illustrated in Fig. 4 (D and G). As optical excitation destroys the original FE topology, both the uniformly polarized and Bloch skyrmion systems gain similar polar string topologies, indicating that this phase change is largely independent of the original PTO FE domain structure.

In addition, we observed that the polarization-erasing waves correspond to the oxygen rotations, whose dynamics for the bulk crystal are illustrated in movie S6. A snapshot of the rotation dynamics for the bulk crystal under optical excitation with each unit cell colored by the parity of the tilt orientation about the c axis is shown in Fig. 4H. The red and green colors represent the two mirrored rotation orientations with respect to original lattice (which are illustrated in the x axis of Fig. 3D), and the blue color represents the nontilted unit cells. While these two mirrored orientations belong to the same symmetry group, they can become frustrated if they meet at a boundary that is not commensurate to the periodicity of the tilt. Figure 4H clearly illustrates that polarization stripes occur at the meeting of two different tilt domains where the octahedral rotations are frustrated, creating a tilt domain wall. A zoomed view on the local atomic structure of the TiO6 cages in the striped region is provided in Fig. 4I, which further illustrates the frustration of the TiO6 cage rotation at the domain boundaries. The frustration topologically preserves the c-axis polarization in these regions as the FE and tilt phonon motions responsible for controlling these order parameters are nonlinearly coupled. The preservation of ferroelectricity under rotation frustration was also confirmed in an artificially frustrated PTO supercell with FOQMD simulation and is discussed section SIX. Similar frustration and symmetry breaking of different tilt domains has been seen at heterostructure interfaces with different tilt phases () and domain walls in improper FEs ().

The effect of simulation size on the formation of topological defects was also investigated with large XS-NNQMD simulations of bulk PTO involving up to a billion atoms with the length of the a and b axis of ~3 μm. After 7.5 ps after photoexcitation, domains on the order of ~10 to 100 nm were formed. Smaller simulations were found inadequate to represent such large-scale topological features. This was quantified through the computation of a structure factor of the unit cells composing the string defects, which is illustrated in section SX.

The wave-like propagation of octahedral rotations and formation of topological defects resembles Higgs condensation of a disordered field (). This type of transition has previously been observed at FE critical point of multiferroic hexagonal manganites (). In our case, instead of using rapid temperature quenching to induce the nonadiabatic phase change, photoexcitation nonadiabatically changes the energy landscape by markedly lowering the FE energy barrier while also creating a deep-tilt potential well. In turn, topological defects emerge in the form of frustration of the tilt domains and polar strings.

DISCUSSION

Our XS-NNQMD simulation incorporating massive light-induced electronic excited state allowed us to study ultrafast polarization control of PTO-based nanodomains. Simulation results reveal how intricate coupling between acoustic tilting and TO polarization controlling phonon modes under optical excitation can lead to unique topological structures in the form of polarization stripes. Photoinduction of octahedral tilts in PTO is similar to pervious experimental study of ETO (), and the observed cross-over of tilting and polar phonon dynamics under photoexcitation is consistent with previous static DFT study of the PTO cubic structure under excitation (). Coupling of octahedral tilting with polar phonon modes in perovskites and manganites has also been attributed to a similar form of spontaneous symmetry breaking in the form of Higgs and Goldstone–like phonon excitations (, ) that have potential analogs to the Kibble-Zurek–like phase transition observed here. Here, we would like to emphasize two aspects. First, in contrast to conventional Kibble-Zurek transitions (, ), where the system is driven from disordered to ordered phases, a transformation analogous to an inverse Kibble-Zurek transition () is seen here, where polarization is driven from an equilibrium-ordered phase to a disordered phase in highly nonequilibrium state, and string-like topological defects arise because of the emergence of a hidden order parameter in the form of octahedral tilting. The uncovered light-induced octahedral rotations demonstrate the power of this simulation framework for exploring ultrafast control of topological phase changes.

Last, the multiscale XS-NNQMD framework used here is not unique to PTO or FEs and can be applied to broader light-induced phase changes for which adequate training data can be generated. By incorporating time-dependent carrier densities through calculation of electron-hole recombination rates, the framework will allow for investigation of more extensive topological defects on the microsecond time scale. Overall, this work has introduced a powerful framework that incorporates the forefront of multiscale quantum simulation and machine learning for exploring nonthermal phase changes, with a specific application to PTO, where a unique topologically protected polar string phase was found. Further investigation of this topological control in PTO and other FE materials could lead to development of novel topological FE devices on ultrafast timescales.

METHODS

QMD simulations

QMD simulation follows the trajectories of all atoms, while computing interatomic forces from first principles in the framework of DFT calculations (, ). To incorporate electronic excitations, FOQMD simulations () were carried out using the highly parallelized plane wave–based software, QXMD (, ). A 4 × 4 × 4 PTO tetragonal supercell was simulated using the gamma point for Brillouin zone sampling. A plane-wave basis was used with a cutoff energy of 30 rydberg (Ry) (410 eV) for the wave functions and 250-Ry (3400 eV) cutoff for the charge density. Vanderbilt-style ultrasoft pseudo-potentials were used, and local density approximation (LDA) was used for the exchange-correlation functional (). MD simulation was carried out in the NPT ensemble using a time step of 1.206 fs.

Phonon calculation

Phonon dispersion was computed using DFT in the Vienna Ab Initio Software Package (VASP) (, ). Calculations were performed on a 2 by 2 by 2 PTO supercell in the tetragonal structure at the experimental lattice constants, and a 8 by 8 by 8 k-point grid using Monkhorst-Pack sampling was used to sample the Brillouin zone (). LDA was used for the exchange-correlation functional (). A plane-wave basis with an energy cutoff of 800 eV within the projected augmented wave vector method was used to calculate electronic states (, ). Ionic optimization was performed until forces were less than 5 × 10−4 eV/Å. The dynamical matrix was calculated using finite differences, which was then used to calculate the phonon band structure and density states using the Phonopy software ().

Neural network quantum molecular dynamics

NNQMD uses feedforward neural network architecture, which consists of an input layer, multiple hidden layers, and an output layer. The nodes between adjacent layers are connected by adjustable weight parameters, which are determined by ground truth QMD training dataset. Atomic trajectory from ground-state QMD and FOQMD was used to train NNQMD model to perform billion-atom XS-NNQMD simulations. The input for NNQMD is a multidimensional feature vector for each atom, which represents its local environment within a cutoff distance in a permutational, translational, and rotationally invariant fashion using radial and angular symmetry functions (, ). The feature vector is fed into a fully connected neural network, which was trained to predict the total free energy including entropic contributions of the system, and atomic forces are then computed by taking derivate of the network with respect to the atomic coordinates. The training and test root-mean-square error values well converges within a few milli-electron volts per atom, and obtained XS-NNQMD model is validated using radial distribution functions and bond-angle distribution. Full details on model generation, training, and validation are provided in section SVII.