The frequency-dependent polarization switching in nanograined BaTiO3 films under high-strength electric field

ABSTRACT The polarization reorientation in ferroelectric nanomaterials under high-strength AC electric fields is intrinsically a frequency-dependent process. However, the related study is not widely seen. We report a phase-field investigation regarding the dynamics of polarization switching and the electromechanical characteristics of a polycrystalline BaTiO3 nanofilm under applied frequency from 0.1 to 80 kHz. The grain boundaries and the in-plane strains are considered in the model. The obtained hysteresis and butterfly loops exhibit a remarkable variety of shapes with the changing frequency. The underlying mechanism for the observed frequency-dependent physical properties was discussed via domain structure-based analysis. In addition, we examined the influence of the kinetic coefficient in the Ginzburg-Landau equation as well as the influence of the electric-field amplitude to the frequency dependency. It was found that a higher value of kinetic coefficient or field amplitude tends to enhance the mobility of polarization switching and to transform high-frequency characteristics to low-frequency ones. GRAPHICAL ABSTRACT


Introduction
Ferroelectric nanofilms are receiving increasing attentions due to their outstanding dielectric property, tunable polarization state, and strong electromechanical coupling, which are essential for application in high-density capacitive energy storage [1,2], nonvolatile ferroelectric memory [3] and microelectromechanical system [4,5]. In particular, the functions of both binary data access and electrical energy storage are achieved via polarization reorientation driven by a high-strength electric field. These processes usually take place in an electrically dynamic environment -a single round of data reading/writing or discharging may only take several to tens of nanoseconds to complete. Hence, the physical response of ferroelectric nanofilms under pulsed voltage or high-frequency electric field is particularly significant in such application scenarios.
The frequency dependence of certain key properties of ferroelectric thin films has been investigated [6][7][8][9][10], e.g. the dielectric constant, the coercive field, the remanent polarization, and the piezoelectric coefficient, etc. In recent studies, attentions are given to the unique physical response of ferroelectric nanofilms under high voltage with high frequency [11][12][13]. It was numerically observed that the unique overall properties of ferroelectric nanofilms under a high-strength AC electric field are attributed to the dynamics of domain-wall motion as well as the frequency-dependent polarization switching [11,[14][15][16][17]. In the meanwhile, a number of special characteristics of the polarization domain structure were investigated for ferroelectric nanofilms. Aguado-Puente and Junkera [18] predicted the structure of ferromagnetic-like closed domains in BaTiO 3 ultrathin-film capacitors via first-principle simulation. Chen et al. [19] revealed various ferroelectric domain patterns in nanofilms due to mechanical loads and flexoelectricity. Jia et al. [20] experimentally verified the existence of vortex domain structures in nano-sized ferroelectric films. In addition, it is found that both the size-dependent vibration [21,22] and lattice misfit strain [23] have significant influences to the performance of nanofilm devices.
In retrospect, there are not many reports regarding the frequency-dependent hysteresis of ferroelectrics. One reason is the limited experimental ability for frequency-dependence testing on bulk materials which requires both high frequency and high electric-field strength. A few existing experimental observations, however, often exhibit points of inconsistency in terms of the frequency-dependent ferroelectric behaviors. Therefore, it is quite necessary to conduct theoretical and numerical investigation on this kind of problems, particularly for ferroelectrics with nano-dimension(s). A frequency-dependence test is more likely to be carried out on such kind of materials. There are very limited reports on the frequency-dependent microstructural dynamics, too, which may be remarkably influenced by the boundaries of the constituent nanograins and by the in-plane strain in the nanofilm. It is a crucial task to understand the underlying mechanism for the effect of electric-field frequency on the domain structures and the correlated physical properties of the nanofilms. With such purpose, we are carrying out phase-field modeling of the frequency-dependent polarization switching in a nanograined BaTiO 3 film under highstrength electric fields. In addition to the frequency dependency, the effect of the electric-field amplitude and the effect of the kinetic coefficient of the governing equation on the overall physical characteristics of the nanofilm are discussed, too.

Method and model
In this section, we describe the phase-field method that is going to be adopted in the modeling of the polarization switching of the nanograined BaTiO 3 films. The phase-field approach has been employed in the study of a wide spectrum of phase-transformation problems [24][25][26][27] and microstructure-evolution problems [14,15,[28][29][30]. The order parameters are set as the components of the electric polarization vector P i i ¼ 1; 2; 3 ð Þ, which was temporally governed by the timedependent Ginzburg-Landau (TDGL) equation, as @P i =@t ¼ À L ij @F=@P i i ¼ 1; 2; 3 ð Þ, where L ij are the components of the positive-definite kinetic coefficients. The total free energy of a system, F, was given by the volume integral of the energy density,ψ, as F¼ ð V ψdV. The free energy density function depends on the strain ε ij , electric displacement D i , the electrical polarization P i and its gradient P i;j . In addition, a general form of three-dimensional space free energy density ψ is given by Su and Landis [31]: The free energy of the system mainly comprises four parts: the first term corresponding to the ferroelectric gradient energy, the next four terms in the first brackets are used to create the Landau free energy, the following four terms in second brackets denote the elastic energy, and the final term represents the energy stored at the free space, and κ 0 is the permittivity of free space. We construct the formula of the virtual work based on the Maxwell equation, the equation of motion and the TDGL equation, and then solve them simultaneously by the finite element procedure. Interested readers are referred to these works [14,15,31] for a more comprehensive description of this approach and the selection of the material parameters. The structure of the studied BaTiO 3 ferroelectric film with nanograins is schematically depicted in Figure 1. The constituent nanograins are regarded as core-shell-type structure [32,33], which is composed of grain boundary (the white zone) and nanograin interior (the multi-color zone). The thin-film unit can be regarded as a representative part of the whole nanofilm material with numerous random grains. The periodic boundary conditions are exerted on the left and right sides of the thin-film unit. The nanograins with varying local crystal-lattice orientations are separated by the grain boundaries of low permittivity [34]. The dielectric permittivity of the grain boundary is set as 35% of the interior grain. The thickness of the grain boundary does not change with the grain size, and it is assumed to be 0.8 nm [35]. As illustrated in Figure 1, the local crystal lattice is along the local coordinates x L . The rotation angle between the local coordinates x L and the global reference coordinates x G is indicated by φ. Each grain is randomly assigned an orientation angle as φ ¼ 7:5 � n where n is randomly selected from (0, 1, 2, . . ., 11). Such specific configuration of the problem is different to the one discussed in the work by Su et al. [14], wherein a nanograined bulk BaTiO 3 ceramic is the subject of interest.
All calculations are performed under the assumption of a plane problem: the z-components are set to zero, i.e.
where � ε is the in-plane strain caused by the lattice mismatch between the nanofilm and the substrate [36,37]. The in-plane strain is assumed to be applied to the film uniformly. It's implemented by means of prescribed boundary displacements, as where l is the total length of the representative film segment. We assume that the top surface of the film is set free and the displacement is set to zero at the bottom. At the left and right sides of the nanofilm, the periodic boundary conditions of electric potential and polarization are given by ϕ 0; y ð Þ ¼ ϕ l; y ð Þ , P x 0; y ð Þ ¼ P x l; y ð Þ, andP y 0; y ð Þ ¼ P y I; y ð Þ, respecively. In order to characterize the ferroelectric hysteresis loops, an external electric field is applied along the where E max is the amplitude of the electric field, f the loading frequency, and t the real time. The amplitude of the applied electric field needs to be high enough to achieve complete polarization switching, we here set E max = 8E 0 , wherein the value of E 0 is 218 kV/cm [11,14,31]. To implement the alternating electric field in the thin-film model, we fix the electric potential to zero at the bottom of the film, and then assign the top surface an electric potential with the value of negative electric field multiplied by the thickness of the film. All the material constants used in this phasefield simulation are referred to the work of Su and Landis [31].

The frequency-dependent hysteresis behaviors
The hysteresis loops and butterfly loops were plotted in Figure 2 (a) and (b) with the frequency ranging from 0.1 to 80 kHz under the high-strength AC field (8E 0 ). The reference electric displacement, electric field, and strain are set at P 0 ¼ 0:26 C=m 2 , E 0 ¼218kV=cm, and ε 0 ¼ 0:82%, respectively. It can be found in Figure 2 that both loops exhibit conventional shapes when the applied frequency is relatively low (less than 20 kHz). Once the applied frequency is increased over 20 kHz, the hysteresis and butterfly loops start to change their shapes. The hysteresis loop turns to an ellipse shape at 20 kHz and above. The butterfly loop evolves to a kidney shape at 40 kHz and then an ellipse shape at 80 kHz. As a consequence, the frequency response of the BaTiO 3 polycrystalline nanofilms can be divided into two phases: 0.1-20 kHz for the low-frequency phase and 20-80 kHz for the high-frequency phase. Previous studies have also revealed such frequency-induced shape change. Liu et al. [38] experimentally found that the shape of hysteresis loops evolves from a thin rhomb-like loop to an ellipse-like loop as the frequency changes from 10 −5 to   [39] numerically simulated the dynamic hysteresis response of the ferroelectric system, and also found that the hysteresis loop in the high-frequency range becomes an ellipse-like loop. Yang et al. [6] experimentally observed the hysteresis behavior of PZT films with a thickness of 100 nm under an AC electric field of 50-2000 Hz. Notably, the electromechanical properties of the above materials are different from those of BaTiO 3 , but the general trend of their frequency-dependent ferroelectric behavior can be regarded as a basic reference, especially in the absence of experimental data on BaTiO 3 nanofilms.
The dielectric permittivity, the coercive field, piezoelectric constant, and the remanent polarization were calculated based on the hysteresis loops in Figure 2. Their values under the selected frequencies were intuitively summarized and provided in Figure 3. The dielectric constant κ r =κ 0 and the piezoelectric constant d 33 =d 0 33 are computed from the slope of the hysteresis loop and the slope of the butterfly loop at the zero electric field, respectively. From Figure 3, it can be found that both the dielectric coefficient and the piezoelectric coefficient monotonically decrease with frequency. Such observation is consistent with the results measured for nanocrystalline BaTiO 3 films [40][41][42]. Both the remanent polarization and the coercive field gradually increase with increasing frequency in the low-frequency region (yellow zone); while it decreases with increasing frequency in the highfrequency region (white zone). The rise-then-drop trend of the remanent polarization and coercive field were observed in multi-domained PZT films [38] and nanograined BaMn 3 Ti 4 O 14.25 polycrystals [43].
When the applied frequency is less than 20 kHz, the thin film can be fully poled in each cycle of electric-field loading. The hysteresis keeps the conventional look within this frequency range. Along the increase in the applied frequency, the polarization switching starts to fall behind the electric loading. That's why we observed stronger coercive field and higher remanent polarization (higher frequency causes severe delay in polarization switching). As the applied frequency further increases over 20 kHz, the thin film cannot get fully poled anymore. The non-fully poled hysteresis then shows the ellipse shape. In this scenario, the ellipse-shaped hysteresis becomes unsymmetrical about the horizontal axis, and higher applied frequency causes lower value of both the nominal remanent polarization and coercive field. Such observation well reflected our theory -the competition between the polarization switching and the electric-field loading leads to the frequency-dependent hysteresis as well as the corresponding material parameters.
In order to provide a better understanding on the origin mechanism of frequency dependence, some real-time calculations for the development of polarization versus the sinusoidal electric field were plotted. The results are shown in Figure 4(a,b), at the loading frequencies of 0.1 and 40 kHz, respectively. At 0.1 kHz, the development of polarization is seen to adequately keep up with the variation of the electric field. As the frequency increases to 40 kHz, the maximum of the P-curve (blue online) is reached long after the E-curve (red online) reaches its maximum, so the development of polarization becomes increasingly falling behind the electric field. Higher frequency leads to a larger phase difference between the polarization and the electric field. The magnitude of the coercive field, which is essentially the value of the electric field when the total polarization is zero, will increase with a stronger phase lag between the two curves. It reflects another point of view of the increased coercive field with the rising frequency. Such phenomenon can be explained through the hypothesis that the frequency dependence of the ferroelectric hysteresis behavior is the result of direct competition between the microstructure evolution and the external loading. In the lowfrequency range, it has sufficient time to evolve toward the equilibrium state. As such, a full 180° polarization switch can be realized during a complete cycle. This is also evident from Figure 4(a) that the amplitudes of the P-curve are equal in both positive and negative directions, revealing a complete 180° reversal. But in Figure 4(b) at 40 kHz, the amplitude of P-curve in the negative direction is lower than in the positive direction; it indicates that the electric polarization did not fully complete the 180° switching.

The underlying domain pattern
The origin of such frequency-dependent hysteresis behavior is revealed through the microstructure-based analysis of the phase-field simulation results in the nanograins. We plotted the domain structure under selected frequencies of 2 kHz and 50 kHz at the diverse state of loading electric field, as illustrated in Figures 5 and 6. At the remanent state, the electrical polarizations are all brought upward, and the film exhibits a monodomain pattern, as shown in Figure 3(a). When the electric field varies from negative maximum to the positive maximum, the polarizations experience a complete 180° switch, as demonstrated in Figure 3(c,e). The polarizations have enough time to reverse to saturation under a low-frequency loading, and the ferroelectric film reaches a fully polarized state. It is noticed in Figure 3(b,d) that the polarizations near the grain boundaries appear to be parallel to the boundaries, and they are particularly stronger than the polarizations in the grain interior at the coercive state. This is due to the fact that the polarization distribution near the grain boundary is nonuniform, collecting a large amount of polarization charge on the surface of the grain boundary and eventually generate a strong depolarization field.
Once the frequency is raised to a higher level the loading speed increases and the polarization switching starts to fall behind, and the remanent polarization decreases substantially. The domain structures under selected electric-field loading states at 50 kHz were plotted in Figure 6(a-d), respectively. It can be clearly observed that the polarization cannot undergo a 180° switching due to the excessive speed of electric loading. As shown in Figure 6(c), a part of the electrical polarizations underwent switching in the nanograins, while the other part did not change the direction (marked by bluedashed box) during a cycle of electric-field loading. The electric dipole-dipole interactions in each nanograin are stronger at high-frequency level, which inhibits the polarization reversal, resulting in a large coercive field. Since the electric dipoles did not have enough time to switch under the high-frequency loading, the changing rate of the polarizations is relatively small. It caused the gradual decrease in the dielectric coefficient and the piezoelectric coefficient with the increasing loading frequency.

The influence from the kinetic coefficient and the field amplitude
The kinetic coefficient relates the loading frequency to the dynamic response of the electrical polarizations of the nanofilm in the TDGL kinetic equation. In this work, the normalized time is written as t ¼ t=t 0 , where the characteristic time is given by σ 0 L 0 , and the characteristic stress is given by The value of the characteristic kinetic coefficient L 0 was determined such that a unit dimensionless time corresponds to 50 ns in real time, i.e. L 0 ¼0:002 C 2 � m 2 Ns. This value was adopted in our polycrystalline phase-field model and the results qualitatively agreed with the experimental observations [43]. It has also been successfully adopted in other phasefield simulations [44,45].
The speed of external loading is represented by the loading frequency, while the speed of polarization evolution is mainly characterized by the kinetic coefficient. Therefore, it is very essential to explore the relation between the kinetic coefficient and the loading frequency. As shown in Figure 7, four sets of data with different frequencies and kinetic coefficient L were selected. The kinetic coefficient is selected as 0.2 L, L, and 5 L, respectively, at the same loading frequency of 5 kHz. The fourth one is selected as 5 L at 50 kHz. It is found from the first three sets of data that the coercive field increases significantly with the decreasing kinetic coefficient at 5 kHz. When the kinetic coefficient decreases to 0:2L, the hysteresis loop becomes elliptical and the tail of the butterfly loop disappears. The results with lower kinetic coefficient are similar to the results at higher frequency. Obviously, a larger kinetic coefficient will raise the mobility of the domain structures and will transform the high-frequency results into low-frequency ones. Comparing the curves of 0:2L at 5 kHz and 5L at 50 kHz, we found that both exhibit the elliptic and kidney shapes of hysteresis and butterfly loops. In previous works, a constant kinetic coefficient was assumed in the phase-field models [46,47]. Therefore, it is reasonable to set the kinetic coefficient as a constant in the computation.
The results in Figure 7 were obtained with E max ¼ 8E 0 . We now consider the influence to the hysteresis and butterfly loops with increased value of the amplitude E max . To illustrate the influence of loading amplitude, we plot in Figure 8 the original curves with 8E 0 at 20 kHz, and two additional results, i.e. 16E 0 at 20 kHz and 16E 0 at 40 kHz. The results with 16E 0 at 20 kHz are found to be similar to the characteristics of the hysteresis loop with 8E 0 at 10 kHz, as shown in Figure 2. It implies that higher amplitude turns the high-frequency results into the low-frequency ones. The results in Figure 8 indicate that the stronger thermodynamic driving force given by a higher field amplitude will speed up the rate of polarization evolution. Such enhanced domain mobility could compete well with the highfrequency loading. However, the last set of results with 16E 0 at 40 kHz point to the fact that, even with a higher field amplitude, the provided driving force is not adequate to allow the domain evolution to catch up with the high loading frequency of 40 kHz, and thus the elliptic and kidney shapes of the hysteresis and butterfly loops are observed again.

Conclusion
We numerically investigated the frequency-dependence of the ferroelectric properties of nanograined BaTiO 3 films via phase-field modeling of the domain evolution and hysteresis behaviors under a high-strength AC field. The considered field frequency ranges from 0.1 to 80 kHz. A conventional hysteresis can be observed at 10 kHz and below. Both the hysteresis and the butterfly loops exhibit an ellipse and a kidney shape, respectively, once the applied frequency increases above 20 kHz. The remanent polarization and the coercive field show a rise-then-drop trend as the frequency increases from 0.1 to 80 kHz, while the piezoelectric coefficient and the dielectric coefficient monotonically decrease in this frequency range.
In the low-frequency range, the polarization has sufficient time to complete a full 180° switching, and thus the speed of polarization reorientation can catch up with the speed of the applied electric field. In the high-frequency range, however, a part of the polarizations could hardly complete a 180° reversing. In addition, a higher value of the kinetic coefficient of the governing equation tends to yield the results similar to the low-frequency ones. Conversely, a lower value of the kinetic coefficient tends to yield the high-frequency results. Likewise, increasing the value of E max , which effectively enhances the domain mobility, can compete well with high-frequency loading, but if the frequency continues to increase, the polarization would return to the unsaturated state again.
In summary of the observed frequency-dependent behaviors, we conclude that a risethen-drop trend of the remanent polarization and coercive field can be obtained in nanograined BaTiO 3 thin films as the applied frequency keeps increasing. Similar phenomena were observed in our past calculated results for BaTiO 3 nanoceramics, too. Such frequency-dependent behavior is crucial for design of related ferroelectric nano-devices, a comprehensive report can be found in the review by Scott and Gardner [48].

Acknowledgments
This work was supported by the National Natural Science Foundation of China under Grant No. 12172046.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
The work was supported by the National Natural Science Foundation of China [12172046]