Monte Carlo Simulation of Ferroelectric Behavior in PZT Films by Using a Stress Dependent DIFFOUR Hamiltonian

In this work the polarization and hysteresis response of Lead Zirconate Titanate (PZT) ferroelectric thin films was studied in relation to the variation on temperature, stress, electric field and the content of non-ferroelectric impurities by using a Monte Carlo simulation. The simulation was based on a DIFFOUR Hamiltonian that takes into account the effect of uniaxial stress, in addition to the nearest neighbor dipoles interaction and the effect of an external electric field. The obtained results for hysteresis loops and polarization curves correspond with reported experimental data for this material.

The study of variables that directly affect domain formation and the quality response of these materials under the influence of electric fields have been of wide interest, because of their potential use in the design of non-volatile memories [11] and the development of multiferroic complex perovskites from the mixture of ferroelectric and ferromagnetic simple perovskites [12]. A big amount of experimental data reports the behavior of several ferroelectric materials such as Lead Zirconate Titanate (PZT) and Bismuth Titanate (BIT) [13], considering effects of temperature, stress and dopants, among other variables in relation with polarization and hysteresis. Experimentally, Wondsamnern [14] characterized ferroelectric hystere-sis loops of PZT ceramics, observing the influence of external stress and temperature in an inverse relation of these variables with the hysteresis area. As important as those analysis that we mentioned before are research proposals trying to contribute to a better understanding of these materials through a theoretical modeling approach. Recently Laosiritaworn et al. [15] studied by Monte Carlo technique, PZT thin films by proposing a uniaxial stress dependent Hamiltonian. They simulated the dynamic response of a ferroelectric thin film, identifying changes in the energy dissipation from hysteresis loops with the stress variation.
In this research, the response of polarization, susceptibility and hysteresis of ferroelectric thin films to the variation of temperature, stress, electric field and the content of non-ferroelectric impurities was addressed through a Monte Carlo simulation. For this purpose, a DIFFOUR Hamiltonian was proposed to model ferroelectric properties in PZT thin films, one of the most studied ferroelectrics nowadays.

Simulation Methodology
PZT is a piezoelectric ceramic material in the form of a solid solution of P bZrO 3 and P bT iO 3 This material exhibits different stoichiometries according to Zr/T i composition ratio, which directly affects its properties as it is evident from its phase diagram [16]. When the Ti content is higher than 0.47, PZT transforms from cubic structure to tetragonal structure at temperatures ranging from 370 to 490 • C. When the Ti content is lower than 0.47, there is a cubic to rhombohedral transition with Curie-temperatures ranging from 230 to 370 • C [17]. The sample construction for simulation takes into account the transition expected for a thin film with a Ti content, x = 0.5, in which the perovskite geometry of the cell can be closely associated with a cubic structure (lattice parameter a ≈ c ≈ 4.07A • ) [18]. This reduces to six the number of possible electric dipole orientations for each cell (−x, +x, −y, +y, −z, z), to be consistent with 90 • or 180 • experimentally observed polarization domain switching. In this manner for the model, an electric dipole associated to a PZT cell is assigned to each point of the cubic lattice, with one of the above orientations.
A single-spin movement Monte Carlo algorithm was used for obtaining equilibrium thermodynamic properties [19]. The lattice for simulation was designed with a thin film geometry of LxLxd dimensions. Periodic boundary conditions in the transverse directions (x-y plane or in-plane) and free along the perpendicular direction z of the film (out-of-plane) were employed. A linear dimension of L = 40 and a thickness of d = 5 (lattice units) were considered. Simulations for larger sizes of L were computed, giving no significant variation for the parameters and variables under consideration.
The system began with a random electric dipole configuration, corresponding to a high disordered or a high temperature system. We used Monte carlo simulation and Metropolis Algorithm as well as a DIFFOUR Hamiltonian [20], [21], Equation (1), to simulate the ferroelectric behavior (looking for a Transition temperature or Curie's temperature) of the system, with a variation range temperature from a high temperature up to a low temperature that correspond to a high ordered system.
In this hamiltonianû i is the polarization unit vector at site i,u iz is the z component of polarization vector at site i, U ij is the magnitud of ferroelectric interaction between pairs of nearest pairs of dipoles and E is the electric field applied in the z direction. The first summation for inplane interactions, the second for out-of-plane interactions and the last one considers the electric field action over dipoles in z direction. For this model, we proposed that the stress dependent interaction intensity took the form of a Lennard-Jones potential: here, U 0 stands for the magnitud of ferroelectric interaction associated to zero strain, P is the uniaxial compressive stress in the out of plane direction. The quantities 1 + αP and 1 − βP are proportional to the augmented |68 Ingeniería y Ciencia and reduced length of distances in xy plane and z axis respectively under the action of P , counting for strain-polarization intensity relationship. That means, α and β values depend on the specific effect of mechanical stresses on the ferroelectric material and parameters to be adjusted with reported experimental data . A particular case, for α = 1/Y and β = ϵ/Y (Young modulus Y , and Poisson ratio ϵ) gives the Hamiltonian previously studied by Laosiritaworn et al. [15].
For simulations U 0 was set as 1, re-defining the unit of temperature T as (U 0 /K B ) (being K B the Boltzmann constant), the electric field E was considered in units of U 0 and stress introduced through the dimensionless ratio P/Y . The parameters used in the simulation are listed in Table 1. The updating procedure for electrical displacement at any site is done as follow: For a given temperature and starting from a random lattice site, a new random orientation is assigned to u i among the six possible directions, then the Hamiltonian is calculated and the Metropolis algorithm is used to update the electric displacement in the site [22]. This procedure is repeated for all sites in the lattice(One Monte Carlo step (MCS)), and repeated as many times as necessary to reach the relaxation of the system. We used 10.000 Monte Carlo steps (MCS) for the relaxation of the system and 5.000 MCS more, to determine the average for polarization along the z direction and susceptibility, Equations (4) and (5); N is the total number of dipoles in the sample. The same algorithm was used to obtain the hysteresis loops, but in this case the electric field was varied instead of the temperature.

Results and Discussion
The phase transition from paraelectric to ferroelectric and the hysteresis phenomena were studied. First by varying the adimentional ratio (P/Y ) 0 to 0.008, secondly considering an electrical field variation in the out-ofplane direction from 0 to 2.0 U 0 units, finally considering non-ferroelectric impurities (sites in the lattice without dipole) from 0 to 100% concentration.
The susceptibility curves as a function of temperature are shown for various stress values in Figure 1. A decrease in intensity and critical temperature with increase of stress can be observed.This could be explained considering that in the Hamiltonian, the stress in the out-of-plane direction favored the dipole flipping into the in-plane, becoming z a hard axis of polarization and then reducing the formation of domains in the out-of-plane direction. The LGD (Landau-Ginzburg-Devonshire) theory estimates that for this kind of ferroelectric, a 100 MPa stress will shift the transition temperature by about 20 • C. In this matter, a systematic study of the ferroelectric properties of PZT nanowires was done by Hong and Fang [23] considering this thermodynamic approach. Two competing mechanisms had been proposed to explain this transition temperature shift. One, the uniaxial stress causes a distance diminution between dipoles and consequently an increase in ferroelectric coupling intensity leading to higher transition temperatures. On the other hand, uniaxial stress causes a movement of the Zr/Ti cation to a more symmetric position in the cell, reducing the dipole magnitude and the ferroelectric coupling in the out of plane direction, favoring in this way, the shift to lower transition temperatures. This is closely related to the variation of hysteresis loops with the applied stress, Figure 2. A reduction in their area with the increase of stress, and thereby, a decrease in the coercive field can be observed due to non easy alignment in z direction, Table  2. Experimental data on ferroelectrics and in particular on PZT systems coincide reasonable well with this simulated behavior [24], [25], [26], [27].  On the other hand, the observed effect of an increasing field in the out-of-plane direction is the horizontal shift of the polarization curves, accompanied with higher values of polarization and minor deviating slopes, meaning a less defined critical temperature peaks for the corresponding susceptibility curves, Figure 3. The presence of the electrical field E term in the Hamiltonian Equation (1) Figure 4 shows hysteresis loops at different non-ferroelectric impurity percentages from 0 to 40%. The increased presence of these impurities in the lattice, causes a transition temperature reduction and the decrease on the hysteresis loop areas Table 2. This is due to the fact that less states of polarization for the system are available with the impurities. Therefore, less dipole-dipole interaction and domain formation results. The effect of these defects in the sample is notorious, becoming an important parameter to control these ferroelectric properties.
At higher percentages of non-ferroelectric impurities the simulated sample presented critical temperatures very close to absolute zero and insignificant hysteresis areas, meaning a paraelectric behavior. In an experimental way, the hysteresis area can be changed with the influence of dopants on the electrical performance of these materials, including imprint effects due to pinned dipoles and material fatigue [28].

Conclusions
Monte Carlo simulation of the ferroelectric phenomenon on PZT thin films was carried out. The effect of outstanding variables on the ferroelectric performance as the stress, electric field and content of non-ferroelectric impurities were evaluated by proposing a DIFFOUR Hamiltonian with a stress dependent ferroelectric coupling coefficient. The results showed a significant influence of these variables in the hysteresis loops and polarization curves. The stress effect in ferroelectric sample, leads to a change in the probability of dipole switching along its axis of application. An increase in stress in the out-of-plane direction diminishes the number of domains in this axis, favoring dipole switching to the in-plane direction. All of the above means, that rising the stress causes lower values in the susceptibility peaks and a shift of them toward lower temperature values. Besides, this out-of-plane domain diminution (also caused by non-ferroelectric impurities) yields to a reduction in the hysteresis area. These effects correspond with the qualitative behavior in the experimental reports observed.