Ultrafast Polarization Switching in BaTiO3 Nanomaterials: Combined Density Functional Theory and Coupled Oscillator Study

The challenge of achieving ultrafast switching of electric polarization in ferroelectric materials remains unsolved as there is no experimental evidence of such switching to date. In this study, we developed an enhanced model that describes switching within a two-dimensional space of generalized coordinates at THz pulses. Our findings indicate that stable switching in barium titanate cannot be achieved through a single linearly polarized pulse. When the intensity of the linearly polarized pulse reaches a certain threshold, the sample experiences depolarization but not stable switching. Our study also reveals that phonon friction plays a minor role in the switching dynamics and provides an estimate of the optimal parameters for the perturbing pulse with the lowest intensity that results in the depolarization of an initially polarized sample.


Introduction
Developing non-volatile memory devices with fast writing and reading operations while minimizing power consumption is a challenge in information storage.However, traditional magnetic storage and flash may not be suitable for future fast devices due to their limited operation speed, which is in the milliseconds range.Thus, this challenge can only be addressed by utilizing different physical mechanisms for writing and reading bits.
A potential physical mechanism for write operation is magnetization switching by an ultra-short electromagnetic pulse of optical or THz range.2][3] Similarly, electric fields can be utilized for ultra-fast polarization switching in ferroelectric materials.Although this possibility has garnered significant attention, it has not yet been observed experimentally.The closest successful result to date, which involved reversible polarization change, was achieved by Mankowsky et al. in their work on lithium niobate. 4Other studies [5][6][7][8][9][10] have also explored the selective excitation of lattice vibrations under ultra-short optical or THz pulses, which is essential for achieving practical polarization switching.
The absence of a predictive model poses a significant obstacle to experimentally observing ultra-fast switching of electric polarization.Such a model could provide optimal pulse parameters and answer a series of questions, such as: which normal mode should receive energy injection, whether energy should be injected directly into the mode that leads to switching or another strongly coupled mode; whether it is beneficial to use a series of pulses; which pulse polarization is optimal for switching; whether pulse shape affects switching; and which ferroelectric material is best suited for ultra-fast switching of electric polarization, among others.
2][13] To calculate material constants of ferroelectrics as oxides and chalcogenides, first principles methods like Density Functional Theory (DFT) are often utilized. 14These methods are effective in determining the structure of stable polarized states, energy barriers, ions' effective charges, polarization values, and the phonon spectrum. 12,15-20Moreover, it is important to highlight that DFT calculations' results are highly dependent on the chosen exchange-correlation functional. 21assical molecular dynamics (MD) simulations enable the examination of ultra-fast polarization switching at an atomistic level 11 and even take into account domain behavior. 22e proposed model aims to investigate ultra-fast polarization switching in ferroelectrics.
The model utilizes a system of ordinary differential equations (ODEs) to represent the time progression of the generalized coordinates within a ferroelectric material's elementary cell.
Radiation interaction is included by incorporating a perturbation force within the ODE, which functions for a specific duration.The potential energy surface (PES) is obtained from DFT calculations.Barium titanate (BTO) is used as a test material in this research, as it is a well-studied, prototypical ferroelectric material.
The proposed model primarily builds upon earlier works, 12,23-26 where a similar approach was employed for polarization switching and structure changes driven by ultra-short pulses.
However, two significant modifications were introduced.First, instead of representing the PES in the form of Taylor's series, we directly interpolate PES using cubic splines.This is because switching results in substantial atomic displacement, leading to high numerical errors in Taylor's series.Second, in terms of generalized coordinates, we consider the polarization mode (q p ), which undergoes the switch, and the normal mode (Q IR ) where radiation is pumped.In contrast to previous studies, 12 both generalized coordinates were normal modes, this approach contradicts the fact that the potential must be scalar, independent of the crystal's symmetry (for more details, please refer to 27 ).The article is structured as follows.
In the methods section, we give details of calculating the PES and constructing the system of ODEs.The results and discussion section presents the data obtained for BTO, along with a discussion on metastable switching, effective friction, perturbation duration, and optimal frequency.The conclusion section provides general observations and recommendations for future experiments.

Computational Details
We take the experimental values of a material's unit cell and relaxing the atomic positions to obtain the equilibrium structure.Both ionic relax and calculations for phonon spectra and energies are carried out using the Vienna Ab initio Simulation Package (VASP) software package, [28][29][30][31] employing a plane-wave basis set.The projector augmented-wave (PAW) pseudopotential with a general gradient approximation PBE 32 and a cutoff energy of 600 eV is utilized in all calculations.Numerical integration over the Brillouin zone is conducted using an 8 × 8 × 8 k-point sampling with a Gamma-centered grid.
The phonon dispersion curves are calculated within the framework of Finite Displacements (FD) using the Phonopy code. 33All corresponding DFT calculations are executed for a perfect 2 × 2 × 2 supercell structure.After identifying the normal modes, the PES is calculated as a function of two independent normal mode generalized coordinates: q p (polarization mode) and Q IR (high-frequency mode).
The individual atomic displacements, associated with the generalized coordinate q p , can be expressed as: Here, U i represents the displacement of the i-th atom, while Z U i and Z D i denote the coordinates of the i-th atom in the direction of polarization, corresponding to equilibrium positions with positive and negative polarization, respectively.
The individual atomic displacements, related to the generalized coordinate Q IR are given by: where U i is the displacement of a i -th atom, m i -atomic mass, η IR i is the corresponding component of the normal mode dimensionless eigenvector.
The PES is interpolated using cubic splines 34 at points where DFT calculations are obtained, allowing us to define the PES continuously as V (q p , Q IR ).The dynamic behavior It is essential to note that the frequency is significantly high, allowing approximately 2 oscillations to fit within the 250 fs envelope.
of the coupled generalized coordinates is characterized by a system of associated nonlinear differential equations of motion: where γ represents the effective friction coefficient, and F (t) is the initial force exerted on the system due to external pulse perturbation.The integration of Eq. ( 3) is performed using the odeint library from the SciPy package. 34We assume F (t) takes the following form: where F 0 is the force amplitude, ω d is the perturbation's driving frequency (assumed to equal ω IR , unless stated otherwise), and τ is the pulse's time length.A graphical representation how the the perturbing force increases with increasing of pulse duration is illustrated in Fig. 1.

Results and Discussion
The ferroelectric state of BTO is present in the crystal structure featuring a lattice with the P4mm space group.We adopt the following experimental crystal unit cell parameters: a = 3.986 Å and c = 4.026 Å. 35 The primitive unit cell is composed of one barium atom, one titanium atom, and three oxygen atoms (refer to Fig. 2).This structure gives rise to 15 normal modes at the Gamma point, including three acoustical and twelve optical branches, which are of particular interest to us.The optical normal modes at the gamma point can be decomposed as Γ = 3A 1 + B 1 + 4E.The initial cubic symmetry Pm-3m of the paraelectric BTO crystal at 130 Celsius goes for transition to ferroelectric state through atomic displacements strictly along the c-axis into tetragonal P4mm symmetry. 36Consequently, the coupling between normal modes and the motion (q p ) responsible for polarization switching is likely to occur with normal modes that possess large c-axis components in their eigenvectors.
In BTO, these modes are 5, 9, and 11, corresponding to frequencies of 5.3, 8.8, and 14.1 THz.The excitation of only three low frequency modes allowed us to avoid the nonlinear coupling between low and high frequency modes that it is known to affect the polarization switching in such material when both are present. 37In this work, we chose to investigate mode 5 because it represents a typical frequency that can be achieved with modern powerful terahertz radiation sources avoiding the presence of second harmonics. 4gure 3: Barium titanate potential energy surface is illustrated in generalized coordinates (q p , Q IR ).The heat map displays energy in eV units, measured from the base value of potential energy at (0, 0).
The PES was computed in the space of two generalized coordinates (q p , Q IR ), with each representing collective displacement of all atoms in the unit cell (refer to eq. 1 and 2).The sampling for q p was performed in the range from -2.0 to 2.0 with a step of 0.05 in Å √ amu, while the sampling for Q IR was carried out in the range from -3.0 to 3.0 with a step of 0.01 in Å √ amu units, resulting in a total of 48000 static DFT calculations.The point representation of PES was interpolated using cubic splines for solving the systems of ODEs.
This method offers a more accurate representation of polarization switching compared to the Taylor series expansion, which is only effective in the local vicinity of the expansion point. 4,12PES cross-section (shown in Fig. 3) along the direction Q IR ∼ 0 Å √ amu enables the examination of the barrier obtained by linearly interpolating the system's atomic coordinates from an upward polarization state to a downward polarization state.For DFT calculations, the barrier height is found to be ∼ 12 meV, which is consistent with other calculations employing the PBE exchange-correlation potential .The time evolution of generalized coordinates under the influence of varying pulse amplitudes is depicted, with the trajectory of generalized coordinates on the potential energy surface shown as a green line.(A) When the perturbation amplitude is relatively small, no switching takes place, and the system remains at its initial minimum; (B) The switching may not be "stable" -the system can momentarily enter a state with opposite electrical polarization, but due to inertia, it may return to and remain in the initial minimum, preventing the switching from taking place; (C) If the perturbation amplitude is large enough, switching occurs, and the system transitions into a state with reversed electric polarization.
To analyze the trajectory of generalized coordinates under a perturbing pulse for differing perturbation amplitudes, a series of calculations was performed (refer to Fig. 4).The effective friction coefficient was set at µ = 0.04 f s −1 .Three distinct scenarios were observed: 1.When the perturbation force is not sufficient, the system remains in the initial minimum, with the trajectory localized nearby (see Fig. 4A).

2.
A scenario not typically addressed by other authors, 4,12 but worth noting, involves the system entering a different polarization state only to return to its initial state after a period of time due to inertia.Thus, even a strong enough perturbation impulse may not alter the final electric polarization (see Fig. 4B).
3. Upon reaching a specific threshold for perturbation amplitude, enough energy is transferred into the system to surpass the barrier between local minima, causing the system to switch to a state with reversed polarization (see Fig. 4C).
A reversible polarization switch was previously observed in a study 11 where lead titanate (PTO) was modeled at the atomic level.Therefore, exposing BTO to a single polarization pulse could lead to irreversible switching if the pulse parameters fall within a narrow range.
However, even with carefully chosen pulse parameters, irreversible polarization switching might not be achieved due to the chaotic nature of polarization switching. 38Further research is needed to investigate this hypothesis in detail.
A crucial fitting parameter in the equations that describe the dynamics of generalized coordinates is the friction coefficient.Estimating this coefficient can be done through calibration experiments.Nonetheless, several factors can impact the friction coefficient, such as: (1) the domain structure's dependency on geometrical dimensions of the ferroelectric material sample; (2) the influence of neighboring unit cells (not considered in this work); (3)   the density of local defects.As a result, we conducted calculations by varying the friction coefficient over a broad range, analyzing its influence on the threshold switching force and switching stability (refer to Fig. 5).Calculations were performed for three pulse duration: 250, 350, and 450 fs, and a set of friction coefficients ranging from 10 −3 to 10 −1 f s −1 .The calculations determined whether a switch occurred and if it was reversible or irreversible.
We observed that different pulse length only with friction coefficient between 10 −3 and Additionally, the frequency of the perturbing pulse was varied in the calculations (see Fig. 6).The lowest threshold amplitude was observed for frequencies in the range of 0.95 ω IR to 1.10 ω IR for the eigenfrequency of the perturbed Q IR mode, while the threshold force amplitude reduction was approximately 1.6 times.The optimal frequency shift observed results from coupling with the high-frequency mode and the presence of a friction term in the equation of motion.
An analysis of the motion equation reveals that for small amplitude excitations, 23,39 coupling effects cause renormalization of the optimal frequency, ω IR .A frequency shift in an underdamped oscillator is a well-studied phenomenon. 40t's also estimate the fluence corresponding to a typical force at which polarization switching occurs.We adopt the smallest noted value (see Fig. 6), which is on the order of F th = 2.5 × 10 −4 Å√ amu/f s 2 .This force (F th ) is equivalent to the acceleration a th = where Ω, S, and h represent the volume, surface area, and the length of the unit cell in the 'c' direction (see Fig 2a) respectively, and ϵ 0 is the vacuum permittivity.This simple estimation yields a value of F th = 250 mJ/cm 2 .Although this is a rather basic analysis, the derived estimation should be approached with caution.For comparison in article 4 studying lithium niobate (LNO), the onset of polarization switching occurred at fluences of 95 mJ/cm 2 , which is on the same order of magnitude as our F th estimate for BTO.

Conclusions
In this study, a model was examined and evaluated to characterize the ultrafast switching polarization in ferroelectric materials using BTO as test case.Analyzing the proposed model indicates that exist an operative range of the friction coefficient where the ultrafast switching polarization has the highest probability to happen.Such probability increases with increas-ing of pulse and the smallest threshold force amplitude necessary for switching is achieved within the range of 0.95 ω IR to 1.10 ω IR , where ω IR represents the normal mode frequency.
Polarization switching has been shown to be reversible, and it is probably a random process, meaning that slight changes in the perturbing pulse parameters might lead to an opposite final polarization.Thus, the complexity of the model in the future should include arbitrary polarization of the perturbing pulse, which may prove difficult to interpret, and possibility to consider multi-pulse cases.For example involving the depolarization potential, which is generated by secondary high-frequency pulses, which inject energy into the electronic subsystem raising the electronic temperature to tens of eV favoring the switching of polarization as seen in previously works.

Figure 1 :
Figure 1: A visual representation of the perturbing force F (t) is shown for two pulse durations (250 and 50 fs) and a frequency related to the high-frequency optical normal mode (5.3 THz).It is essential to note that the frequency is significantly high, allowing approximately 2 oscillations to fit within the 250 fs envelope.

Figure 2 :
Figure 2: (A) An atomic illustration of the tetragonal ferroelectric phase in barium titanate (BTO) P4mm is provided.Primarily, electric polarization switching is linked to the motion of the titanium atom along the c-axis: UP (U), same direction to the c-axis; NEUTRAL (N), no polarization; DOWN (D), opposite direction to the c-axis.The figure also illustrates the displacement patterns of the generalized coordinates denoted by q p and Q IR .(B)The energy barrier for BTO divides the two stable states related to the nominal downward and upward electrical polarization.The barrier's height, as calculated from first principles calculations, and it is approximately 12 meV, which agrees well with results from similar studies.

Figure 4 :
Figure4: The time evolution of generalized coordinates under the influence of varying pulse amplitudes is depicted, with the trajectory of generalized coordinates on the potential energy surface shown as a green line.(A) When the perturbation amplitude is relatively small, no switching takes place, and the system remains at its initial minimum; (B) The switching may not be "stable" -the system can momentarily enter a state with opposite electrical polarization, but due to inertia, it may return to and remain in the initial minimum, preventing the switching from taking place; (C) If the perturbation amplitude is large enough, switching occurs, and the system transitions into a state with reversed electric polarization.

Figure 5 :
Figure 5: A series of computations were performed for threshold amplitude of the perturbing pulse (F th ) from 0.0 to 0.3 Å √ amu / f s 2 at three different pulse lengths: (A) 250 fs, (B) 350 fs, (C) 400 fs.For each computation set, the friction coefficient (µ) was modified over a wide range of values, from 10 −3 to 10 −1 f s −1 .In each calculation, the presence or absence of polarization change was noted: red circles represent calculations where polarization switching did not occur; blue circles indicate instances where polarization shifted but eventually returned to its original state; and green circles denote calculations where the polarization switched to its opposite value.

Figure 6 :
Figure 6: The relationship between the amplitude of the perturbing pulse, which causes switching, and the frequency of the perturbing pulse is demonstrated.The graph indicates that the lowest threshold force amplitude falls within the range of 0.95 ω IR to 1.10 ω IR .A continuous line is included merely to serve as a visual guide.Pulse duration 250 f s.