Computational Techniques for Modeling Time-Fractional Dynamics of Polarization Switching in Ferroelectrics ? ??

Abstract

Ferroelectric polarization switching under external exposure is a complex phenomenon, the study of which can be enhanced by the application of mathematical modeling and computer simulation techniques. The paper extends the Landau{Ginzburg{Devonshire{Khalatnikov approach to simulate polarization switching in ferroelectrics as materials exhibiting time memory e ects. The proposed mathematical model is expressed by a time-fractional semilinear partial di erential equation. An implicit iterative nite di erence scheme based on an approximation of the Caputo fractional derivative was constructed and then implemented in Matlab. The computational performance of derived algorithm is provided by a numerical analysis of a test-problem solution and demonstrated by results of computer simulations of polarization hysteresis in ferroelectrics under applied periodical eld.

Keywords

• Landau{Ginzburg{Devonshire{Khalatnikov model Ferroelectric polarization switching Time-fractional partial di erential equation Caputo fractional derivative Implicit nite di erence scheme

-

In recent years, numerous studies have been extensively conducted to explore physical systems, indicating irregularity, spatial scale invariance, self-similar behavior, hereditary properties and time memory e ects. One of the theoretical frameworks describing complex physical phenomena is given by the fractional di erential theory. Fractional di erential equations can be used for the mathematical modeling of time-dependent responses in complex-structured media under non-equilibrium conditions. Such dynamics is also referred to as non-classical or anomalous processes [ 1, 3, 29 ]. In the present study, we consider ferroelectric materials as model objects examined with an application of the apparatus of fractional di erential equations.

Ferroelectrics are classi ed as a subclass of promising polar dielectrics, which have a wide range of applications in microelectronics, acoustics, radio engineering, optoelectronics, piezo- and pyrotechnics, etc. The most relevant applications of ferroelectrics in science and technology are associated with general mechanisms of polarization switching and domain structure dynamics stimulated by an external exposure [ 22, 31 ].

Signi cant attentions have been given over the years to theoretical investigations and experimental justi cations of the domain boundaries moving, restructuring of domain con gurations and kinetics of polarization switching in ferroelectrics. For instance, the Kolmogorov{Avrami model has been developed to model the polarization reversal processes in ferroelectrics based on statistical approach [ 20 ]. Also, the non-Kolmogorov{Avrami model has been proposed to simulate ferroelectric domain structure dynamics in polycrystalline disordered thin lms [ 28 ]. A modi cation of the Kolmogorov{Avrami model applied to polarization reversal process stimulated by electron irradiation has been designed in our previous works [ 13, 18 ]. A wide spectrum of models constructed with approximation techniques have been reported in di erent studies, in particular, the numerically stable Preisach model [ 4 ].

For this study, it is of interest to model of polarization switching process in view of the importance that ferroelectrics possess self-similar domain structures and exhibit time memory e ects in the process of registration of dynamic responses. Among the various models, we can emphasize models based on fractional di erential approach used for description of polarization switching processes. Speci cally, the study [ 17 ] has been reported a fractional di erential modi cation of the Kolmogorov{Avrami model. Also, the fractional di erential analogue of the model in the injection mode has been proposed in [ 18 ]. An approximation model of polarization-electric eld hysteresis dependence has been presented in studies [ 10, 17 ] using computation of time-fractional derivative of polarization. In all these studies, fractional-di erential models of complex domain dynamics and polarization switching have been considered, which are justi ed by indicating time memory e ects, fractal properties and self-similar behavior of ferroelectrics. Furthermore, the simulation results based on fractional di erential models have demonstrated a better agreement with the experimental data than the integer analogs.

In addition, a fundamental phenomenological theory of the Landau{Ginzburg{ Devonshire can be applied for mathematical simulation of polarization switching processes in ferroelectrics [ 22, 26, 30 ]. Moreover, the Landau{Khalatnikov equation allows one to examine ferroelectric state and polarization dynamic behavior under external eld [ 24 ]. The our previous study has been done to investigate the Landau{Ginzburg{Devonshire{Khalatnikov model from both theoretical and computational points of view [ 14 ]. Mathematically the generalized Landau{Khalatnikov model is governed by an initial-boundary value problem for a time-dependent nonlinear partial di erential equation of reaction-di usion type. To give some idea of the bene ts of this approach the current study assumes the fractional di erential modi cation of generalized Landau{Khalatnikov model to describe time memory e ects arising in ferroelectrics during polarization switching.

The mathematical formulation of this model represents a whole class of initial-boundary value problems for time-dependent nonlinear fractional partial di erential equations of the reaction-di usion type. As far as analytical methods can be applied for a su ciently limited range of problems, numerical methods play a crucial role in mathematical modeling and computer simulation of anomalous reaction-di usion systems. In particular, to solve numerically these problems, one can use computational schemes based on the nite di erence approach [ 2, 5, 7, 12, 15, 16, 19, 21, 23, 25, 27, 32 ]. The application of nite di erence schemes for solving fractional di erential equations have some special features, which are attributed to providing a satisfactory approximation order and reducing rather intensive computations costs. This primarily depends on the chosen de nition of the fractional derivative and on the method of approximation. In practice, the de nitions of Grunwald{Letnikov [ 15, 16, 23 ] and Caputo [ 12, 23 ] are widely used. Therefore the construction and implementation of e ective numerical schemes is particularly important and decisive.

Hence, the present study was undertaken to develop a model of polarization switching in ferroelectrics by means of constructing a time-fractional modi cation of the generalized Landau{Khalatnikov equation and a nite di erence scheme for further computer simulations. 2

Mathematical Model

The phenomenological Landau{Ginzburg{Devonshire theory of ferroelectricity enables one to examine the relation between the polarization P and the external electric eld E [ 22, 30 ]. In addition, the Landau{Khalatnikov approach can be used to formulate a time-dependent model of polarization dynamics in ferroelectrics. We will therefore apply the generalized Landau{Khalatnikov equation based on the Landau{Ginzburg{Devonshire{Khalatnikov theory to modeling polarization switching processes in ferroelectrics taking into account time memory e ects. The fundamental thermodynamic model is described by a semilinear parabolic partial di erential equation. By inductive assumption, we can introduce a time-fractional modi cation of the Landau{Khalatnikov equation to model complex polarization dynamics and polarization switching process.

To be de nite, we also assume that an uniaxial ferroelectric crystal is under consideration and the polarization P is the order parameter describing its state. This implies that polarization P depends on one space variable x and sample surfaces are associated with the coordinates x = 0 and x = L as illustrated in Fig. 1.

Also, suppose that the polarization reversal process in a ferroelectric crystal is realized due to the application of the sinusoidal electric eld. Notice that the polarization and the intensity of applied electric eld are vector quantities. In these terms, only two feasible states for polarization " P and # P are simulated (1) (2) (3) after the whole 180 switching event related to the orientation of applied eld " E and # E respectively.

Whence, we can arrive at the mathematical problem statement described by an initial-boundary value problem for a time-fractional semilinear parabolic partial di erential equation: P jt=0 = 0; 0 x

L; =

P ; =

P ; 0 t

T T ; where is the order of the time-fractional Caputo derivative (de ned below), 2 (0; 1); t is the dimensionless time; T is the characteristic time in s; T is the observation time in s; P = P (x; t) is the polarization distribution in C/m2; D = kT ; k is the thermodynamic restoring force in F=(m s); is the gradient coe cient in m3/F; E = EekT ; Ee=E0 sin(!t) is the electric eld applied along the polar axis in V/m; a = eakT , b = ebkT and c = eckT ; ea in m/F, eb in m5/(C2 F), ec in m9/(C4 F) are the thermodynamic constants; is the extrapolation length in m.

By construction, note that the characteristic time T is introduced to the model to adjust the dimensions of time and distance (see, e.g., [ 19 ]). In addition, we consider the model of polarization reversal process in ferroelectrics with rst order phase transitions. In the case of ferroelectrics with the second order phase transition, equation (1) can be reduced to a more simple form speci ed as a cubic time-fractional partial di erential equation with parameter b < 0.

Since the proposed model is expressed by a nonlinear fractional-di erential equation, the design of numerical schemes for such problems is of particular importance.

Outline of Algorithm

The accuracy of numerical solutions of fractional di erential equations is directly related to the de nition of the fractional derivative and the method of its approximation. A fractional derivative in contrast to a classical one, does not have an unambiguous de nition. In our case, we consider the Caputo time-fractional derivative [ 23 ], which is de ned as follows:

DC f (t) = (n 1

Z t where is the Gamma function.

In recent decades, the theory of fractional di erential problems has stimulated many investigations. Among the applications a special place is occupied by the numerical approaches for solving di usion-type equations [ 2, 5, 12, 15, 16, 21, 25, 27, 32 ]. Numerous studies have been devoted to numerical solutions of the anomalous di usion equation with a time-fractional derivative [ 2, 5, 12, 21, 25, 32 ] or with space-fractional derivatives [ 15, 16, 21, 27 ] or both of them [ 21 ]. In these studies numerical solutions of di usion problems have been obtained using the de nitions of Riemann{Liouville, Caputo, Grunwald{Letnikov, Ritz, etc.

The order of the approximation of fractional derivatives in most of the reported studies are less than the second. More accurate approximations have been derived for advection-di usion equations (see [ 2, 12 ]) and time-fractional subdi usion equation [ 5 ], as well as by constructing an algorithm for the numerical solution of the di usion-wave equation [ 7 ].

In order to obtain a numerical solution of the problem (1){(3) we derive a computational scheme based on the Caputo de nition of fractional derivative (4), an implicit nite di erence method and an iterative procedure. Here we apply the approach to the nite di erence approximation of fractional derivative proposed in study [ 12 ], in which numerical solution of the advection-di usion equation with a time-fractional derivative has been found, and also stability and convergence of the obtained di erence scheme have been investigated.

Hence, for the Caputo fractional derivative of order 0 < < 1 we can introduce an approximation on a nite-di erence grid = tk = k ; k = 0; N : + (5) d f tk dt = (3 ) k 1 X j=0 wk j f j+1 f j 1 + sk j f j+1

2f j + f j 1 +O 3 ; k = 1; N ; where is the time step; f (t) is a di erentiable function; the weight functions w and s are given as follows: wk j = 2 2 (k j)1 (k j 1)1 ; sk j = (k j)2 (k j

The fractional derivative (5) is approximated with the order of O( 3 ). However, if j = 0 in (5), we will have the function value f 1 at the dummy node, which is de ned outside of the computational interval. We can approximate this value by following f 1 = f (0) + O( ). As a result, the general scheme (called also L1-2 formula [ 8, 12 ]) is much more e ective and more accurate than the L1 formula, which is derived using direct nite-di erence approximation of the Caputo derivative [ 6 ].

Using the idea of constructing an implicit nite-di erence scheme in conjunction with an approximation of the Caputo time-fractional derivative [ 12 ], and a nite di erence approximation of the second-order space derivative, we can formulate a nite-di erence analogue for the di erential equation (1) on the space-time grid h = xi = ih; i = 0; M ; tk = k ; k = 0; N :

Pij 1 + sk j Pij+1 2Pij + Pij 1 i = (3 k 1 X hwk j P j+1

i j=0 ) D h2 =

Pik+1 2Pik + P k i 1 + aPik + b(Pik)3 c(Pik)5 + Ek; (6) where k = 1; N , i = 1; M 1.

Since we solve the cubic-quintic partial di erential equation, we arrive at a system of nonlinear di erence equations on each time layer. In this way we can use of an iterative procedure, which will allow us to solve a system of linecaornvaelrggeibnrgaitco ePquikatfioornse.acSho ttihmaet wmeo mfoernmt atk,sekqu=enc1e; Nof, aupspinrogxitmheatiroenlastiPoni(qs): (Pi(s))3 (Pi(q 1))2Pi(q); (Pi(q))5 (Pi(q 1))4Pi(q), where q = 1; 2; ::: is the number of iteration.

The iterative algorithm at the k time step starts by estimating an initial vPaiklue1.oTfhpeolcaorimzabtiinoantiuosninogf thneitveadluie efrreonmcethscehpermeveisowusitthimane sitteepra, ttihveatpirsoPcei(d0u)r=e enables one to solve applied problems without losing the accuracy of the general computational scheme (see, e.g., [ 9, 11, 14 ]).

Also, we use the initial condition Pi0 = 0, i = 0; M and the asymmetric nite di erence approximation for Robin boundary conditions (3) for k = 0; N : 3P0k + 4P1k 2h

P k 2 =

P k 0 ; 3P Mk 4P Mk 1 + P Mk 2 = 2h

P Mk : (7)

The general system of algebraic equations is solved by the Gauss method, which guarantees solving the problem with the error determined by the accuracy of computations. The design of the computational algorithm emphasis the relation between the applied framework of fractional derivatives and the formalized long time memory physical process. The polarization value at a given time moment depends not only on the behavior of the polarization in the vicinity of this point, but also on the values from the entire range of the temporal variable. Hence we proposed the implicit iterative nite di erence scheme based

This problem has an analytical solution de ned as u(x; t) = x3t3. In order to verify this solution, one can use a simple relation for calculation of fractionalorder derivative of a power function v = tp [ 15 ]: d v dt = (p (p + 1) + 1) tp t : (8) (9) (10) on an approximation for the Caputo fractional derivative to solve semilinear time-fractional partial di erential equation. This algorithm was implemented in Matlab software which allows high processing rates working with matrices. 4 4.1

Computer Simulation Results

Test Problem and Numerical Analysis Let us rst consider a simple test problem speci ed as a time-fractional di usion equation with imposed initial and boundary conditions:

(4) (3:15) x3t2:15 3xt3; 0 < x < 1; 0 < t 1;

We conducted simulations based on the proposed numerical scheme (6), (7) and using an implicit nite di erence scheme based on the Grunwald{Letnikov approximation of the fractional derivative [ 15, 16 ]. The latter is one of the most popular computational schemes used in di erent applications. The order of approximation of this scheme is estimated to be O( +h2). The graphs are visualized as coordinate pro les of function u(x; t) calculated at the last time moment of observation of the process.

In addition, we performed a numerical analysis of solutions of the problem (8){(10). The accuracy of the results were estimated using the maximum norm = ku uek1 = kuk1, where ue is the numerical solution, u is the exact solution calculated for the last time moment t = 1. Figure 3 demonstrates the values of the relative errors in the double log scale with the variation of the number of nodes M , N along the x and t axes, respectively.

A comparative analysis of the errors suggests that the described scheme provides a converged solution and an acceptable accuracy. These observations indicate that we can apply the derived computational scheme to get accurate numerical simulation for this class of problems. It should be pointed out that the algorithm can be characterized as resource-intensive and time-consuming. As an example, for number of nodes equal to N = M = 160 the accuracy corresponds to 10 6 and the computation time approximately equals 2000 s. 4.2

Simulation of Polarization Hysteresis Loop In this section, we perform numerical simulation of the hysteresis dependence between polarization and applied electric eld arising in ferroelectrics. For instance, let us consider the following numerical example as a mathematical model describing a time-space distribution of polarization, which changes with the external sinusoidal electric eld: x = 100P; = 100P; 0 t 1:5: (12) (13)

Here, a set of model parameters are normalized arbitrarily to promote comparability of simulated quantities (in particular, in view of the existence of the P (E) hysteresis dependence) approximately corresponding to ferroelectric crystals with rst-order phase transitions.

The result of computer implementation of the model (11){(13) is presented in Fig. 4. The hysteresis dependencies P (E) are plotted with varying the order of the time-fractional derivative for values of space step h = 0:0143 and time step = 0:0214. These ndings suggest that a decrease in the order of the timefractional derivative leads to a narrowing of the ferroelectric hysteresis loop while maintaining its shape. In other words, polarization hysteresis loop has a more narrow shape for crystals with a signi cant time-memory e ect.

Thus, the use of a time-fractional derivative in modeling polarization switching process potentially allows one to control the results of computations due to a variation of the order of fractional derivative or adjust this parameter as a numerical characteristic of the time memory e ect in ferroelectrics to provide the better agreement of simulations results with experimental data. 5

Conclusion

In this article we proposed a time-fractional modi cation of the thermodynamic model of polarization hysteresis in ferroelectrics and numerical scheme for its computer implementation.

The modi cation of the generalized Landau{Khalatnikov model was described by an initial boundary value problem for a time-fractional semilinear partial differential equation. We derived an implicit iterative nite di erence scheme based on an approximation of the Caputo fractional derivative. The combination of nite di erence schemes with an iterative procedure allowed applied problems to be solved without losing the accuracy of the general computational scheme. The formalized computational algorithm reveals the key peculiarities of simulation of long time-memory dynamic processes. The polarization value at a given time moment depends not only on the behavior of the polarization in the vicinity of this point, but also on the values from the entire range of the temporal variable. This algorithm was implemented in Matlab software and validated using test-problem.

We performed numerical simulations of the hysteresis dependence between polarization and applied electric eld arising in ferroelectrics. Our ndings indicate that the use of the fractional di erential apparatus provides a more " exible tool" due to a variation of the order of fractional derivative. This can be used to provide the better agreement of simulations results with experimental data.

References

1. Bogolyubov , A.N. , Potapov , A.A. , Rehviashvili , S.Sh. : An approach to introducing fractional integro-di erention in classical electrodynamics . Moscow Univ. Phys . 64 , 365 { 368 ( 2009 )

2. Cao , J. , Li , C. , Chen , Y. : High-order approximation to Caputo derivatives and Caputo-type advection-di usion equations (II) . Fract. Calc. Appl. Anal . 18 ( 3 ), 735 { 761 ( 2015 )

3. Chen , W. , Sun , H. , Zhang , X. , Korosak D. : Anomalous di usion modeling by fractal and fractional derivatives . Comput. Math. Appl . 59 ( 5 ), 1754 { 1758 ( 2010 )

4. Dragosits , K. , Hagenbeck , R. , Selberherr , S. : Transient simulation of ferroelectric hysteresis . Proc. 3rd Int. Conf. on Modeling and Simulation of Microsystems , pp. 433 { 436 ( 2000 )

5. Dimitrov , Y. : Three-point compact approximation for the Caputo fractional derivative . Comm. App. Math. Comp. Sci . 31 ( 4 ), 413 { 442 ( 2017 )

6. Dimitrov , Y. : Approximations for the Caputo derivative (I) . Mathematics, 1 { 39 ( 2016 )

7. Ding , H.: A high-order numerical algorithm for two-dimensional time-space tempered fractional di usion-wave equation . Appl . Numer. Math. 135 , 30 { 46 ( 2019 )

8. Gao , G.H. , Sun , Z.Z. , Zhang, H.W.: A new fractional numerical di erentiation formula to approximate the Caputo fractional derivative and its applications . J. Comput. Phys . 259 , 33 { 50 ( 2014 )

9. Grenkin , G.V. , Chebotarev , A.Yu. , Kovtanyuk , A.E. , Botkin , N.D. , Ho mann, K.-H.: Boundary optimal control problem of complex heat transfer model . J. Math. Anal. Appl . 433 ( 2 ), 1243 { 1260 ( 2016 )

10. Guyomar , D. , Ducharne , B. , Sebald , G. : High frequency bandwidth polarization and strain control using a fractional derivative inverse model . IOP Publ. Smart Mater. Struct . 19 ( 4 ), 045010 ( 2010 )

11. Kovtanyuk , A.E. , Chebotarev , A.Y. , Botkin , N.D. , Ho mann, K.-H.: Optimal boundary control of a steady-state heat transfer model accounting for radiative e ects . J. Math. Anal. Appl . 439 , 678 { 689 ( 2016 )

12. Li , C. , Wu , R. , Ding , H.: High-order approximation to Caputo derivatives and Caputo-type advection-di usion equations . Commun. Appl. Ind. Math. 6(2) , 1 { 33 ( 2015 )

13. Maslovskaya , A.G. , Kopylova , I.B. : Analysis of polarization switching in ferroelectric crystals in the injection mode . J. Exp. Theor. Phys . 109 , 90 { 94 ( 2009 )

14. Maslovskaya , A.G. , Moroz , L.I. , Chebotarev , A.Yu. , Kovtanyuk , A.E. : Theoretical and numerical analysis of the Landau{Khalatnikov model of ferroelectric hysteresis . Commun. Nonlinear SCI . 93 , 105524 ( 2021 )

15. Meerschaert , M.M. , Tadjer , C. : Finite di erence approximations for fractional advection-dispersion ow equations . J. Comput. Appl . Math. 172 ( 1 ), 65 { 77 ( 2004 )

16. Meerschaert , M.M. , Tadjer , C. : Finite di erence approximations for two-sided space-fractional partial di erential equations . Appl . Numer. Math. 56 ( 1 ), 80 { 90 ( 2006 )

17. Meilanov , R.P. , Sadykov , S.A. : Fractal model for polarization switching kinetics in ferroelectric crystals . Tech. Phys. 44 , 595 { 596 ( 1999 )

18. Moroz , L.I. , Maslovskaya , A.G. : Hybrid stochastic fractal-based approach to modeling the switching kinetics of ferroelectrics in the injection mode . Mathematical Models and Computer Simulations 12 , 348 { 356 ( 2020 )

19. Moroz , L.I. , Maslovskaya , A.G. : Simulation of nonlinear pyroelectric response of ferroelectrics near phase transition: fractional di erential approach . Mater. Sci. Forum 992 , 843 { 848 ( 2020 )

20. Orihara , H. , Hashimoto , S. , Ishibashi , Y. : A theory of DE hysteresis loop based on the Avrami model . J. Phys. Soc. Jpn . 63 ( 3 ), 1031 { 1035 ( 1994 )

21. Podlubny , I. , Chechkin , A. , Skovranek , T. , Chen , Y. , Jara B . Matrix approach to discrete fractional calculus II: Partial fractional di erential equations . J. Comput. Phys . 228 ( 8 ), 3137 { 3153 ( 2009 )

22. Rabe , K.M. , Ahn , C. , Triscone , J.-M.: Physics of ferroelectrics: a modern perspective . Springer, Berlin ( 2007 )

23. Samko , S.G. , Kilbas , A.A. , Marichev , O.I. : Fractional Integrals and Derivatives . Theory and Application . Gordon and Breach Science Publishers, USA ( 1993 )

24. Song , T.K.: Landau{Khalatnikov simulations for ferroelectric switching in ferroelectric random access memory application . J. Korean Phys. Soc . 46 ( 1 ), 5{ 9 ( 2005 )

25. Sontakke , B.R. , Shelke., A.S.: Approximate scheme for time fractional di usion equation and its applications . Glob. J. Pure Appl . Math. 13 ( 8 ), 4333 { 4345 ( 2017 )

26. Starkov , A.S. , Pakhomov , O.V. , Starkov , I.A .: E ect of thermal phenomena on a second-order phase transition in the Landau{Ginzburg model . J. Exper. Theor. Phys. Lett . 91 ( 10 ), 507 { 511 ( 2010 )

27. Szekeres , B.J. , Izsak , F. : A nite di erence method for fractional di usion equations with Neumann boundary conditions . Open Math. 13 ( 1 ), 553 { 561 ( 2015 )

28. Tagantsev , A.K. , Stolichnov , I. , Setter , N. , Cross , J.S. , Tsukada M.: NonKolmogorov{Avrami switching kinetics in ferroelectric thin lms . Phys. Rev. B . 66 ( 21 ), 214109 ( 2002 )

29. Uchaikin , V.V. : Anomalous di usion and fractional stable distributions . J. Exp. Theor. Phys . 97 , 810 { 825 ( 2003 )

30. Wang , C.L. , Zhang , L. , Zhong , W.L. , Zhang , P.L.: Switching characters of asymmetric ferroelectric lms . Phys. Lett. A . 254 , 297 { 300 ( 1999 )

31. Whatmore R. : Ferroelectric materials. Springer, Cham ( 2017 )

32. Xiaoting , L. , Sun , H. , Zhang, Y. , Fu. , Z-J. : A scale-dependent nite di erence approximation for time fractional di erential equation . Comput. Mech . 63 ( 3 ), 429 { 442 ( 2019 )

33. Zhang , B. , Ducharne , B. , Guyomar , D. : Inverse model of the piezoelectric ceramic polarization under wide bandwidth mechanical excitations with fractional derivative consideration . Opt. Quant. Electron . 46 ( 1 ), 103 { 110 ( 2014 )