Determining thermal activation parameters for ferroelectric domain nucleation in BaTiO₃ from molecular dynamics simulations

The mechanisms, which control polarization switching in ferroelectric materials, are not only of fundamental scientific interest but also of considerable technological relevance. This switching process, which involves the reorientation of electric dipoles within the material, is pivotal in applications ranging from nonvolatile memory devices¹ to piezoelectric actuators.² As the performance of ferroelectric materials is intricately tied to the kinetics and mechanisms of polarization switching, understanding the underlying atomic-scale parameters is of paramount importance for ferroelectric materials design.

Ferroelectric materials respond to electromechanical loading through two primary mechanisms: the growth of existing domains with energetically favorable lattice orientations or the nucleation of new domains, resulting in the emergence of complex domain structures. The interplay of these mechanisms is intricately linked to the activation enthalpies associated with the processes of domain nucleation and growth, both of which depend sensitively on the magnitude of the applied electric field.^3–5 

Understanding these activation parameters provides insights into the enthalpy barriers that must be overcome for polarization switching to occur. This allows to assess the contribution of thermal activation to the nucleation of domains, recognizing that the effective values of the macroscopic coercive fields may depend on the effects of atomic-scale thermal fluctuations, which can play a substantial role in initiating switching.^6,7 Moreover, tuning the activation parameters for domain nucleation, e.g., by addition of defects, might offer a pathway to modify the macroscopic behavior of ferroelectric materials.

Molecular dynamics (MD) simulations provide a powerful tool to investigate the process of polarization reversal, revealing insights into the stochastic nature of polarization switching driven by thermal fluctuations.^8,9 In this study, we focus on calculating the activation parameters that govern domain nucleation in BaTiO₃ at room temperature based on the analysis of coercive fields, obtained using MD simulations. From the statistics of coercive fields in repeated field sweeps, we determine activation energies and volumes for domain nucleation in perfect, defect-free crystals as well as in crystals with monovacancies (V_Ba, V_Ti, V_O), and first- and second-neighbor divacancies (V_Ba-V_O, V_Ti-V_O).⁹ To this end, we refer to MD data published previously by Tsuzuki et al.⁹ These MD simulations were carried out using a shell model interatomic potential developed in Ref. 10 that reproduces the multiple phase transitions of BaTiO₃ in good agreement with experimental data. We note that another widely used interatomic potential, ReaxFF, which captures bond breaking and formation, has been adapted for the study of ferroelectric behavior in BaTiO₃. While interesting for simulating the effects of vacancy re-orientation and migration,¹¹ this potential is not suitable for the investigation of thermal fluctuation effects at ambient temperature because it seriously underestimates the ferroelectric phase transition temperature of 393 K, instead predicting that above 240 K the system is in its cubic, non-ferroelectric phase.

Molecular dynamics simulations of polarization change processes in BaTiO₃ single crystals at temperature T = 300 K were performed using a constant field electric rate of $E ̇ = ± 0.05$ MV/ps. The field was aligned with one of the crystallographic axes of the BaTiO₃ unit cell (henceforth z axis). The details of the simulation setup and simulation methodology can be found in Ref. 9.

Hysteresis curves resulting from the simulations are depicted in Fig. 1, top. Polarization change was observed to initiate in strongly localized regions of columnar shape aligned with the applied electric field, see Fig. 1, bottom where such nuclei appear at point defects. In different field sweeps, the polarization was observed to switch at randomly varying values of the coercive field. Since the polarization switch occurs over a very narrow range of E values, we identify for simplicity the coercive field with the field at which a critical nucleus forms. The corresponding cumulative distribution is shown in Fig. 2(a). Re-plotting according to Eq. (4) shows an almost perfectly linear behavior [Fig. 2(b)], where the linear fit curve represents the actual data with a Pearson correlation coefficient r = 0.98 and a coefficient of determination $R 2 = 0.96$⁠. We conclude that the thermal activation theory provides a very good description of the statistics of coercive fields in the simulations.

From the parameters of the linear fit, we can now deduce important parameters of the activation process by means of Eqs. (2) and (5). With $T = 300 K$⁠, an estimated value of the attempt frequency $ν 0 = 10 13$ s⁻¹ and $P 0 = 0.127$ C/m², we obtain with $a = − 15.8$ and $b = 0.38 × 10 − 6$ m/V the activation parameters $V N = 6.23 × 10 − 27$ m³ and $H N 0 = 9.14 × 10 − 20$ J = 0.57 eV. Note that the obtained nucleus size is in good agreement with nucleus sizes observed in the simulations, see Fig. 1, bottom.

We now investigate how the picture changes in the presence of defects. Defects studied include monovacancies on the Ba, Ti, and two nonequivalent O sites as well as Ba–O and Ti–O divacancies on nearest and second nearest neighbor sites. Such defects lead to a local charge re-distribution within the defected crystal, and therefore may exhibit both a permanent charge and a permanent dipole moment. The investigated defects and the associated local charges and dipole moment vectors are compiled in Table I. For a visualization of the defect configurations, the reader is referred to Ref. 9.

Figure 3 shows two typical scenarios. For a system containing a Ti monovacancy, the coercive field data are slightly shifted to lower fields as compared to a defect free crystal (full line), indicating a lowering of the enthalpy barrier. This shift is approximately independent of the direction of the field sweep (upward or downward). For a first-neighbor Ti-O divacancy, on the other hand, we observe a strong asymmetry since there is a pronounced shift to lower coercive fields during an upward sweep, but only a very slight shift during a downward sweep.

To analyze the reasons for this behavior, we formulated a defect model, which considers the defect as a spherical inclusion of radius $r D$ carrying a homogeneously distributed charge q and dipole moment $p$⁠. We consider three different locations of the defect relative to the nucleus. For details of the model and the related calculations, the reader is referred to the supplementary material.

1. The nucleus is located away from the defect. In this case, the interaction between defect and nucleus is negligible and the energy of nucleation is unchanged, $δ H = 0$⁠.

2. The defect is located in the center of the nucleus [see Fig. 1(b), left]. The nucleation enthalpy is then shifted by

   $δ H C = − 2 s P 0 p D 3 k,$

   (6)

   where k denotes the dielectric constant of the material, and s the sign of the polarization change (s = 1 for an upward change during an upward field sweep, and s = −1 for a downward field sweep).

3. The defect is located at the endpoint of the nucleus [see Fig. 1(b), right]. The nucleation enthalpy is then shifted by

   $δ H E = − P 0 ( r D | q D | 2 k + s p D 6 k ) .$

   (7)

Table I contains enthalpy shifts calculated for the different scenarios (center/end location and upward/downward sweep) from the charges and dipole moments of the different investigated defects. We make the assumption that the actual activation enthalpy is controlled by the defect location that produces the largest downward shift, $δ H U = min ( δ H C U , δ H E U , 0 )$ for an upward sweep and $δ H D = min ( δ H C D , δ H E D , 0 )$ for a downward sweep. Values of the enthalpy shifts were evaluated using the values $r D = a / 2 = 2 × 10 − 10$ m, where a is the lattice constant, and $k = 1.15 × 10 − 9$ C/Vm as taken from the work of Zgonik et al.¹³ 

To understand why divacancies affect polarization reversal much more strongly than monovacancies, we first note that Ba and Ti vacancies are negatively charged, while O vacancies carry positive charge. Monovacancies, thus, have significant (positive or negative) net charge but only small electrical dipole moment. The opposite is true for Ba–O and Ti–O divacancies. If the dipole moment of such divacancies is parallel to the ferroelectric polarization change, their presence may significantly reduce the energy of a domain nucleus provided that the nucleus forms around the defect. Monovacancies embedded centrally into a nucleus, on the other hand, create no net energy change. A comparatively small energy reduction is only effected when they are placed near one of the nucleus endpoints such that their electrical field enhances the external field within the nucleus. The bottomline of our findings is, thus, that polarization reversal is mainly affected by the electrical dipole moments, rather than by the electrical charges of defects. Note that this observation may become even more significant when dealing with multiple defects: the effects of multiple divacancies embedded into the volume of a nucleus are additive, but the effects of end-located monovacancies, for which there is essentially only one favorable location, are not.

In this study, we have conducted a comprehensive statistical analysis of coercive fields derived from hysteresis curves obtained through MD simulations. Our investigation covered multiple field sweeps in both pristine and defected BaTiO₃ crystals, featuring monovacancies, first-neighbor divacancies, and second-neighbor divacancies.

We find that the statistics of coercive fields is consistent with results expected for the electric field-assisted thermally activated overcoming of localized energy barriers. The corresponding barrier configurations can be envisaged as critical nuclei of domains of reverse polarization. Our statistical analysis reveals key activation parameters for the formation of such nuclei. Moreover, we observe that defects in the form of mono- and divacancies can significantly reduce the activation barrier and, hence, the critical fields required for polarization reversal. In particular, divacancies, which carry a significant electrical dipole moment, introduce an asymmetry into the hysteresis curves as a switch of polarization in the direction of the defect dipole moment is energetically favorable as compared to a reverse switch. Monovacancies, on the other hand, lead to only a small, symmetric shift of the coercive fields to slightly lower values.

Our approach represents a paradigm change as compared to previous studies of the influence of defects on ferroelectric hysteresis. Works, such as that of Arlt and co-workers on internal bias fields¹⁴ consider the coercive field as the critical field at which a minimum of the macroscopic enthalpy function disappears. Here, by contrast, we envisage the change of polarization as a nucleation-and-growth phenomenon controlled by the field at which thermal fluctuations can create a critical domain nucleus. This difference has important consequences: in the traditional perspective, defects change the coercive field by modifying the macroscopic enthalpy function; their effect is, thus, dependent on their concentration. By contrast, in the present viewpoint, even single defects can have a significant influence on the coercive field as they facilitate domain nucleation; at the same time, kinematic considerations, such as the question of domain wall mobility, become much more important for understanding the switching process.

From a methodological viewpoint, the present approach based on the analysis of cumulative distribution functions has significant advantages over the consideration of coercive fields in terms of mean value and standard deviation. Thus, in the previous work of Tsuzuki et al.,⁹ because of the large scatter, a statistically significant influence of defects on the mean coercive fields could only be established for first-neighbor Ba–O and Ti–O divacancies, where the effects are largest. In our study, on the other hand, defect-induced energy shifts as low as $δ H ≈ 0.5 k B T$⁠, corresponding to energy differences of 0.015 eV or field shifts of about 1 MV/m, can be established in a statistically significant manner.

As it stands, our investigation focuses on the effects of single defects. At higher defect densities the cumulative action of many defects may become an essential feature that enables polarization reversal at strongly reduced fields. The present work, which has established quantitatively the nucleation energy shift caused by different types of single point defects, serves as an essential starting point for the study of many-defect phenomena. The importance of point defects in enabling domain nucleation is corroborated by experimental observations (for an overview, see Ref. 15), which indicate that, at fields above about 1 MV/m, domain nucleation is no longer confined to macroscopic defects or surfaces but may occur diffusely in the bulk of BaTiO₃ samples. At the same time, point defects may exhibit a dual role in polarization switching because, while facilitating domain nucleation, they may act at the same time as pinning centers that inhibit the motion of domain walls at high defect concentrations. This idea is consistent with the behavior of BaTiO₃ after vacancies are introduced by γ irradiation:¹⁶ with increasing irradiation dose, the coercive field first decreases (domain nucleation is facilitated) but then increases again (domain wall motion and thus nucleus expansion is impeded).

While these findings highlight the potential significance of defect density in influencing material behavior, it is essential to emphasize that a comprehensive understanding of the collective impact of defects demands a detailed investigation beyond the scope of our current study. Such an investigation, which must cover the dual aspects of domain nucleation and domain wall pinning by the collective action of point defects, remains an important task for further studies.