Depolarizing field in temperature-graded ferroelectrics from an atomistic viewpoint

An atomistic approach for computing the depolarizing, or internal, electric field in materials with inhomogeneous polarization is developed. Application of the approach for studying the depolarizing fields in technologically important (Ba0.7Sr0.3)TiO3 ferroelectric alloy with temperature gradients has revealed the intrinsic features of these fields as well as their role in the establishment of polarization gradients. It is found that the depolarizing fields are inhomogeneous and can be tailored to yield both zero and non-vanishing potential drops. Such findings could pave the way to unusual thermoelectric materials, photovoltaics and locally conducting materials, all of which are at the frontier of current research.

From an electrostatic point of view, the variation of polarization leads to an appearance of an electric field, which is often referred to as the depolarizing, or internal, field. The depolarizing field in ferroelectrics is typically quite large [36] owing to the large values of spontaneous polarization [34]. As a result, its effect on the polarization is dramatic and believed to cause the polarization offsets in polarization-graded ferroelectrics [43], unusual domain patterns in nanoscale ferroelectrics [29] and enhanced photocurrent [15,42] in ferroelectric thin films. In the presence of any free charge this depolarizing field may be reduced or even completely annihilated, leading to spatially inhomogeneous free charge distribution and even local conductivity. Owing to its contributions to the properties of ferroelectrics, the depolarizing field emerges as a key milestone on the way toward understanding and tailoring the properties of ferroelectrics and generally dielectric materials. Interestingly, despite its key role, the depolarizing field remains elusive for computations in many cases of practical importance. One of them is the depolarizing field associated with gradient of polarization.
Consider a polarization-graded ferroelectric sample in which the polarization changes along the z-direction. The requirement of the electric displacement field D continuity dictates that along the sample. Here E(z) is the internal, or depolarizing, field that exists in the sample, while P(z) is the polarization. Differentiating equation (1) with respect to the position z within the sample allows the determination of the electric field gradient along the sample Integration of this equation yields the electric field along the sample where E(z 0 ) is the reference field. In this paper, we will consider the case of a sample sandwiched between two electrodes that provide a compensating surface charge of where z 0 and z 2 indicate the positions of the left and right interface with the electrodes, respectively.
There are several challenges in using equation (3) to determine the internal field in the sample with a polarization gradient. The first one is that the integration in equation (3) requires a rather precise knowledge of P(z) under the condition of non-zero electric field E(z). Experimentally, this is quite difficult to achieve. Computationally, at the atomistic scale the polarization profiles may exhibit some rapidly varying features which could complicate evaluation of the field from equation (3). The second challenge is that a knowledge of the reference field E(z 0 ) is required. As a result, while the internal field in polarization-graded material represents a clear physical concept, its practical determination either from experiment or computations remains challenging.
Alternatively, we can start from the relation P(z) = P sp (z, E = 0) + ε 0 (ε(z) − 1)E and rewrite equation (1) in a different form where ε(z) is the dielectric constant of the material that may vary along the sample. P sp (z, E = 0) is the spontaneous polarization in the material under the condition of zero electric field. Repeating the above derivation, we obtain the following expression for the internal field: However, this equation does not offer much advantage compared to equation (3), since one needs to know the spontaneous polarization profile in the sample under the condition of zero electric field, which is a rather artificial condition for many cases of polarization-graded ferroelectrics. In addition, a dielectric constant along the sample must be known. As in the case of equation (3) one needs to determine the reference field E(z 0 ) as well.
Therefore, one would like to have a computational approach that would allow us to bypass the use of equation (3) or (6) and obtain the depolarizing field directly from atomistic simulations. In this paper, we (i) develop such an atomistic approach to calculate depolarizing field due to polarization gradients; (ii) use the approach to compute the depolarizing fields associated with the existence of temperature gradient in technologically important ferroelectric (Ba 0.7 Sr 0.3 )TiO 3 alloy; and (iii) reveal the intrinsic features of such depolarizing fields as well as some possible ways to tailor them. More specifically, we will show that the depolarizing fields due to polarization gradients could be tailored to yield both zero as well as the finite potential drops. The latter may have important technological consequences. The proposed computational approach can be generalized to other cases of inhomogeneous polarization fields.

Methods
In this paper, we study the depolarizing fields in temperature-graded (Ba 0.7 Sr 0.3 )TiO 3 alloy.
(Ba x Sr 1−x )TiO 3 alloys are very attractive for their exceptional compositional tunability. For example, their Curie temperature varies almost linearly with Ba concentration ranging from 0 up to 405 K [44]. Our choice of (Ba x Sr 1−x )TiO 3 alloy with x = 0.7 is motivated by the facts that (i) its Curie point is very close to the room temperature; (ii) it exhibits a rich phase diagram in a relatively narrow range of temperatures; and (iii) some experimental data are available for the temperature-graded alloy with x = 0.75 [1]. The alloy is modeled by a 12 × 12 × 31 supercell (22 320 atoms) which is periodic along all three directions to simulate the bulk system. We chose x-, yand z-axis to lie along the (100), (010) and (001) crystallographic directions, respectively. The total energy of this supercell is given by the first-principles-based effective Hamiltonian [45] that depends on the local soft mode (which is proportional to the electric dipole moment), inhomogeneous and homogeneous strain variables, and alloy configuration. It includes a local mode self-energy (harmonic and unharmonic contributions), a long-range dipole-dipole interaction, a short-range interaction between local modes, an elastic energy, the interaction between the local modes and strains, and the interactions responsible for alloying effects. This Hamiltonian correctly reproduces the complex sequence of phase transitions in (Ba x Sr 1−x )TiO 3 alloys for a wide compositional range 2 and has been used to study many properties of ferroelectric alloys yielding results in good agreement with both experiment and first-principles calculations (see [45,47,48]).
To simulate the temperature gradient, we use the direct method within non-equilibrium molecular dynamics (MD) [16,49]. The supercell used in the calculation is divided into N = 31 slabs along the direction of the measurement (z-axis in our case). This supercell is schematically shown in figure 1(b). One of these slabs is used as a heat source, whereas another is used as a heat sink [16]. For every MD step, each particle velocity in the source/sink region is rescaled to generate a heat flux along the z-direction. The system is then allowed to reach steady state with the resultant temperature gradient controlled by the magnitude of the heat flux. We begin our simulations by equilibrating the system at the desired temperature T eq using 500 000 MD steps within NPT ensemble. The heat flux is then turned on and the simulations are continued for another 4 000 000 MD steps using non-equilibrium MD [16]. Last 1 000 000 MD steps are used for computing averaged values of polarization, temperature and other properties.
Once the temperature gradient is established, we proceed with the calculation of the depolarizing field as follows. The local field at a point r of the sample is where E app (r) is the applied electric field induced by external sources, while E dip (r) is the field produced by all the dipoles inside a sample or simulation supercell. On the other hand, we can write the local field using the macroscopic fields as Here E dep (r) is the depolarizing field and E Lorentz (r) = P(r)/3 ∞ 0 is the Lorentz field due to the charge density P on the surface of the Lorentz sphere. Note that ∞ in the denominator reflects the fact that the field is produced in the medium with a dielectric constant of ∞ . Here we assume that outside the Lorentz sphere the polarization field is homogeneous and is equal to P(r). E sph (r) is the field produced at the point r by all the dipoles inside the Lorentz sphere. Because of the cubic symmetry of the perovskite crystal E sph (r) = 0. Combining equations (7) and (8), we obtain This equation relates the atomistic electric field E dip (r) due to internal dipoles to the macroscopic electric fields. Note that in the absence of an applied electric field E dip (r) = E loc (r) and one recovers the familiar relationship between the local and depolarizing fields in longitudinal optical modes of the ionic crystals [50].
To derive the expression for the depolarizing field due to polarization gradient, we write equation (9) for two cases: with and without polarization gradient, referred to as PG and NPG cases, respectively, In the first one of equations (10), we took into account the fact that when no polarization gradient exists in an infinite sample, the depolarizing field is zero. From the last two equations the depolarizing field can be computed as It is convenient to rewrite the equation in a short form where indicates the difference in the fields between the PG and NPG cases. Equivalently, we can rewrite the latter equation as Equation (13) allows one to calculate the depolarizing field due to the polarization gradient from the difference in the dipole and polarization fields which are associated with nonequilibrium (PG) and equilibrium (NPG) conditions. Note that equation (13) is general and applies to other cases associated with the existence of depolarizing field. For example, it was proposed in [51] that the depolarizing field in ferroelectric nanostructures can be calculated using the difference in the dipole field E dip (r) associated with different boundary conditions. It is straightforward to see that this depolarizing field is a special case of equation (13) when P(r) = 0. It will be instructive to verify that the depolarizing field obtained using such an atomistic approach satisfies the condition of continuity of the electric displacement field in the supercell. To test this condition, it is convenient to compute the difference in the electric displacement fields D between the PG and NPG cases. Since D must be constant along the z-direction (direction of polarization) for both the NPG and PG cases then the corresponding difference D z must also remain constant along the sample. Such a difference is given by Note that for the cases considered, the depolarizing field exists only along the z-direction and is denoted as E dep .

Results of the simulations
We first carry out simulations in the tetragonal state of the alloy near equilibrium temperature of T eq = 245 K which is 40 K below the Curie point. The supercell is polarized along the positive z-direction. The temperature gradient established along the z-direction leads to the associated polarization gradient along the same direction (longitudinal case). All the fields are calculated at the lattice sites r and then averaged in the (x, y) planes. The use of our computational approach (equations (7)-(13)) to study ferroelectric phases requires some clarification. The approach was derived under the assumption of a cubic symmetry. While ferroelectric phases are generally associated with non-cubic symmetries, the actual distortions from the cubic structures are very small [52]. As a result, the contribution of these distortions to the electrostatic energy is of higher order and, therefore, omitted in the effective Hamiltonian [52]. In other words, the dipole-dipole interaction term of the Hamiltonian 'sees' the cubic structure. This justifies the use of equations (7)-(13) for the depolarizing field calculations in ferroelectric phases.  Figure 2 shows the computational data obtained for some temperature gradients established near T eq . The left column gives the data for the smallest temperature gradient considered. The temperature changes from 240 K at the point of the heat sink to 250 K at the point of the heat source (see figure 2(a)). The entire temperature range is associated with the tetragonal phase of the alloy under equilibrium conditions. Figure 2(d) reports the depolarizing field profile along the supercell as obtained using the atomistic approach described above. The atomistic simulations confirm the existence of a depolarizing field due to the polarization gradient as predicted by the macroscopic equation (3). Under given conditions the field varies linearly along the supercell changing its sign from positive to negative upon crossing the boundary between warm and the cool parts of the supercell. Such boundaries are shown by dashed vertical lines in figure 2. The depolarizing field profile indicates that in the warm area of the sample, the depolarizing field acts to support the polarization (and should probably be called 'polarizing' rather than depolarizing), while in the cool part of the sample the field acts against the polarization. Such a behavior is a consequence of the competition between the internal forces acting to maintain the spontaneous polarization associated with a given temperature profile and electrostatic forces acting to conserve the electric displacement field along the supercell. The magnitude of the depolarizing field reaches up to 2 MV m −1 in the hottest and coldest areas of the supercell. The depolarizing field is zero at the boundaries between the warm and cool parts of the supercell. Interestingly, if free charge is present in the sample (such as extrinsic carriers, e.g.) it may be driven along the field lines toward the point of zero electric field which will act as a carrier sink or a source. This could lead to local conducting channel(s) which are potentially very attractive for many technological applications, including non-destructive ferroelectric memory [17,18]. Under the given conditions the average depolarizing field inside the sample vanishes. Figure 2(g) shows the change in the z component of polarization P z associated with the establishment of polarization gradient. We define the polarization offset with respect to the equilibrium profile as the average P z along the supercell. Such an offset is well approximated by the P z at the boundary between the warm and cool regions (dashed vertical line). The polarization offset is negative and indicates that the depolarizing field is minimized through a reduction in the average magnitude of polarization.
We have used equation (14) to calculate the change in the electric displacement field D z inside the supercell associated with the establishment of steady-state non-equilibrium profiles given in figure 2. The data are shown in black circles in figure 2(g) and demonstrate that D z remains constant along the supercell as dictated by the electrostatics.
In the second and third columns of figure 2, we report the data obtained for larger temperature gradients. As T increases the depolarizing field also increases in magnitude, while maintaining its linear profile with the zero average field value. We have compared the data for the depolarizing fields obtained from our atomistic approach with the depolarizing fields predicted by equation (3) and shown by solid black lines in figures 2(d)-(f). It is found that once the reference field is known the depolarizing field from equation (3) fits the fields given in figure 2 with the accuracy of 4%. The polarization offset increases as a function of T reaching 150 nC cm −2 in magnitude for T = 100 K.
Note that for all the cases reported in figure 2, the alloy remains in the tetragonal phase, even when the temperatures of the hottest and coldest regions of the supercell are outside of the temperature range associated with the tetragonal phase of the alloy under equilibrium conditions. Interestingly, our data suggest that under the given conditions switching the positions of the heat source and sink should not change the sign of the polarization offset since the overall effect of the depolarizing field is to reduce the polarization magnitude. Indeed, this was confirmed in additional calculations where the positions of the heat source and sink were reversed. Furthermore, it appears that under such conditions the depolarizing field will not offset the hysteresis loops since it acts to reduce the polarization magnitude and will, therefore, produce symmetric shifts of polarization. We expect that under these conditions, the loops should shrink along the D-axis but remain symmetric with respect to the origin. To check this, we reversed polarization direction in the supercell while keeping the heat sink and source in their places. Indeed, we found that P z changed its sign suggesting symmetric offsets.
Next, we apply a much larger temperature gradients to induce formation of multiple phases in the supercell. The data from such simulations are given in figure 3 and the data indicate that in this case the temperatures inside the supercell sample the entire family of phases associated with the (Ba 0.7 Sr 0.3 )TiO 3 alloy under equilibrium conditions. There are two notable features that differentiate these data from the data obtained for smaller temperature gradients. Firstly,  we indeed observe the coexistence of multiple phases as previously reported [16]. The lowsymmetry monoclinic phases are formed in the cold region of the supercell (shaded region in figure 3). Secondly, we observe that the average depolarizing field no longer vanishes, but exhibits a finite value which increases with T (see figures 3(d)-(f)). In fact, for the smallest temperature gradient reported in figure 3(d) the average depolarizing field reaches 5 MV m −1 in magnitude. This non-zero average depolarizing field, and, therefore, non-vanishing potential drop, is a consequence of the phase 'inhomogeneity' associated with the given conditions. Such 'inhomogeneity' can be described as a rotation of polarization away from the z-direction in the cool part of the supercell. As a result the z-component of the polarization experiences further decrease which is partially absorbed by the increase of the depolarizing field. Interestingly, if free charge is present in the sample the depolarizing field will then drive the positive charge from the hot part of the sample toward the cold part of the sample. This potentially could offer an unusual approach to thermoelectrics. Moreover, the ability of a material to produce such an internal field is very advantageous for photovoltaic applications, where the electric field is required to separate photocarriers and decrease their recombination rates. Figures 3(g)-(i) give the polarization offsets that reach very large values under such conditions.
We have carefully examined the dipole patterns associated with different temperature gradients. We found that for the smallest gradients considered (see figure 2), the simulated sample remains in a monodomain state. For largest temperature gradients (see figure 3), on the other hand, the local dipole vector gradually changes from (0, 0, v) (near heat source) to (u, u, v) (near heat sink). This dipole pattern could, in principle, be classified as a polydomain state. However, the polarization vector rotates smoothly, so no domain boundaries or walls could be defined.

Conclusions
In summary, we have developed and used an atomistic approach to compute depolarizing fields associated with inhomogeneous polarization fields. Application of this approach to study temperature-induced polarization gradients in technologically important (Ba 0.7 Sr 0.3 )TiO 3 alloy revealed some intrinsic features associated with these depolarizing fields. For the temperature gradients considered, we found that the depolarizing field changes linearly inside the simulated sample and can reach rather large magnitudes. Moreover, the field changes its sign in such a way that it acts to support polarization in the warm areas of the sample and to reduce polarization in the cold areas of the sample. This change in the sign of the depolarizing field could be utilized to create locally conducting channels. We also found that for the largest temperature gradients considered, the simulated sample develops monoclinic phases which lead to the appearance of a non-zero average depolarizing field and, therefore, a non-zero potential drop along the sample. Our results suggest that the depolarizing fields in the temperature-graded ferroelectrics can be tailored via manipulation of the temperature profiles. These could find potential applications in unusual thermoelectrics, locally conducting materials as well as in photovoltaics.