Investigation of domain structure in ferroelectric thin films by means of the Ising model

The domain structure in thin ferroelectric films was investigated in the framework of the modified Ising model with short-range exchange and long-range dipole-dipole interactions. The depolarizing field was also taken into account for the study of size effects in thin films. This field is the dominant factor influencing on the existence of phase transitions in the films. We used Monte Carlo method with the standard Metropolis algorithm to generate configurations on three-dimensional cubic lattices. The domains sizes were calculated in the longitudinal and transverse directions relatively to the spontaneous polarization direction. Temperature dependences of the domain sizes are obtained for different values of interaction constants. It was proved that the introduction of the depolarizing field leads to the appearance of a dead layer on the film boundaries. The dead layer thickness increases with increasing temperature. At a certain temperature (not equal to the Curie temperature), the dead layer is destroyed. The temperature dependences of the dielectric susceptibility at different film thicknesses have been calculated. It was shown that these dependencies have two maxima. The first maximum corresponds to a phase transition, the second one exists only for only thin films and corresponds to the destruction of the dead layer.


Introduction
One of the most important theoretical tools for studying phase transitions is the Ising model. It allows to explore different properties of many magnetic [1] and non-magnetic systems [2] on standard lattices (square, triangular, hexagonal) and more complex ones (with geometric frustration, Bethe, recursive); with different spins at the sites of lattices and interactions between them. The classical model predicts only well the behavior of uniaxial bulk ferroelectrics [3]. However, it is not applicable for description of thin films, because cannot predicts the existence of a critical film thickness at which a phase transition does not exist. It has been experimentally found that a decrease in film thickness or grain size leads to a decrease in dielectric constant, remanent polarization, dielectric breakdown field, an increase in loss tangent, coercive field, band-gap energy, and diffuseness of the phase transitions [4,5]. In theory, however, there is no consensus on the specific mechanisms of the size effects, but it is known that the influence is exerted by electrical and mechanical conditions at the substrate-film boundary, incomplete screening of the depolarizing field, and impurities [6]. The study of thin films is relevant, since they are an integral part of storage devices, microwave electronic components, microsensors, actuators etc. Most of the important applications of ferroelectric materials are largely  [7]. The aim of this work is to study the influence of size effects on changes in the domain structure of a ferroelectric thin film during the phase transition.

Model
In contrast to the classical Ising model, we took into account the energy of long-range dipole-dipole interactions and the energy of dipoles in a depolarizing field, since bound charges in real ferroelectrics create an electric field that tends to change the direction of spontaneous polarization near the surface. We suppose that the magnitude of the depolarizing field is proportional to the long-range orientational order parameter that is determined by the formula: are the numbers of lattice nodes in the direction of the x,y,z axes of the Cartesian coordinate system, respectively. We suppose that the z axis is directed perpendicular to the free surface of the film. The quantity i S takes values 1 ± . The configuration energy of the dipoles is calculated by the following formula: where the constants J , D characterize exchange and dipole-dipole interactions, respectively, the constants E and λ are determined by the number of free charges in the film and its electrical conductivity, the value d is the lattice constant, the quantity ij r characterizes the distance between i-th and j-th dipoles, ij z is the z-axis projection of the distance ij r , i z is the z coordinate of the i-th dipole, The first term in equation (1) characterizes the exchange interaction, this sum is calculated over the nearest neighbors of the i-th dipole. The second term in equation (1) defines dipole-dipole interactions, the sum is calculated over the neighbors located inside the ball radius with the i-th dipole at its center. The third term in equation (1) is calculated over all lattice dipoles and it is equal to the potential energy of the dipoles in the depolarizing field. As in the work [8], we assumed that the depolarizing field inside lattice decreases exponentially from both free surfaces of the film due to the compensation by different screening mechanisms.
The calculations were performed by the Metropolis algorithm for periodic boundary conditions in the direction of the x and y axes. Free boundary conditions were taken along the z-axis direction, that is, the influence of electrodes was not considered.  anisotropy of the domain sizes increases, i.e. the value Δ increases. This behavior is explained by the difference in the spin configurations with the minimum energy when the ratio D J changes [10]. As the inverse temperature decreases, the point at which anisotropy begins to grow shifts to the right.

Results and Discussion
The temperature dependences of the susceptibility χ of a ferroelectric thin film at different values of its thickness are shown in figure 2. The quantity χ was calculated by the formula: where P is the maximal value of the film polarization equal to the value H N N 2 1 μ . The depolarizing field strongly affects the convergence of the simulation results, especially nearly the phase transition point. In real experiments, the heating or cooling of samples occurs at a finite rate. Therefore, we could use the simulation with not very large number of Monte Carlo steps, namely, 10 steps per dipole. To obtain satisfactory results, it is necessary perform the averaging is over 300 steps.
As can be seen from the figure 2, the dependence χ on β has two peaks. The large peak corresponds to the Curie point. The small peak has also been found experimentally [11]. It is usually assumed that it arises from the inhomogeneity of the film material. But theoretical studies of bulk solid solutions have shown that, there is one susceptibility peak, which is shifted depending on the impurity concentration. To eliminate this contradiction, we investigated the domain structure of the film.
Calculations have shown that the introduction of the depolarizing field leads to the formation of a layer near the film boundaries, where the polarization of the dipoles is opposite to the spontaneous polarization. Such a layer is called dead one and its existence in thin films was discussed in a number of theoretical [8] and experimental [12] works. Figure 3 shows the dependences of the ratio of the average dead layer size h along the z axis, to the film height H. With increasing temperature, the value of H h/ increases due to the growth of the near-surface layer. At a certain temperature, the domain structure is rearranged and the dead layer thickness reaches its maximum value. In this case, a diffuse structural transition occurs and a jump of the susceptibility is observed. The magnitude of this jump depends on the film thickness and is not observed in bulk samples.
With a further increase in the temperature, the dead layer is destroyed, so the value H h/ decreases. At the Curie point, the ratio H h/ also experiences a maximum due to the fluctuations, which grow infinitely at the phase transition point.

Conclusion
An analysis of the domain structure of the thin ferroelectric films showed that the temperature dependence of the susceptibility has an additional jump that is not associated with a change in the crystal structure of the substance. It is shown, that at this point ( d β ), there is a transition from a single-domain structure to a multi-domain one. A more detailed study of the system is required, since in this case d β β < additional depolarizing fields arise in the film, the magnitude of which depends on the domains size.