Molecular dynamics study of ferroelectric domain nucleation and domain switching dynamics

Ferroelectric materials contain domains of ordered electric dipoles, separated by domain walls, that can undergo polarisation switching under externally applied electric fields. The domain switching dynamics in ferroelectric materials plays an essential role in their application to electronic and electro-optic de- vices. Previous studies suggest that the switching occurs largely through domain wall motion which is explained from the viewpoint of statistical physics on surface growth as the behaviour of a pinned elas- tic interface. We perform molecular dynamics simulations to investigate the domain switching process and quantitatively estimate the switching speed of anti-parallel 180° domains in ferroelectric, tetragonal BaTiO3 perfect single crystals at room temperature using the core-shell model. We observe an unprece- dented, non-linear increase in the domain switching speed caused by the nucleation of new domains within the switching domain. We determine the strength of the electric field to evoke nucleation of new domains and show that the nucleated domains diffuse into nearby favourable domains when the electric field is removed. Furthermore, we discuss the prominence of domain nucleations during ferroelectric switching.

Scientific RepoRts | 7: 806 | DOI: 10.1038/s41598-017-01002-0 where U is a characteristic energy barrier, k B is Boltzmann's constant, E C0 is a critical field at which depinning occurs at 0 K and μ is the dynamical exponent determined by the nature of the defects in the material sample. Equation 1 can be seen as a generalized formulation of the Merz's law 2, 10 , a where E a = UE C0 /(k B T) is the activation field and μ = 1. Due to rich and pronounced finite size effects, ferroelectric behaviour needs to be probed and understood at different length scales. There have been efforts to develop ferroelectric material models at the macroscale, using finite element methods based on non-linear continuum mechanics, and at the mesoscale, using phase field methods based on the Landau-Devonshire theory of phase transition [11][12][13] . However, such models essentially rely on empirical parameters and thus do not allow re-examination of assumptions in existing theories of domain switching dynamics. Furthermore, phase field methods need certain physical parameters, such as domain wall energies, domain wall thickness and domain wall speed etc, which are usually computed by atomistic simulations 14 .
At present, full first-principles calculations of ferroelectric crystals are limited to a few hundred atoms per supercell at zero Kelvin. This has stimulated interest in the development and use of effective atomistic models, from which a portion of the degrees of freedom have been integrated out 15 . Atomistic models of perovskite ferroelectrics, such as isotropic and anisotropic core-shell models 16 , approaches based on effective Hamiltonians 17,18 and bond-valence models 19 , have been developed and shown to be suitable to probe a certain range of length and time scales that are currently beyond the reach of first-principles calculations. At the atomistic length scale the electric dipole of a single unit cell divided by its volume gives the electric polarisation of the unit cell. Thereby using molecular dynamics the domain switching processes can be studied with atomistic detail by computing the polarisation of each unit cell in the simulation system at any instant of time.
While domain nucleations inside the reversal domain (as opposed to nucleations at the domain wall) have been observed in epitaxial Pb(Zr,Ti)O 3 and BaTiO 3 capacitors through experiments 20,21 , atomistic simulations with domain nucleations have not been observed or reported, and accounted for. We perform core-shell molecular dynamics simulations of domain switching processes in ferroelectric, tetragonal BaTiO 3 single crystals (see Methods). In particular we investigate the switching processes and the domain wall movement in uncharged anti-parallel 180° domains. We apply a range of electric fields and identify the strength of the electric field (E N ) necessary to evoke domain nucleation inside a domain. We find that the ferroelectric switching is characterized not just by the domain wall motion but more importantly by the nucleations and growth of new domains within the reversal domain. We explain the behaviour of domain switching speed in response to the applied electric field.

Results and Discussion
To investigate the ferroelectric switching processes in BaTiO 3 single crystals, we consider a set-up of anti-parallel domains under periodic boundary condition separated by 180° domain walls as shown in Fig. 1. The growing and the reversal domains are shown in red and blue colours, respectively. The growing domain has its polarisation along the +Z direction and the reversal domain along the −Z direction. The growing domain, as its name suggests, grows at the expense of reversal domain in the presence of an electric field energetically favouring the growing domain. In an atomistic simulation such anti-parallel domains set-up can be created by applying anti-parallel electric fields. Such poling electric fields are removed when the anti-parallel domain set-up is deemed stable. In other words, both the growing and the reversal domains undergo an equilibration process after which the domain walls remain stationary in the absence of external electric fields. It is imperative that the domains created by poling electric fields are stable (see Methods).
After the anti-parallel domains set-up is achieved an external electric field is applied along the +Z direction, favouring the growing domain. The ferroelectric switching process is then studied under various external electric fields. To visualize the switching processes, a fixed colour map with red to blue gradation is employed to indicate the Z component of the polarisation of each unit cell in the system, as shown in Fig. 2. Weak electric fields. When the magnitude of the externally applied electric field is around or below 20 MV/m (in the case of BaTiO 3 single crystals), the growing domain increases steadily and in a quantifiable manner. This is also to say that the domain walls gradually move into the reversal domain. The domain patterns obtained from the atomistic simulations at various instants during the ferroelectric switching process are shown in  Though the propagation of domain walls is identifiable, the contribution of domain wall propagation towards ferroelectric switching is considerably low. The domain patterns obtained from the atomistic simulations at various instants during the ferroelectric switching process under strong electric fields are shown in Fig. 4. Figure 4a shows four domain nucleations inside the reversal domain.
Two important inferences can be drawn upon studying the ferroelectric switching process under strong electric fields. First, the nucleation of new domains in the reversal domain appear primarily in the early stages of the switching process. Second, the nucleation of new domains create new domain walls which would further propagate through the reversal domain contributing to the switching process. Eventually however, apparent from Fig. 4d, the nucleated domains merge leaving behind fewer domain walls (or lesser domain wall area). The switching process occurs through propagation of the domain walls thereafter.
Since the switching processes under weak and strong electric fields are significantly different, one expects that the switching kinetics are also different. The switching kinetics can be obtained by pursuing the time-dependent evolution of the fraction of growing domains (or reversal domains) or by pursuing the change in overall polarisation during the switching processes. Figure 5 shows the evolution of the fraction of growing domains under a strong electric field, E = 30 MV/m. The initial stages essentially include overcoming the activation energy required for the nucleation of new domains. Nucleated domains and their growth are then responsible for speeding up the switching process until they merge. After most of the nucleated domains merge, the switching process slows down as the subsequent switching occurs primarily through propagation of the domain walls. On the contrary, the switching kinetics under a weak electric field remain almost linear.
It can be observed that the switching kinetics under a strong field resembles that of mono-domain switching. Experimentally, the switching kinetics of mono-domains have been explained using the Kolmogorov-Avrami-Ishibashi (KAI) model 22  Diffusion of nucleated domains. In an earlier theoretical study through atomistic simulations, it is concluded that the ferroelectric domain walls do not exhibit significant intrinsic inertial response 9 . In other words, the evolution of domain pattern revealed that the domain walls stop moving when the electric field is turned off and show that the domain wall speed is solely determined by the strength of the electric field at any given time. The simulations performed in this study are in agreement with the above under weak electric fields; on removal of the electric field the domain walls remained stationary. However, the evolution of domain pattern upon application and removal of strong electric fields is slightly different. This is due to the sustained nucleations of new domains.
The nucleated domains grow as long as the (strong) electric field is not removed. When the electric field is removed, the nucleated domains diffuse into the nearby nucleated domain or into the growing domain so as to reduce the total surface area of the domain walls as shown in Fig. 6. This shows that while there is no inertial response, there is an intrinsic response from the nucleated domains. The intrinsic response is such that it minimizes the domain wall energy in the system. Flat domain walls, as seen in Fig. 6c, indicate that the total interface area between the domains is minimized.

Domain switching speed.
To account for the influence of domain nucleations on the switching process under strong electric fields, a better descriptor of the switching process as opposed to the domain wall speed is necessary. Figure 7 shows the switching speed, which is defined as the ratio of transversal length of the reversal domain to the total switching time, under a range of electric fields. A large jump in the switching speed is seen between 20-30 MV/m. We attribute this jump in switching speed to the nucleation of new domains and their growth.
It is known that the density of domain nucleations increases with the applied electric field. We have not observed domain nucleations in simulations performed with the reversal domain of less than 20 unit cells (transversal length). Moreover, the critical electric field of 25 MV/m required to observe nucleations is specific to the reversal domain (30 unit cells) employed in this study. We hypothesize that in atomistic simulations of ferroelectric switching under weak electric fields, nucleation of new domains are unlikely to be observed due to the limited size of the reversal domain. Perhaps, it is possible to observe domain nucleations under a weaker critical field if one considers larger reversal domains for atomistic simulations. In other words the nucleation of new domains has pronounced effects on the polarisation switching dynamics in larger reversal domains. Further atomistic simulations are needed to study the size effects of the reversal domain on the critical field for nucleation of new domains.
Previously reported experimental values of the electric field to cause dielectric breakdown in BaTiO 3 single crystals are in the range of 50-365 MV/m [24][25][26][27] . Nucleation and growth of ferroelectric domains have been studied up to field strengths of at least 45 MV/m in bulk BaTiO 3 single crystals and 1.3 GV/m in thin films 3,25 . In this study, we observed domain nucleations at a much weaker field strength of 25 MV/m. We therefore rule out the  thereby. We establish that the domain nucleations can be observed, even within a limited size of the reversal domain, under strong electric fields using atomistic simulations. If the electric field is removed, the nucleated domains diffuse so as to minimize the total domain wall area. We reason that in atomistic simulations of ferroelectric switching under electric fields less than the established E N nucleation of new domains are unlikely to be observed due to the limited size of the reversal domain.

Molecular dynamics simulations.
We carried out isothermal-isobaric (NPT) molecular dynamics simulations to study the ferroelectric switching process in tetragonal BaTiO 3 single crystals at atmospheric pressure and close to room temperature (310 K). We use the isotropic, an-harmonic core-shell model potential for BaTiO 3 (presented below) that has been fitted to results of first-principles density functional theory calculations with a modified generalized gradient approximation 28 . All the simulations are performed with 60 × 10 × 10 unit cells under periodic boundary conditions. The thermostat and barostat relaxation times are both set to 0.1 ps. All shells are assigned a mass of 2 a.u. and their respective cores are assigned the remaining atomic mass 29 . The times step used is 0.4 fs. The growing and the reversal domains are generated through a two step process. The unit cells in the system are first poled using an external electric field for 4 ps. The left half of the system is poled in the +Z direction and the right half is poled in the −Z direction. The external field is then removed and the system is equilibrated for another 4 ps. Now, if the system is kept under no external bias, the growing domain and the reversal domain maintain the same average polarisation and the domain walls remain stationary. This ensures the stability of both the domains and the domain walls generated thereby. After the domains are prepared, the time (t) is set to zero and an external electric field is applied in the +Z direction. We then follow the switching process by computing the polarisation of each lattice unit cell in the system. In this study, we compute the polarisation of a unit cell using the core-shell charges 30,31 .
The core-shell model. The core-shell model has been used in molecular statics and molecular dynamics simulations to gain insights into the properties of ferroelectric materials [32][33][34][35] . In the core-shell model every ion is represented in terms of a charged (atom) core and a charged (electron) shell, linked by a harmonic or an-harmonic spring. This introduces electronic Polarizability in the ions. In the conventional core-shell model the shells are assumed to be massless. Because the shells have no mass, they respond instantaneously to the motion of the cores i.e., at any instant during a molecular dynamics simulation the shells are repositioned so as to be in the minimum potential energy configuration with respect to the current cores' configuration. In this work we adopt an alternative approach by assigning a small mass to the shells and treating their motion dynamically similar to that of the cores since this approach has been shown to be computationally more efficient 36 . We therefore treat cores and shells as point mass particles interacting classically through potentials. BaTiO 3 isotropic, an-harmonic core-shell model potential. The interaction between cores and shells is described by three different interaction potentials. The interaction between the core and the shell of an atom is treated exclusively by an isotropic, an-harmonic spring with the potential parameters k 2 and k 4 , The position vectors of particles i and j are given by r i and r j , respectively. The norm of the distance vector between two interacting particles i and j is given by |r ij | = |r j − r i |.
The interactions between shells of different atoms, which describe the electron cloud repulsion and the Van der Waals attraction, are modelled using the Buckingham potential, The variables A, ρ, and c denote potential parameters. Finally, the electrostatic interactions among the charged particles viz., cores and shells are considered through the Coulomb potential, The particle i's charge is given by q i , ε 0 is the electric permittivity of free space and the total number of particles is denoted by N. Note that there is no electrostatic interactions between an atom's core and its respective shell in the core-shell model. Figure 8 schematically illustrates the interaction between two atoms using the core-shell model.
Accumulation of the electrostatic interactions is one of the most computationally demanding tasks because of the long-range character of the Coulomb potential. Conventionally, Ewald summation is used to accumulate electrostatic interactions in a system with periodic boundary conditions with high accuracy and precision. For large system sizes, however, the Ewald summation method suffers from high computational costs. We therefore substitute the Coulomb potential by the Wolf summation approach, , α is the damping coefficient and R c is the cutoff radius to truncate the pairwise sum. The damping coefficient α and cutoff radius R c must be carefully chosen and calibrated such that the accumulated electrostatic interactions equals with that of the Ewald summation and that there is no artificial pressure developed in the system. To this end, we have studied the effect of the damping coefficient and the cutoff radius on the average Coulombic energy and pressure per ion, captured in the Fig. 9. A higher damping coefficient allows the use of a lower cutoff radius with the consequence of artificial pressure built-up in the simulation system.
It is worth mentioning that the only other molecular dynamics study to quantitatively estimate domain wall speed, to the best of our knowledge, was conducted using the bond-valence model potential 10,19 as opposed to the core-shell model potential, in PbTiO 3 perfect single crystals. The potential parameters used in this study are taken from ref. 23. They are obtained from performing least square minimization of the energy differences in the potential energy surface between the core-shell model and first-principles density functional theory calculations. The potential parameters are shown in Table 1 Table 1. BaTiO 3 core-shell model parameters used in this study 28 .