Thickness Dependence of Switching Behavior in Ferroelectric BiFeO3 Thin Films: A Phase-Field Simulation

A phase-field approach to the analysis of the thickness effects in electric-field-induced domain switching in BiFeO3 thin films has been formulated. Time evolutions of domain switching percentage for films with different thicknesses were explored to reveal the primary switching path and its dependence on film thickness. In addition, hysteresis loop for these films were calculated to obtain their coercive fields. Results show a nonlinear thickness dependence of coercive field for ultrathin films. A parametric study of the interactions between film thickness, coercive field, current-voltage (I-V) response, and polarization switching behavior is herein discussed, which could provide physical insights into materials engineering.


Introduction
Nowadays, progresses in materials science and engineering have shown increasing reliance on the advantages of materials properties that originate from multi-field couplings. In ferroelectric materials, these properties include electrical-mechanical coupling (i.e., piezoelectricity [1,2] and flexoelectricity [3][4][5]), electrical-thermal coupling (i.e., pyroelectricity [6]), magneto-electricity [7] coupling, dielectricity [8,9], and so forth. All of the above properties involve the dynamics of ferroelectric polarization as well as reactions to externally applied fields. Electric-field-induced polarization switching in ferroelectric materials has attracted extensive attention from experimental [10] and theoretical [11][12][13][14] perspectives, owing to their unique features which can potentially provide improvements in data storage [15,16], spintronics [17], tunnel junctions [18], and non-volatile magnetoelectric devices [19]. As less voltage would be needed to generate the same electric field in thin films than that in thick ones, the basic understanding of underlying physics for thickness effects in switching behaviors is of great importance to the applications of ferroelectric thin films in miniaturized devices.
One of the most studied ferroelectric materials, BiFeO 3 (BFO), simultaneously possesses ferroelectric properties and antiferromagnetic properties at room temperature, with Neel temperature T N ∼ 643 K and Curie temperature T C ∼ 1103 K [20][21][22][23]. The crystal structure of BFO below T C is a rhombohedrally distorted perovskite structure with its space group being R3c, [20] resulting in eight possible variants for spontaneous polarization in BFO, which can be denoted as R ± 1 = ±[111], R ± 2 = ± 111 , R ± 3 = ± 111 , R ± 4 = ± 111 . Polarization switching arises when an external electric field is applied to the film. Due to the specialty of rhombohedral structures, only three paths are symmetrically allowed for polarization switching. According to the rotation angle, these switching paths can be defined as 71 degree switching, 109 degree switching, and 180 degree switching, respectively. To discover the optimal conditions for BFO films in applications, Singh et al. successively studied the effect of film thickness, annealing treatments, processing temperatures, deposition substrates, and doping on the electric properties, particularly the polarization and leakage current density, of BFO thin films [24][25][26][27]. Lee et al. studied ferroelectric properties including polarization-electric-field (P-E) loop and leakage behavior for BFO films varying from 85 to 280 nm [28]. It was found that the coercive field decreased with increasing thickness. Similar conclusions were drawn by Tang et al. in a study about thickness-dependent multiferroic properties of BFO thin films, with film thicknesses ranging from 210 to 830 nm [29]. Theoretically, a thermodynamic model based on the Landau-Devonshire free energy expression and the Landau-Khalatnikov dynamic equation was employed to analyze the effects of film thickness and surface parameter on switching time and coercive field [30]. Results showed that the sign of surface parameter would determine the trend of coercive field changing with film thickness. Dead-layer effects in polarization switching has also been studied in ultrathin BFO films using both experimental methods and statistical analysis [31]. Chen and Lynch [12] proposed a micromechanics model to study polarization switching in tetragonal and rhombohedral structures and found that it is easier to move 90 • (for tetragonal) and 70.5 • (for rhombohedral) domain walls than 180 • domain walls.
In our previous work [14], we studied the process of electric-field-induced polarization switching with a phase-field model, and the results show that 71 degree switching is the primary switching path when BFO films are subjected to electric fields. To fundamentally understand the effect of film thickness on polarization switching in BFO thin films, we employed a phase field model to simulate the switching behavior in films with various thicknesses. Further, corresponding hysteresis loops (P-E loops) and voltage-dependent leakage current (I-V responses) were also calculated, revealing the relationship between film thickness and coercive field.

The Phase-Field Model
The present phase-field model describes ferroelectric domain configuration and evolution by polarization distributions P = (P 1 , P 2 , P 3 ). According to the Landau-Ginzburg-Devonshire theory, total free energy density for a ferroelectric single crystal is given by [14,32,33] (1) in which F Landau , F grad , F elas , F elec represent Landau energy, gradient energy, elastic energy, and electrostatic energy, respectively. Landau energy is expanded as a 4th order polynomial of the polarization components P i (i = 1, 2, 3): in which α 1 , α 11 , α 12 are Landau coefficients, and only α 1 is assumed to be temperature-dependent. Gradient energy, which represents the energy arising from the polarization gradient across domain walls, here is written as in which G 11 is gradient energy coefficient, and x i is the ith component of the Cartesian coordinate system. Elastic energy, which is related to elastic strains originating from phase transitions or lattice accommodations, is written as in which c ijkl is the elastic stiffness tensor. e ij = ε ij − ε 0 ij is the elastic strain, where ε ij is total strain, and ε 0 ij is stress-free strain, induced by spontaneous polarization via the electrostrictive effect: where Q ij is the electrostrictive coefficient. Electrostatic energy is obtained from (6) in which E i is electric field, and ε b is the background dielectric constant. Temporal evolution of the polarization distribution is governed by the time-dependent Ginzburg-Landau equation: in which L is the kinetic coefficient related to the domain-wall mobility. In the model setting, we consider an epitaxial grown BFO film on a rigid substrate in the (001) c direction, as shown in Figure 1.
Dimensions of the film are set as 128∆ × 128∆ × h f ∆, in which h f is the grid number for the film thickness, and ∆ is the grid size in the simulation. In the present work, ∆ is set as 1 nm. A periodical boundary condition is applied along the x 1 and x 2 axes. A thin film mechanical boundary condition is adopted in this simulation [33]. As for the electric boundary condition of the film, a short circuit condition is applied in this model [34]. The top and bottom electrodes used in this model are shown in Figure 1. It is believed that electrode materials also have effects on the properties of films, because of their differences in chemical and thermal stabilities, especially the Schottky barrier height, which has an influence on leakage current [35][36][37]. However, in this work, with the purpose of generally studying the thickness effect in electric-field-induced multi-domain switching, we did not take the effect of electrodes as the primary concern of this study.
On the film-substrate surface, we have and on the top surface, we have Electric field can be obtained through the following: . The values of the coefficient for BFO in the simulation are as follows [14,32]: α 1 = 4.9 × (T − 1103) × 10 5 C −2 m 2 N, T = 298K, α 11 = 6.5 × in which is electric field, and is the background dielectric constant. Temporal evolution of the polarization distribution is governed by the time-dependent Ginzburg-Landau equation: in which is the kinetic coefficient related to the domain-wall mobility. In the model setting, we consider an epitaxial grown BFO film on a rigid substrate in the (001) direction, as shown in Figure  1. Dimensions of the film are set as 128Δ 128Δ ℎ Δ, in which ℎ is the grid number for the film thickness, and Δ is the grid size in the simulation. In the present work, Δ is set as 1 nm. A periodical boundary condition is applied along the and axes. A thin film mechanical boundary condition is adopted in this simulation [33]. As for the electric boundary condition of the film, a short circuit condition is applied in this model [34]. The top and bottom electrodes used in this model are shown in Figure 1. It is believed that electrode materials also have effects on the properties of films, because of their differences in chemical and thermal stabilities, especially the Schottky barrier height, which has an influence on leakage current [35][36][37]. However, in this work, with the purpose of generally studying the thickness effect in electric-field-induced multi-domain switching, we did not take the effect of electrodes as the primary concern of this study.
On the film-substrate surface, we have and on the top surface, we have Electric field can be obtained through the following:

Results and Discussion
To study the effects of film thickness on the domain switching pathways in the BFO thin film, we performed three cases of different film thicknesses, where h = 32, 16 and 8 nm, using the above phase-field model. The BFO films are assumed to be epitaxial (001) oriented on a rigid substrate. Firstly, the system was evolved, without any external electric field, from a paraelectric phase to an energy minimizing domain configuration in a ferroelectric rhombohedral phase. Then, the electric field was uniformly applied as 500 kV/cm to the pre-evolved film in the direction perpendicular to the film surface for all three cases. The switching behavior was characterized by switching percentage, which was defined as the ratio of the number of grids that took switch action in a period of time via one specific switching path to the total grid number in the system. These percentages were calculated every 100 steps.
In all three cases, it was observed that domains switch to electric-field-favorable orientations mainly through a 71 degree switching path and that 180 degree switching is barely observed due to

Results and Discussion
To study the effects of film thickness on the domain switching pathways in the BFO thin film, we performed three cases of different film thicknesses, where h = 32, 16 and 8 nm, using the above phase-field model. The BFO films are assumed to be epitaxial (001) oriented on a rigid substrate. Firstly, the system was evolved, without any external electric field, from a paraelectric phase to an energy minimizing domain configuration in a ferroelectric rhombohedral phase. Then, the electric field was uniformly applied as 500 kV/cm to the pre-evolved film in the direction perpendicular to the film surface for all three cases. The switching behavior was characterized by switching percentage, which was defined as the ratio of the number of grids that took switch action in a period of time via one specific switching path to the total grid number in the system. These percentages were calculated every 100 steps.
In all three cases, it was observed that domains switch to electric-field-favorable orientations mainly through a 71 degree switching path and that 180 degree switching is barely observed due to the relatively high energy bias, as also reported in our previous work [14]. Figure 2a presents the typical time evolutions of the switching percentage for three paths in the 32 nm BFO film. It could be seen that, in the switching process, 180 degree switching, which consumed the most energy, remained at zero, and the percentage of the 109 degree switching path was even lower than 1%. The 71 degree path was the dominant one for polarization switching in the rhombohedral BFO films. To study the effect of thickness, switching percentages for this primary switching path in films with different thicknesses were compared (Figure 2b). For all samples, with an application of the constant electric field, polarization switching percentages burst into a maximum value in the beginning and drop to zero as time goes on. Films with different thickness have different decreasing tendencies, indicating different switching times.
Appl. Sci. 2017, 7, 1162 4 of 9 the relatively high energy bias, as also reported in our previous work [14]. Figure 2a presents the typical time evolutions of the switching percentage for three paths in the 32 nm BFO film. It could be seen that, in the switching process, 180 degree switching, which consumed the most energy, remained at zero, and the percentage of the 109 degree switching path was even lower than 1%. The 71 degree path was the dominant one for polarization switching in the rhombohedral BFO films. To study the effect of thickness, switching percentages for this primary switching path in films with different thicknesses were compared (Figure 2b). For all samples, with an application of the constant electric field, polarization switching percentages burst into a maximum value in the beginning and drop to zero as time goes on. Films with different thickness have different decreasing tendencies, indicating different switching times. Here, we plotted the detailed volume fractions for upward variants R + (R , R , R , R ) and downward variants R -(R , R , R , R ) in these samples in Figure 3. As the system evolved under an applied downward electric field, energy-unfavorable variants (in this case, the upward variants R + ) disappeared over time. We used the vanishing point of R+ as the end of electric-field-induced polarization switching. In the 32 nm film, the volume fraction of R+ dropped to zero around Step 1000. In the 16 nm film, this drop occurred at Step 3000, and, in the 8 nm film, at Step 7000. With the decrease of film thicknesses, the time that was needed for the system to finish the switching process became longer.  Here, we plotted the detailed volume fractions for upward variants R + (R + 1 , R + 2 , R + 3 , R + 4 ) and downward variants R -(R − 1 , R − 2 , R − 3 , R − 4 ) in these samples in Figure 3. As the system evolved under an applied downward electric field, energy-unfavorable variants (in this case, the upward variants R + ) disappeared over time. We used the vanishing point of R+ as the end of electric-field-induced polarization switching. In the 32 nm film, the volume fraction of R+ dropped to zero around Step 1000. In the 16 nm film, this drop occurred at Step 3000, and, in the 8 nm film, at Step 7000. With the decrease of film thicknesses, the time that was needed for the system to finish the switching process became longer. , 7, 1162 4 of 9 the relatively high energy bias, as also reported in our previous work [14]. Figure 2a presents the typical time evolutions of the switching percentage for three paths in the 32 nm BFO film. It could be seen that, in the switching process, 180 degree switching, which consumed the most energy, remained at zero, and the percentage of the 109 degree switching path was even lower than 1%. The 71 degree path was the dominant one for polarization switching in the rhombohedral BFO films. To study the effect of thickness, switching percentages for this primary switching path in films with different thicknesses were compared (Figure 2b). For all samples, with an application of the constant electric field, polarization switching percentages burst into a maximum value in the beginning and drop to zero as time goes on. Films with different thickness have different decreasing tendencies, indicating different switching times. Here, we plotted the detailed volume fractions for upward variants R + (R , R , R , R ) and downward variants R -(R , R , R , R ) in these samples in Figure 3. As the system evolved under an applied downward electric field, energy-unfavorable variants (in this case, the upward variants R + ) disappeared over time. We used the vanishing point of R+ as the end of electric-field-induced polarization switching. In the 32 nm film, the volume fraction of R+ dropped to zero around Step 1000. In the 16 nm film, this drop occurred at Step 3000, and, in the 8 nm film, at Step 7000. With the decrease of film thicknesses, the time that was needed for the system to finish the switching process became longer.  The difference in the performance of domain switching between films with distinct thicknesses can be attributed to the influence of the coercive field. To study the relationship between film thickness and coercive field, ferroelectric hysteresis loops for films with different thicknesses-h = 32, 16, 8, 4 nm-were calculated. Selected corresponding domain configurations were presented as well.
Hysteresis loops were calculated by applying dynamic electric potentials on top of the film, following the order 0 V → 5 V → −5 V → 5 V . Consequently, electric fields were generated in directions downward/upward perpendicular to the film surface. Applying positive potentials on top of the film would result in a downward electric field, which further gave rise to downward polarizations (represented with negative signs in Figure 4). On the contrary, negative potentials led to the expansion of upward variants at the expense of downward ones. Domain configurations close to saturation polarizations (Points A and C) and coercive fields (Points B and D) were shown alongside the P-E loops. At Saturation Point A, domains were mainly formed by the four downward variants: As the period of electric field gradually decreased from positive maximum to zero, polarization dropped slowly, suggesting high remnant polarization, especially in thick films. When the reverse electric field grew, particularly near the coercive field (Point B), opposite variants started to emerge around the domain walls, and their areas expanded over domains with the strengthening of electric potential, which was the result of the electric-field-induced polarization switching. As opposite polarizations reached a maximum (Point C), domains were mainly formed by the four upward variants: In the process of a gradually changing electric potential from −5 V (Point C) to 5 V (Point A), polarizations in films switched from upward to downward via the turning Point D. As film thickness decreased from 32 nm in Figure 4a to 4 nm in Figure 4d, domain width decreased remarkably, which is a size effect-also reported in the literature [38,39].
Another important effect of the thickness is the changed coercive field. In Figure 4, it is clear that the coercive potential, at which total polarizations become zero, decreases as film thickness decreases. In order to obtain the critical effect of thickness, coercive field is derived by dividing the coercive potential by the film thickness, as follows: E c = U c h , in which E c is the coercive field, U c is coercive potential obtained from the P-E loop, and h is the corresponding film thickness.
From Figure 5, the relation between coercive fields and film thicknesses can generally be expressed as E c ∝ h − 3 2 , which has a similar trend as the well-known Kay-Dunn law (E c ∝ h − 2 3 ) [40], but a greater value in power. This is because the Kay-Dunn law was derived from films with millimeter scale thickness, in which a switching process originated from nucleation. Yet in our work, the systems we studied are ultrathin films, and the polarization switching mechanism can be quite different due to the size effect [31,41]. We can see from the calculated domain configurations in Figure 5 that the domain width of the 4 nm film is much smaller than the width of the 32 nm film, which means that there are more domain boundaries in the 4 nm film. As seen in the figure, domain boundaries would block the process of domain switching, so thinner film will have a higher coercive field. The results that the coercive field decreases with increasing film thickness, and that thicker films have faster switching processes, confirm the general relation between switching time and coercive field [42]: t s ∝ 1 E−E c , where t s is switching time, E is applied electric field, and E c is corresponding coercive field.
Additionally, I-V response is also an important character for ferroelectric films when they are subject to electric potentials. Therefore, we further studied the leakage current evolution under varying applied potentials [43]. The calculated I-V response for films with different thicknesses are shown in Figure 6. The leakage current density in BFO films increases markedly in samples with small thicknesses, indicating the limitations of application for ultrathin films.
means that there are more domain boundaries in the 4 nm film. As seen in the figure, domain boundaries would block the process of domain switching, so thinner film will have a higher coercive field. The results that the coercive field decreases with increasing film thickness, and that thicker films have faster switching processes, confirm the general relation between switching time and coercive field [42]: ∝ , where is switching time, is applied electric field, and is corresponding coercive field.  Additionally, I-V response is also an important character for ferroelectric films when they are subject to electric potentials. Therefore, we further studied the leakage current evolution under varying applied potentials [43]. The calculated I-V response for films with different thicknesses are shown in Figure 6. The leakage current density in BFO films increases markedly in samples with small thicknesses, indicating the limitations of application for ultrathin films.  Additionally, I-V response is also an important character for ferroelectric films when they are subject to electric potentials. Therefore, we further studied the leakage current evolution under varying applied potentials [43]. The calculated I-V response for films with different thicknesses are shown in Figure 6. The leakage current density in BFO films increases markedly in samples with small thicknesses, indicating the limitations of application for ultrathin films.

Conclusions
The effects of film thickness on polarization switching paths were investigated using phase-field modeling. Results indicate that 71 degree switching plays the most important role in the polarization switching process. Analysis of time-dependent switching percentages in different films shows that polarization switching time increases as film thickness decreases. Calculations of the hysteresis loop for films with different film thicknesses reveal the dependence of a ferroelectric coercive field on film thickness, as ∝ ℎ , in ultrathin BFO films, which could potentially provide guidance in future theoretical research and material design and discovery.