Polarization Switching in Ferroelectric Thin Films Undergoing First-Order Phase Transitions

The main switching properties in ferroelectrics undergoing first-order phase transitions are simulated within the framework of the extended Ishibashi dipole-lattice model including the dipole-dipole interaction in a two-dimensional case for ferroelectric nanoscale objects. The peculiarities of the temperature dependence of the switching rate and the pyroelectric coefficient are discussed in the range of coexistence of the metastable states. The used coefficients of the long-range and short-range interactions between the dipoles are taken from the dielectric and structure measurements in BaTiO3.


Introduction
Nanoscience and nanotechnology is one of the important frontiers in scientific research.The scientific community is addressing a wide range of research topics from fundamental physical, biological, and chemical phenomena to materials science at the nanoscale [1][2][3].A particularly important field of fundamental and applied research is the science of ferroelectric thin films.Ferroelectric thin films have a wide range of properties which are used for applications in nonvolatile ferroelectric random access memories, microelectronic mechanical systems, nonlinear optical devices, and sensing applications [1,2].Since many ferroelectric substances, including nanometer thin films [4], reveal firstorder phase transitions, there is a range of metastability in which the paraelectric and ferroelectric phases coexist.The kinetics of polarization switching in the region of the existence of first-order phase transitions characterized by the appearance of metastable states is of great interest.Usually switching properties are studied without examining the special features appearing in the metastable state.
They proposed simulation of ferroelectric properties using a one-dimensional lattice of dipoles.The discretization enables to estimate correctly the relevant physical properties of thin films.In most of these publications, the interaction between dipoles is represented by the interaction between two neighboring atoms.This means that Ishibashi model takes into account only the short-range interaction and ignores the long-range dipole-dipole interactions which are needed for understanding of the ferroelectric instability.Simplified estimations for the one-dimensional case taking into account the dipole-dipole interactions are published in [13,14].We have recently extended Ishibashi model including the dipole-dipole interactions which are necessary for the description of the ferroelectric ordering for secondorder phase transitions in the two-dimensional case [15].Since the extended Ishibashi model modified by us includes the long-range and short-range interactions between dipoles in the two-dimensional case, it is more realistic for the description of the polarization switching in real ferroelectric films than the known approach.In this communication, we simulate switching processes in the range of the existence of metastable states.We calculate the rate of polarization switching, switching current, and pyroelectric coefficient in the range of the phase coexistence.In contrast to the previous simulations made in the framework of Ishibashi dipole lattice model [5][6][7][8][9][10][11][12][13][14], we use parameters taken from dielectric and structure measurements in BaTiO 3 .Thus, we essentially improve the procedure of simulation carried out so far in the framework of Ishibashi dipole lattice model.

Model
We examine the ferroelectric system, which is the system of dipoles in the nodes (x i , y j , z k ) of a three-dimension network, where each lattice site has its own dipole moment p i, j,k and all dipoles are oriented along the X-axis.For the system undergoing a first-order phase transition, the free energy has the following form: Here E is an applied electric field, oriented along X, α > 0, β < 0, and η > 0 are coefficients of the Landau-type expansion, κ is the coefficient of the interaction between neighbouring dipoles, and Δ is the distance between two neighbouring dipoles.
The dynamics of each dipole moment is described by Landau-Khalatnikov equation [1] where γ is the kinetic coefficient.
We will consider the case when the polarization depends on two indices only, that is, p = p(x i , y j , t) = p i, j (t).We will name such a system two-dimensional one.In this case, (2) is ( Here ).Let us introduce dimensionless variables: the polarization will be measured by using the polarization unit p un = Δ 3 −β/η; the time is given in units γη/β 2 .The division of the equation by β 2 /η leads to Here ν = κη/β 2 Δ 2 is the coefficient of the interaction between nearest neighbours, D = η/β 2 is the coefficient of the dipole-dipole interaction, and E = (η 3/2 /|β| 5/2 )E is the dimensionless external field.Coefficient a is determined by the relation where T * is the stability limit of the ferroelectric phase and T 0 is the stability limit of the paraelectric phase in zero electric field.It is convenient to introduce the dimensionless temperature δ = (T − T 0 )/(T c − T 0 ), where T c is phase transition temperature.Then the coefficient a can be rewritten as a = (3/16)δ.The global polarization is the average of dipole moments over all the sites The switching current can be written as Note that the equations describing the polarization evolution for the system undergoing a second-order phase transition can be obtained by inserting the expression ap i, j − p 3 i, j − p 5 i, j instead of the first three terms in right-hand part of (4).In this case, coefficient a has a form a = αη/β 2 = (α (T c − T)η)/β 2 = a (T c − T), where T c is the phase transition temperature.For the systems undergoing secondorder phase transitions, one can introduce the dimensionless temperature θ, so as coefficient a would have the same form as in the case of first-order phase transitions.Really, if we choose the temperature T 1 according to formulae T 1 = T c − (3/16a ) and determine the dimensionless temperature θ = (T c −T)/(T c −T 1 ), then a can be rewritten as a = (3/16)θ.
In the homogeneous case in zero external electric field, equations for all p i, j will coincide with the equation for the global polarization in the absence of the dipole-dipole interaction and external electric field The maximum of absolute values of roots of the equation g(P) = 0 gives the following polarization value in the ferroelectric phase: The value of the maximum of g(P) gives the value of the coercive field E c .The positions of the maxima are ±P 1,2 = ± 0.3 ± √ 0.09 − 0.2a, and the dimensionless coercive field If the external electric field is larger than the coercive field E c , then there is only one global minimum of the free energy, and local minima do not exist.We consider the evolution of the ferroelectric twodimensional system described by (4) in square area, which consists of N × N dipoles.We suggest that the solution is periodical in X and Y directions with a period N and use the periodical boundary conditions.The numerical solution of ( 4) is based on exploiting the discrete Fourier transform.The integration of the equations for Fourier components is performed by explicit Euler technique.The similar approach was used in [16,17] for the simulation of the domain structure evolution in ferroelectric systems undergoing second-order phase transitions.

Results and Discussion
First of all, we investigate the transitions between metastable and stable states in ferroelectric systems undergoing firstorder phase transitions in the two-dimensional (2D) case.The calculations are carried out in two temperature ranges: (a) the two minima corresponding to the ferroelectric phase in two domains (metastable domains) are above the minimum corresponding to the paraelectric phase at zero electric field, δ = 1.330;(b) the two minima corresponding to the ferroelectric phase in two domains (stable domains) are below the minimum corresponding to the paraelectric phase at zero electric field, δ = 0.7407.The number of dipoles is N × N = 64 × 64 and N × N = 128 × 128, that is we consider the systems of a nanometer size.Since, experimentally, barium titanate (BaTiO 3 ) is the most intensively studied ferroelectric material, the parameters of our model are chosen in correspondence with those of ) and from the first stable ferroelectric state to the second stable ferroelectric state (for δ = 0.7407, curve 2) for the 2D system as a function of the dimensionless applied electric field E/E c ; E c is the coercive field.The number of dipoles is N = 64 × 64.
The initial state corresponds to the ferroelectric state with a nucleus of the second ferroelectric domain with the dipoles oriented in the opposite direction (the size of the nucleus is 4 × 4).Parameters of the calculation are chosen in accordance with the data for BaTiO 3 : the coefficient of the dipole-dipole interaction is D = 200; the coefficient of the interaction between nearest neighbours is ν = 40.These coefficients are taken from the dielectric and structure measurements in BaTiO 3 [17].
BaTiO 3 .In particular, we choose β = 10.8 • 10 −13 CGSE, η = 2.28 • 10 −22 CGSE [17], and κ = λ 2 • π/15, where λ = 4 • 10 −8 cm is the lattice constant [18].We assume that the distance between two neighbouring dipoles Δ coincides with the lattice constant.Then for dimensionless variables D and ν, we have D ≈ 200, ν ≈ 40.In all calculations of the polarization switching and polarization reversal, we use the initial state describing the situation when the ferroelectric system is in one phase and there is the nucleus of another phase inside.The size of the nucleus is usually equal to 4 ×4.
One can estimate the influence of the nucleus size on the switching time.It is possible to show [14] that the switching time at large D and E < E c is mainly determined by the velocity of the one-dimensional motion of a domain wall, V .Under these conditions, the switching time is approximately proportional to the width of the area, where the polarization is close to −P * .Taking into account that this width is equal to N − n, where n is the size of the nucleus, one can conclude that τ ∝ (N − n)/V = τ 0 − δn.Here τ 0 is the switching time for the system with negligible small nucleus.Thus, the increase in the nucleus size reduces the switching time by the value proportional to the size.
Figure 1 represents the switching rate, 1/τ, from the first metastable state to the second metastable state (case (a)) and from the first ferroelectric stable state to the second Figure 2: Switching rate, 1/τ (in units of β 2 /(γη)), from the metastable state to the stable state for the 2D system as a function of the dimensionless applied electric field E/E c ; E c is the coercive field.The number of dipoles is N = 64 × 64.The initial state corresponds to the ferroelectric state with a nucleus of the domain of the paraelectric phase (the size of the nucleus is 4 × 4) for δ = 1.330 (curve 1) and vice versa for δ = 0.7407 (curve 2).The calculations are carried out using D = 200, ν = 40.These coefficients are taken from the dielectric and structure measurements in BaTiO 3 [17].
Figure 4: Dependence of the switching current j on time t (in units of γη/β 2 ) during the polarization reversal process in the systems undergoing first-order and second-order phase transitions (transition from the first stable state to the second one).The values of dimensionless temperature for the system undergoing a firstorder phase transition are δ = 0.9835 (curve 1) and δ = 0.7407 (curve 2).The value of θ for the system undergoing a second-order phase transition is θ = 1.6 (curve 3).The number of dipoles is 64 × 64.The initial state corresponds to the ferroelectric state with a nucleus of the second domain of polarization of the opposite polarization orientation (the size of the nucleus is 4 × 4).The value of the applied electric field is Figure 6: Dependence of the switching rate, 1/τ (in units of β 2 /γη), from the first stable ferroelectric state to the second one (or from the first metastable state to the second one) for the 2D system undergoing a first-order phase transition at E = E c (δ = 1.17) on the dimensionless temperature δ.The number of dipoles N = 64 × 64, ν = 0.1 and D = 0.05 (curve 1), D = 0.5 (curve 2), D = 5 (curve 3).The initial conditions describe the nucleus of the ferroelectric state inside the domain with polarization of the opposite direction.One can see that the behavior of the dependence is nonmonotonical.
depends on the number of dipoles: the larger the number of dipoles, the smaller the switching rate.Figure 2 demonstrates the dependence of the switching rate, 1/τ, from a metastable ferroelectric state to the second ferroelectric stable state for cases (a) and (b) on the applied electric field E/E c under the same initial conditions.The different signs of the values of the switching rate for δ = 0.7407 and δ = 1.330 reflect the fact that in the first case the transition occurs from the paraelectric phase to the ferroelectric phase, while at δ = 1.330-in the opposite direction.
A property which is especially sensitive to the presence of metastable phases and phase coexistence is switching current.We investigate the dynamics of the switching current j at the temperature range in the vicinity of phase transitions within the framework of the 2D model by using the system of (4).The calculations were carried out for δ = 1.017 and δ = 0.984.At δ = 1.017, the two minima corresponding to the ferroelectric metastable domains are above the minima corresponding to the paraelectric phase at zero electric field, while at δ = 0.984 the situation is reversed.The number of dipoles is 64 × 64.The initial state corresponds to the nucleus of a ferroelectric domain inside the ferroelectric phase of the second domain; the size of the nucleus is 4 × 4. The dependences of the switching current on time at D = 200, ν = 40, E = 0.95 • E c , and E = 1.05 • E c , are given in Figure 3.One can see that the dependence of the form of j(t) on temperature (in the vicinity of the phase transition temperature T c ) and E (in the vicinity of E c ) is rather weak.We investigate the influence of the coefficient of the interaction between the closest neighbours ν on the transition rate.The calculations show that the increasing of  There is a singularity in the vicinity of θ = 0 of the order of 1/2.The pyroelectric coefficient for systems with the first-order phase transition is described by curve 2. There is a discontinuity of the dependence (dP * /dδ)(δ) in the vicinity of δ = 1.
ν leads to decreasing of the switching time.When the value of ν is small, there is a significant difference between the dependencies j(t) at E = 0.95 • E c and E = 1.05 • E c .The variation of the temperature in the vicinity of T c does not change the behaviour of j(t) significantly.
When the temperature decreases and the value of δ becomes smaller than 1, then the minima corresponding to the two ferroelectric stable domains are placed below the minimum corresponding to the paraelectric phase at zero electric field, and the form of the free energy for the first-order phase transition is similar to the one for a second-order phase transition.It is interesting to compare the form of the dependencies j(t) for first-and secondorder phase transitions.Such a comparison is presented in Figure 4.One can see that the form of the time dependence of the current in the process of the polarization reversal for the first-order phase transition differs from the one for second-order phase transitions.According to the results of the calculations, the time dependence of the switching current has two peaks at second-order phase transitions.Unlike this, the time dependence of the switching current for first-order phase transitions has three maxima.Two of them are localized at the same points as in the case of the second-order phase transition.The third maximum lies between initial and final moments of the switching process and is connected with the existence of the paraelectric phase (stable or metastable).The calculations show that in the temperature range T 0 < T < T c , when T tends to T 0 , the form of the time dependence of the switching current for first-order phase transitions becomes similar to the one for second-order phase transitions.The position of the internal peak tends to the moment of completion of the reversal process, and the limiting form of the time dependence of the current has two maxima.The appearance of the third internal maximum should occur in the narrow temperature range.For this reason, this effect is difficult for observation.
The temperature dependence of the switching rate for the system undergoing first-order phase transitions in the presence of an applied electric field has some interesting peculiarities.To understand these peculiarities, let us consider the process of switching in more detail.In the presence of the applied field, the mechanism of the switching is determined by the value of the field.If the applied field E is larger than the coercive field E c , then the switching rate is determined by the polarization switching due to reorientation of all the dipoles via the interaction with the applied field which is independent of the interaction between the dipoles.On the other hand, for E < E c the switching is caused by the interaction between neighboring dipoles.In this range of field values, the polarization reversal occurs due to domain wall type motion caused by the interaction between the dipoles.The switching rate is significantly smaller in this case.Now we can calculate the dependence of E c on temperature for the systems undergoing first-order phase transitions according to (9).This gives This dependence is represented in Figure 5.The nonmonotonic behaviour of the E c (δ) is explained by the fact that the value of the coercive field at small δ is determined by the minima in the vicinity of P * = − 0.5 + √ 0.25 − a, while at rather large values of δ the main role plays the minimum in the vicinity of zero.Note that the value (3/80/0.09)• δ is rather small at δ ∈ [0, 4/3] and expression (10) can be well approximated by the first two terms of power series expansion, that is, by piece-linear function.
Taking into account this behaviour of E c (δ), one can find the regimes which provide the nonmonotonic change of switching rate.Really, if we choose the value of the external field E within the interval [E 1 , E 2 ], where E 1 = min(E c (δ)), E 2 = E c (4/3), then the whole range of δ variation will be subdivided into three intervals: [0, δ 1 ], [δ 1 , δ 2 ], and [δ 2 , 4/3].In the first and third intervals, the applied field is smaller than the coercive field, while in the second one, on the contrary, E > E c (δ).In this interval, the switching is caused by the interaction of dipoles with the applied field, and one can to expect that the switching rate is rather high.In the interval [0, δ 1 ], the value of the applied field is smaller than E c (δ).Note that in this interval E c (δ) coincides with E + (δ) and E > min(E c (δ)) > E − (δ).The inequalities E < E + (δ) and E > E − (δ) mean that the dimensionless free energy F = (a/2)P 2 − (P 4 /4) + (P 6 /6) − EP has a local minimum in the vicinity of P * = − 0.5 + √ 0.25 − a and has no a local minimum in the vicinity of zero.Therefore, the switching from the ferroelectric phase with the polarization close to P * = − 0.5 + √ 0.25 − a to the ferroelectric phase in the second domain with the polarization close to P * = 0.5 + √ 0.25 − a cannot be produced by the applied field only.The switching occurs due to domain wall type motion caused by the interaction between the dipoles.At small ν, the rate of the polarization reversal is not high in this case.Finally, in the interval [δ 2 , 4/3] the inequality E − (δ) > E > E + (δ) is true and the free energy has a local minimum in the vicinity P * = 0 while the local minimum in the vicinity of P * = − 0.5 + √ 0.25 − a does not exist.Therefore, the polarization of dipoles changes from the vicinity of P * = − 0.5 + √ 0.25 − a in the vicinity of P = 0 by the action of the applied field and the further switching occurs due to domain wall-type motion.The rate of this switching is not high if the value of the coefficient of interaction between neighbours is not large.If the value of ν is rather large, the influence of the applied field will not be strong and there will not be a difference in the transition rates for various intervals on δ.
We can see all these peculiarities in Figure 6 where the switching rate is presented as a function of temperature.
Here the applied field is equal to E c (δ = 1.17).The increase in the long-range dipole-dipole interaction does not change the behaviour of 1/τ(δ) dependence.At the same time the increase in the coefficient of the interaction between neighbours enhances the switching rate in intervals [0, δ 1 ] and [δ 2 , 4/3] so that the sharp maximum in 1/τ(δ) dependence is smoothed out. Figure 7 demonstrates that the 1/τ monotonically depends on δ when the parameters corresponding to BaTiO 3 are used, but at small value of ν the dependence 1/τ(δ) has a maximum.
Note that there are no such peculiarities in the behaviour of the 1/τ(θ) dependence for the systems undergoing secondorder phase transitions.This is connected with the fact that the coercive field in this case is determined by the formulae According to this formula, E c is (θ 3/2 /32) • (1 + o(θ)) at θ → +0, remains equal to zero at θ < 0, and depends on θ monotonically.One can see that the behaviour of E c for the first-and second-order phase transitions is quite different in the vicinity of T c .Such a difference can be easily explained if we consider the variation of the form of the free energy density in this range.For the first-order phase transition, the free energy density has three minima at T < T c and at T > T c .When T = T c , the depths of all the minima are the same.
It is necessary to apply some field to obtain the free energy density with only one minimum.This field determines the coercive field.On the contrary, for the second-order phase transition the free energy density has two minima at T < T c whose positions tend to each other.To obtain the free energy density with one minimum, it is necessary to apply a small field whose value tends to zero when T → T c − 0. At T > T c the free energy density has one minimum; it is not necessary to apply any external field and the coercive field is zero.The monotonic character of the dependence (11) allows to subdivide all the temperature range into two intervals.In the first one, the applied field is larger than E c and polarization reversal is determined by the interaction of the dipoles with the field.In the second interval conversely E < E c (θ) and the switching process is provided by the dipoledipole interaction.Calculations show that the switching rate decreases monotonically with the temperature growth at arbitrary fixed applied field.
The proposed model allows to investigate the pyroelectric effect in systems undergoing first-and second-order phase transitions.For the system undergoing a first-order phase transition in the homogeneous state in the absence of an applied field, the global polarization in the ferroelectric phase has the form (8). Taking into account that a = (3/16)δ, it is easy to find the pyroelectric coefficient Note that the ferroelectric state becomes metastable at δ = 1, and the stable state is paraelectric one where P * = 0.The pyroelectric coefficient in this case is evidently equal to zero.
For the system undergoing a second-order phase transition, the polarization in the ferroelectric state is where a = (3/16)θ.The calculation of the pyroelectric coefficient then gives For θ < 0, the free energy F = −(a/2)P 2 + P 4 /4 + P 6 /6 has a single minimum at P * = 0 (the ferroelectric state does not exist) and the pyroelectric coefficient is zero.The dependences of the pyroelectric coefficient on temperature for first-and second-order phase transitions are given in Figure 8.Note that the behaviour of the pyroelectric coefficient for the second-order phase transition in the vicinity of θ = 0 differs from the one for the first-order phase transition in the vicinity of δ = 1.In the first case, the pyroelectric coefficient has a singularity of the form dP * /dθ = √ 3/8θ 1/2 in the vicinity of θ = 0.For the systems with first-order phase transitions, the pyroelectric coefficient has a discontinuity in the vicinity of δ = 1.

Summary
We have investigated cases of polarization switching between the metastable and stable ferroelectric domains and the polarization reversal in the ferroelectric thin film simulated by using a 2D ferroelectric nanosystem.We have calculated the electric field and temperature dependences of the switching rate for the system undergoing a first-order phase transition in the region of the presence of metastable states.
We have found the time dependence of the switching current for the systems undergoing first-order phase transitions and have shown that it differs from the one for second-order phase transitions.We have predicted the appearance of the third internal maximum in the dependence of the switching current on time in the region of metastability.
We have obtained a maximum or nonmonotonic behaviour (depending on the value of the parameter of the short-range interaction between the dipoles ν) of the temperature dependence of the switching rate in the case of first-order phase transitions.In the case of secondorder phase transitions, the temperature dependence of the switching rate is described by a monotonic function.In contrast to second-order phase transitions, in which the pyroelectric coefficient exhibits singularity at the phase transition temperature, we have found a discontinuity in the pyroelectric coefficient at the first-order phase transition.

Figure 3 :
Figure3: Dependence of the switching current j on time t (in units of γη/β 2 ) at polarization reversal process (transition from the stable state to the second stable state or from the first metastable state to the second one) in the 2D system.Calculations are carried out for temperatures in the vicinity of the phase transition (δ = 1.017, δ = 0.9835) and for the values of the applied field in the vicinity of coercive fieldE c (E = 0.95 • E c , E = 1.05 • E c ): (1) δ = 1.017,E = 0.95•E c ; (2) δ = 0.9835, E = 0.95•E c ; (3) δ = 0.9835, E = 1.05•E c .At δ = 1.017, the two minima corresponding to the ferroelectric metastable domains are above the minima corresponding to the paraelectric stable phase at zero electric field, while at δ = 0.9835 the situation is reversed.The number of dipoles is 64 × 64.The initial state corresponds to the ferroelectric state with a nucleus of the second domain of polarization oriented in the opposite direction (the size of the nucleus is 4 × 4).The values of coefficients D and ν correspond to those for BaTiO 3 , D = 200, ν = 40.

Figure 5 :
Figure 5: Dependence of the coercive field E c on the dimensionless temperature δ = (T − T 0 )/(T c − T 0 ) for the system undergoing a first-order phase transition.

1 E 200 Figure 7 :
Figure7: Dependence of the switching, 1/τ (in units of β 2 /γη ), from the first stable ferroelectric state to the second stable ferroelectric state (or from the first metastable state to the second one) for the 2D system undergoing a first-order phase transition at E = E c (δ = 1.17) on dimensionless temperature δ.The number of dipoles N = 64 × 64.The initial conditions describe the nucleus of the ferroelectric phase inside the domain with polarization of the opposite direction.The results are given for the values of D and ν corresponding to those for BaTiO 3 (D = 200, ν = 40, curve 1) and for D = 200, ν = 0.4 (curve 2).The large value of the coefficient of the interaction between nearest neighbours results in the smoothing out of the maximum.At small value of ν the maximum appears.

Figure 8 :
Figure8: Dependence of the pyroelectric coefficient on the dimensionless temperature for systems undergoing first-and secondorder phase transitions.The curve 1 describes the behavior of the pyroelectric coefficient for the second-order phase transition.There is a singularity in the vicinity of θ = 0 of the order of 1/2.The pyroelectric coefficient for systems with the first-order phase transition is described by curve 2. There is a discontinuity of the dependence (dP * /dδ)(δ) in the vicinity of δ = 1.