Time-fractional Landau–Khalatnikov model applied to numerical simulation of polarization switching in ferroelectrics

Ferroelectrics are complex structured media that exhibit fractal properties and memory effects in terms of a number of dynamic characteristics. The most relevant applications of ferroelectrics in science and technology are associated with the primary mechanisms of polarization switching and domain structure dynamics induced by external exposure. In this study, we propose a time-fractional modification of the Landau–Ginzburg–Devonshire–Khalatnikov model to describe the dynamics of ferroelectric polarization switching. To solve the time-fractional cubic-quintic partial differential equation numerically, a computational scheme is derived. The technique combines an iterative procedure and an implicit finite-difference scheme based on an approximation of the Caputo derivative. A series of computational experiments are presented to visualize polarization switching characteristics on the example of thin films of barium titanate. A variation of the order of fractional derivative as a numerical characteristic of the memory effect allows one to “adjust” suitable regimes of simulated dynamical system. The obtained findings indicate that the generalized time-fractional model can provide better reproduction of experimental data in comparison with the classical analogue.


Introduction
Nowadays fractional differential equations are widely used for mathematical modeling of dynamic responses of complex structured physical media arising under non-equilibrium conditions [1][2][3]. Such processes can be referred to as non-classical or anomalous, and they are accompanied by significant gradient changes in the analyzed characteristics or a very long waiting time for aftereffects. Generally, memory effects in physical systems can be modeled using time-fractional partial differential equations (PDEs), whereas spacefractional PDEs are applied to describe dynamical processes in materials with nonhomogeneous structure and multiphase composition. Solid structures, disordered, random, and amorphous media, colloid materials are examples of such systems, which exhibit a complex geometry of bulk structure and surface topography, as well as a fractal behavior of characteristics of dynamical responses [4].
Specifically, ferroelectrics as a special class of promising polar dielectrics demonstrate complex scaling of domain configurations, self-similarity of nucleation processes and memory effects during the process of polarization switching, and common fractal properties of dynamical responses of crystals. Fractal features of domain structures, dielectric responses, and characteristics of polarization switching in ferroelectrics have been found by a number of independent researchers. Fractal analysis can be considered as a methodology for the study of the nature of ferroelectric phenomena and their most important properties, such as dielectric permittivity, light scattering, dielectric losses, hysteresis, acoustic emission, etc. Numerous works have established that domain configurations of ferroelectric systems possess both fractal structures and a self-similar character of nucleation and development of domains (see [5,6] and references therein). In particular, the application of methods of fractal and multifractal analysis of raster images allows one to explore the selfsimilarity of ferroelectric domain structure geometry, the complexity of domain configurations, and irregularity of domain boundaries, including those under nonequilibrium conditions [7][8][9][10][11]. Also, fractal and multifractal analysis of time series underlies the estimation of scaling characteristics of polarization switching dynamics and dielectric responses [12][13][14].
Despite the fact that in recent decades the properties of self-similarity of the geometry of domain configurations and the memory effects in ferroelectrics have been intensively studied using the methods of fractal analysis, the application of the apparatus of fractional differential calculus is not so developed. At the same time, it is obvious that the differential approach allows one to explore and predict the dynamical characteristics of the physical system. Certain fractional differential models have been proposed by independent research groups earlier to examine dynamic responses in ferroelectrics. For instance, the modification of the Kolmogorov-Avrami model for the estimation of polarization switching current in ferroelectrics has been developed on the basis of analytical calculation of fractional derivative in [15]. The fractional relaxation and oscillation equations have been introduced for ferroelectrics in [16]. Moreover, the fractional analogue of the Cole-Cole equation describing the dielectric relaxation has been considered in [17]. This approach of solving the problem has involved the use of finite-difference approximations of fractional derivatives. The similar numerical technique has been applied to implement the quasi-static model of ferroelectric hysteresis in [18,19]. The fractional-order derivatives have been introduced to describe the nonlinear mechanical and electric behavior of ferroelectric polymer composites in [20]. Moreover, a numerical solution of an fractional ordinary differential equation underlies the model of polarization switching current induced in ferroelectrics by electron beam irradiation [11,21].
The most significant applications of ferroelectrics in science and technology are associated with the general mechanisms of polarization switching and dynamics of domain structures. A challenging problem which arises in this field is predicting polarization responses induced by external exposure in ferroelectric materials. The theoretical description of characteristics of ferroelectric polarization switching is supported by the framework of the Landau-Ginzburg-Devonshire phenomenological theory. The thermodynamic formalism for ferroelectrics have been reported in many reviews [22]. The Landau theory of phase transitions has been previously used for modeling of polarization reversal processes and ferroelectric hysteresis loops in [23]. The corollaries of the Landau-Ginzburg-Devonshire theory can be employed to estimate the ferroelectric characteristics by means of the Ising and Heisenberg models implemented with the use of Monte Carlo simulation and Metropolis algorithm [24,25]. Furthermore, the time-dependent Landau-Khalatnikov equation has been proposed to describe the ferroelectric state and dynamics of polarization reversal under external field [22]. This model has attract more attention in the field of study of switching characteristics and hysteresis phenomena in ferroelectric thin films [24,26], the partial switching and external stress effect [27], the peculiarities of ferroelectric phase transitions in the view of electrocaloric and pyroelectric effects [28], etc.
Generally, the mentioned approach is based on the mathematical apparatus of differential equations. In the present study, we will focus on the Landau-Ginzburg-Devonshire-Khalatnikov model, described by the initial-boundary value problem for cubic-quintic partial differential equation (PDE). Our previous study has been devoted to the exploring this model from both theoretical and computational points of view [29]. As shown above, the use of the fractional differential con-cept provides a generalization of integer models and significantly expands the range of functional capabilities of the mathematical modeling of physical systems. Since ferroelectrics possess complex polarization structures and demonstrate memory effects, we propose a time-fractional modification of the Landau-Ginzburg-Devonshire-Khalatnikov model to describe the polarization switching processes in ferroelectrics. This approach gives a significant advantage because it allows us to use an additional "flexible tool" for modeling of the space-time distribution of polarization in ferroelectrics and predicting dynamic responses. Obviously, analytical solutions of semilinear fractional PDEs meet serious problems.
In the context of specific applied problems, numerical methods are usually used, in particular, computational schemes based on the finite-difference approximations of fractional derivatives [30,31]. The preliminary result has been reported in [32], where we applied the simplified computational scheme and solved numerically the test-problem related to the generalized Landau-Khalatnikov model (without the using units of physical quantities).
Since exact methods can be used for a sufficiently limited range of fractional problems, numerical methods have received substantial interest for their potential use in different applications. Presently, various numerical methods are being developed. Specifically, one can apply computational schemes based on the finite-difference approach [33][34][35][36][37], random walk models [38], the finite element method [39], Monte Carlo simulation [40], the matrix transform method [41], the method of Adomian decomposition [42], the timespace spectral method [43], etc. In the present study, we will focus on the application of a finite-difference approach to derive a computational algorithm and implement the mathematical model.
The current study was undertaken to develop the Landau-Ginzburg-Devonshire-Khalatnikov model and the computational algorithm for its implementation with application to the study of polarization switching characteristics and hysteresis in ferroelectric materials. The theoretical basis and governing equations of the generalized Landau-Khalatnikov model describing the polarization states of ferroelectrics will be presented in Sect. 2. Section 3 will be devoted to the construction of computational scheme for solving cubic-quintic fractional PDE. The Sect. 4 will focus on the computer simulation of polarization switching characteristics for typical ferroelectrics with the first-order phase transition. Here the performing full cycle of mathematical modeling and computer simulation to solve multidisciplinary problem also underlies the main significance of the present study. Finally, in Sect. 5, we will summarize main findings and present concluding remarks arising from this study.

General description of the mathematical model
In this section, the mathematical problem statement designed for modeling of polarization switching in ferroelectrics is presented. The emphasis is placed on case of time-fractional modification of thermodynamic model describing dynamical responses in ferroelectrics.
A theoretical analysis of the properties of polarization switching and phase transitions in ferroelectrics can be performed on the basis of the Landau-Ginzburg-Devonshire thermodynamic theory [22]. The role of the thermodynamic approach in the study of ferroelectrics can hardly be overestimated. The Landau-Ginzburg-Devonshire theory enables describing the dielectric, pyroelectric, piezoelectric, electrocaloric and other properties of ferroelectric materials, and in many cases it gives a satisfactory agreement with the experimental results.

Thermodynamic framework for modeling of polarization switching in ferroelectrics
According to the concepts of thermodynamic theory, the equilibrium state of ferroelectrics can be completely specified by a fixed set of quantities: polarization, electric field, temperature, mechanical deformation and stress. The state of the ferroelectric material is determined by the free energy function and the values of the dependent variables correspond to the minimum free energy for the equilibrium state of the system. In a particular case, when there are no applied mechanical stresses, the polarization P is directed along one of the crystallographic axes only. For uniaxial ferroelectrics, the free energy density F can be expanded in powers of a single polarization component given by the following ansatz [22]: where P is the spontaneous polarization in C/m 2 ; E is the field intensity in V/m; A = A 0 (T − T CW ), A 0 , B, C are the thermodynamic constants ( A in m/F; B in m 5 /(C 2 · F); C in m 9 /(C 4 · F); T is the current temperature and T CW is the Curie-Weiss temperature in K). The total free energy is calculated as G = V FdV in J. The change in the free energy between the paraelectric phase at T > T CW and the ferroelectric phase at T < T CW depends on the coefficient B. If B > 0, the second-order phase transition occurs, whereas B < 0 is related to the first-order phase transition. For the second-order phase transition, the term containing the sixth power in the expression for the free energy density (1) can be neglected. The values of the thermodynamic parameters of A 0 and C are positive for all known ferroelectrics [22]. Hereinafter for definiteness, we will consider ferroelectrics with the first-order phase transitions.
Further, we assume that polarization switching is induced by applying a periodical electric field where ω = 2π f is the radial frequency of the applied field in rad/s and f is the frequency of field oscillations in Hz. Isaak Khalatnikov has developed the thermodynamic approach to describe the dynamics of ferroelectric polarization under an external field [22]: where δ is the kinetic coefficient in m·s/F. The proposed equation is known as the nonstationary Landau-Khalatnikov equation. The general Landau approach is based on the self-consistent field theory, and as a consequence, the relation (1) does not consider the fluctuations of polarization. As proposed by Ginzburg, including the gradient term or correlation energy ψ P into the expression for free energy leads to the generalized analogue of the Landau-Khalatnikov equation given as follows: where D = ψ/δ is the thermodynamic parameter as combination of the gradient coefficient ψ and the kinetic coefficient δ, with the sense and the dimension of the diffusion coefficient in m 2 /s; σ = ν/δ; ν is the scaling factor; A = −A/δ, B = −B/δ, C = C/δ are the positive constants. The ferroelectric hysteresis loop can be estimated for certain values of the temperature T as a polarizationelectric field dependence calculated at successive times. In this way, the generalized Landau-Khalatnikov equation (4) and the periodical field intensity defined by (2) underlie the mathematical model of ferroelectric polarization switching in uniaxial ferroelectrics. A suitable set of boundary conditions and initial condition should be posed to close the mathematical problem. As will be shown, the specific form of the boundary conditions is due to the physical meaning of corresponding particular problems.

Mathematical problem statement
As mentioned above, ferroelectrics exhibit the complex nature of the domain structure rearrangement and memory effects observed during polarization switching. The mathematical description of dynamical responses in ferroelectrics due to time memory effects can be given with the use of the fractional differential approach as a way to generalize the possible variations of dynamical modes. As far as changes in polarization as a key dynamical characteristic of polarization reversal process in ferroelectrics can be analyzed by means of differential models, including the Landau-Ginzburg-Devonshire framework, here we develop a time-fractional modification of the Landau-Ginzburg-Devonshire-Khalatnikov equation. Such generalization of the thermodynamic approach is bearing in mind the variation of dynamical regimes as well as their model representations that are associated with time dispersion processes arising in ferroelectric materials.
Let us consider the electric field applied parallel to the general direction of polarization P. In this case, the polarization P is an order parameter, which depends on one space variable x. The surfaces of ferroelectric sample are specified by x = 0 and x = L. Note that the polarization and the applied electric field are vector quantities. This means that only two possible states of polarization ↑ P and ↓ P are modeled during a complete switching of polarization by 180 • , associated with the orientation of the applied field ↑ E and ↓ E, respectively.
Although classical differential equations allow one to describe the dynamic processes in different systems, they do not formalize for such properties as a memory effect. Fractional integro-differential calculus gives rise to describing processes in systems for which nonlocal properties in time and space should be taken into account. Interpreting fractional derivatives as a way of considering the memory effect (the proper time nonlocality) and spatial correlations (the coordinate nonlocality) has led to their widespread use in natural sciences. Nonlocality in time is contained in the kernel of the integral operator. In addition, many authors have noted that the dynamical contribution integer time derivative results in an overestimation of calculated values compared to experimental data attributed to ferroelectric materials [15,18,19,44]. This problem is corrected by introducing in model derivatives of noninteger order. Therefore, we propose to proceed to the generalization of the Landau-Ginzburg-Devonshire-Khalatnikov model by introducing a time-fractional derivative that includes the case of integer differentiation as a particular.
The time-fractional analogue of the Landau-Ginzburg-Devonshire-Khalatnikov model is represented as an initial-boundary value problem for the following semilinear time-fractional partial differential equation: where are the combinations of model parameters; τ = t \ t * is the dimensionless time; t * is the characteristic time of the process in s; q and g are model parameters; θ is the time of process observation in s; L is the thickness of the film in m; ∂ α P ∂τ α is the left-hand Caputo derivative; α is the order of time-fractional derivative, 0 < α < 1.
We stress that, initially, dimension units in fractional models are not agreed due to using fractional operators of different orders. For example, the dimension of the diffusion coefficient in the time-fractional equation has to be is the dimension of distance and [t] is the dimension of time. In this way, solving problems from specific subject areas requires additional unit coordination. Here we used the normalization of the equation with respect to time and the value of the diffusion coefficient D with the standard The characteristic time t * is introduced to provide an agreement of the physical dimensions and units of measure of variables and parameters of the model. In order to analyze the results of computer simulations, the reverse transformation is required (from the dimensionless time τ to the dimensional analogue t).
Within the framework of the current study, the mathematical problem statement (5)-(7) allows us to consider two principle problems. First, the model of polarization switching in ferroelectric thin films can be specified at zero initial condition P 0 (x) = 0 and the boundary conditions (7) should be transformed to the following view [22]: where λ is the extrapolation length in m.
On the other hand, we can also formulate the model of motion of 180 • domain wall separating two domains with opposite polarization orientation. Here we can assume that polarization switching is induced by a linearly increasing or trapezoidal electric field. The model needs the initial value of polarization distribution to be specified, for instance, with the use of approximation P 0 (x) = P 0 tanh(x \ η), where η is the parameter related to the thickness of ferroelectric domain boundary.
Notice also that different variations of the fractional time-dependent Landau-Ginzburg equation (the socalled TDLG equation) have been in the focus of several research efforts. Such modification was first proposed by Weitzner [45] and then developed by Tarasov [46]. In a study by Li et al. [47], a theoretical analysis of the existence and uniqueness of a weak solution to the TDLG equation has been carried out, and the existence of global attractors has been shown. Basically, the detailed study of the Landau-Ginzburg-Devonshire-Khalatnikov equation presents a separate research subproblem and it is beyond the scope of present study.
In general, the fractional cubic-quintic equation (5) can be referred to as semilinear equations of anomalous diffusion, specifically, subdiffusion equations. The construction of computational algorithms for the implementation of the mathematical model, due to their specificity, requires the application of numerical methods.

Numerical techniques
Current section presents computational aspects of model implementation with a focus on construction of the numerical algorithm based on the joint application of a finite-difference method and an iterative procedure to solve a time-fractional semilinear PDE.
In order to derive a computational algorithm, let us reformulate the mathematical problem statement in the following generalized view: where u = u(x, t) is unknown function, ϕ(u) = au + bu 3 − cu 5 is semilinear reaction part of the time-fractional PDE, f (x, t) is the source function, a, b, c, κ, q, g are positive constants, and α is the order of the time-fractional derivative in the Caputo sense.
To be more specific here, we consider the left-hand side Caputo fractional derivative defined as is the integer part of α, 0 ≤ {α} < 1 is the fractional part of α, and (·) is the Gamma function.
Further, we define a finite-difference grid t = t j = j t, j = 0, M . In practice, the L1 approximation is most often used as an approximation of the Caputo derivative. The order of the approximation is estimated as O( t 2−α ). At the same time, schemes of higher order can be derived based on the L1 approximation. For example, a new fractional numerical differentiation formula (the L1-2 formula) has been obtained at [48] to approximate the Caputo fractional derivative of order 0 < α < 1. The computational efficiency and numerical accuracy of the new L1-2 approximation are better than for the L1 approximation. The study [49] has also suggested a modification of the L1 approximation, in which the first three weight coefficients have been modified by using the Riemann zeta function: where 0 < α < 1 and the weight functions μ l are given as: ζ is the Riemann zeta function.
Let us remark that the application of the Caputo fractional derivative (unlike the Riemann-Liouville fractional derivative as an example) in constructing finitedifference schemes leads to an additional requirement for smoothness of the function. We follow here the approach widely accepted elsewhere [35,48,49], in which the function is assumed to satisfy appropriate smoothness condition. According to [48,49], suppose f (t) ∈ C 2 [0, t M ] to provide an accurate approximation.
The considered approximation has been applied previously to solve the initial-boundary problem for linear time-fractional subdiffusion PDE in [49]. Here we expand this approach to solve the semilinear problem (9)- (11). Let the space-time grid t x = {x i = i x, i = 0, N , t j = j t, j = 0, M} cover the computational domain, where x is the coordinate step and t is the time step.
In contrast to the continuous problem with the solution defined on the whole domain, the finitediscretization results in a discrete solution specified only at a finite number of grid points. To deal properly with these functions, we introduce the discrete function space of grid functions U h , which is isomorphic to finite-dimensional Euclidean spaces. Further, to perform the numerical convergence analysis of a finite-difference method, the space of grid functions is considered to be equipped with an appropriate discrete norm (for instance, l 2 norm is further used).
Similar to (12), we approximate the time-fractional partial derivative of the function of several variables u(x, t) as below: Based on the introduced approximation of temporal derivative and taking into account the centered finitedifference approximation of the spatial derivative, we arrive at the general approximation for equation (9): where i = 1, N − 1, j = 1, M. Therefore the implicit scheme for the numerical solution of the problem (9)-(11) can be written as: − κ where i = 1, N − 1, j = 1, M.
As far as we solve the cubic-quintic PDE, we obtain a nonlinear system of equations on each time layer. In this case, we can combine the finite-difference scheme with the iterative procedure. This means that we form a sequence of approximations u .. is the number of iteration. The iterative algorithm at the j time step starts by estimating an initial value of function using the value from the previous time step, that is u . Note that the combination of finite-difference method with an iterative scheme allows one to keep the level of accuracy corresponding to the order of approximation of the applied numerical method [50,51].
To approximate the boundary conditions (11) with the second order of accuracy, the asymmetric finitedifference formulas can be used where j = 1, M. The initial condition (10) is also included into the computational scheme. The general system of algebraic equations is solved by the Gauss method for each time layer.
Using the approximation, we can write the local truncation error R j i of the difference scheme (15): Presently, researchers are paying more and more attention to finding a solution to the fractional Landau-Ginzburg equation using numerical methods, in particular, grid methods, weighted residuals, finite elements, etc. [52][53][54]. Analysis of a numerical solution based on the implicit scheme (12) for linear fractional subdiffusion equation has been shown in [49]. Also, the absolute stability of implicit finite-difference schemes based on the Caputo definition has been examined elsewhere [35,55,56]. In general, the theoretical analysis of the finite-difference schemes derived for nonlinear fractional differential problems meets essential difficulties. The theoretical analysis of finitedifference schemes for the initial-boundary values problems for cubic-quintic time-fractional Landau-Ginzburg-Devonshire-Khalatnikov equation requires further multifaceted consideration and serious independent study. In the framework of present study, we apply numerical approaches to analyze stability, the resulting errors and the rate of convergence.
In addition, we note that our approach involves modeling of a complex dynamics of physical system based on numerical solution of the time-fractional PDE representing a modification of the thermodynamic model, whereas the studies [18,19] have proposed the fractional model of ferroelectric hysteresis based on the numerical approximation of fractional derivatives.

Computer simulation results and discussion
In this section, we apply the constructed numerical scheme for the computer simulation of polarization switching in ferroelectric thin films with the first-order phase transition.

Computational experiments setup
The constructed computational technique was implemented using Matlab programming platform. The application program interface is intended for modeling of the space-time distribution of polarization and polarization-electric field dependence during the process of ferroelectric switching. To perform simulations, we need to initialize characteristics of the material (thermodynamic parameters, the thickness of the ferroelectric sample), the observation time, and the applied field. The number of nodes of the space-time grid covering computational domain is also need to be specified.
The most studied case of polarization switching is 180 • polarization inversion, in which an external field is applied along the spontaneous polarization axis and its orientation changes directly by 180 • . For example, 180 • polarization switching is observed for almost all uniaxial ferroelectric crystals, such as triglycine sulfate (NH 2 CH 2 COOH) 3 H 2 SO 4 , potassium dihydrogen phosphate KH 2 PO 4 , lithium niobate LiNbO 3 , etc. Polarization switching process can be realized by means of two relevant mechanisms, namely nucleation and growth of 180 • domains with spontaneous polarization parallel to applied field as well as lateral expansion of oppositely orientated domains.
Although the generalized Landau-Khalatnikov model (5)-(8) is most suitable for uniaxial crystals, it can also be used to describe polarization switching in multiaxial crystals. In particular, barium titanate (BaTiO 3 ) has both 90 • and 180 • oriented types of domain structures in the tetragonal phase at room temperature. Using an appropriate external field applied along the c-oriented domains of BaTiO 3 , 180 • polarization switching can be observed [22,57,58]. In the current numerical experiments, we consider switching of ferroelectrics under the periodic electric field (2). The simplified scheme of domain structure and 180 • polarization switching in c-domains of BaTiO 3 is illustrated in Fig. 1. The initial polarization is associated with the original equilibrium state of ferroelectric, which equals zero at the start time moment t = 0. Here we assume that the complete 180 • switching event is induced by an applied field. The intermediate polarization states during the switching process match to range of field variation (E 0 , E 0 ). Table 1 lists the values of thermodynamic parameters for barium titanate at temperature T (in K) [22,59] used for computer simulations. The parameters of the computational experiments correspond to the experimental data described in [60]. In the tests, the parameters are set to be: the amplitude of the field intensity is E 0 = 2.5 × 10 5 V/m, the field frequency is f = 50 Hz (ω = 2π · f ), the gradient coefficient is ψ = 5 × 10 −8 m 3 /F, the kinetic coefficient (which determines the velocity of movement of domain walls) is δ = 2 × 10 5 m·s/F [61]. The temperature needs to be specified: "The observation time is θ = 0.03 s, the temperature is T = 300 K, the characteristic time of the process is t * = 1 s, and the sample thickness is varied.
By the above, to visualize the polarization-electric field dependence for each time moment, we calculate the averaged value of the polarization over the crystal thickness:

Numerical study of convergence and stability
In order to verify the obtained simulation results, we applied several approaches. First we performed the comparison of numerical solutions with exact solutions for test linear fractional problems. Further, when solving nonlinear fractional problems, the verification of the simulation results was realized by means of comparison of the simulation result with the data of an integer analogue of the model in the limiting modes (e.g. at α = 0.999). Finally, as we will show further on, the adequacy of the results is established by comparing of simulation data with experimental observations of polarization hysteresis in ferroelectrics. In addition, we conducted a study of the practical convergence of the scheme for numerical solution of the fractional Landau-Ginzburg-Devonshire-Khalatnikov equation. For this purpose, we applied the testing scheme by the double-counting procedure in two versions. The numerical convergence study is based on refining the spatial and temporal steps of the grid covering solution domain separately. This approach is widely used for numerical testing schemes in terms of the study of their convergence rate [62,63].
In the first case, we estimate the spatial relative error using l 2 -norm as follows: whereP(x, t) ∈ U h is the polarization calculated at the spatial partitions number N and the polarizatioñ P(x, t) ∈ U h corresponds to the doubled number of spatial partitions 2N with the number of time partitions M = 500.
The test results are shown in Table 2. The computations were performed for the set of parameters mentioned above, the simple thickness is defined to be 1 μm, and the order of the fractional derivative equals 0.9. The observation period is set to be θ = 0.03 s.
Hence, we can conclude that the error decreases rapidly when the spatial grid is refined. A sufficiently acceptable accuracy (of the order of 0.01 %) can be obtained using a grid with the number of spatial partitions N = 100 − 200.
Further, the time relative error is calculated by the formula: whereP(x, t) ∈ U h is the solution obtained with the number of time partitions M andP(x, t) ∈ U h is associated with the doubled number of time partitions 2M. Computations were performed at N = 200, while the observation period is also specified as θ = 0.03 s. The dependence of relative error on the number of time partitions is illustrated in Table 3. We stress that a decrease in the time step also leads to a decreasing the error, but in terms of accuracy the numerical solution turns out to be more sensitive to a change in the time step compared to the spatial step.
As a result, although the use of a refined grid (for example, with the partitions number N = 1600) corresponds to an error of the order of 10 −6 , does not provide benefit in terms of accuracy, because the total error will be limited by the error caused by the time sampling. The mesh refinement is always associated with an increase in the computational time, for an already resourceintensive fractional differential problem. In computational experiments when controlling the prior accuracy (for example, 0.01 %), it makes sense to take the optimal ratio of steps in space and time as M = 2000, N = 200.
The iterative procedure included into general algorithm was also examined to verify the numerical convergence. In practice, three or four iterations on each time layer were enough to get accuracy sufficient compared to the total accuracy of the numerical solution (in this case ratio error is approximately equal to 0.001 %).  Moreover, let us demonstrate the result of a numerical study of the stability of a computational scheme. Since the initial polarization as the initial condition and the thermodynamic parameter λ in boundary condition are strictly defined, here we examine the stability of the computational scheme by means of variation of applied electric field. The key characteristic is presented by the amplitude of the electric field E 0 . Figure 2 shows the time dependencies of an electric field successfully increasing the field by 20%, 30%, 40%, 50% with respect to the initial value E 0 = 2.5 × 10 5 V/m. In order to estimate the perturbations of grid functions under the variation of the electric field amplitude, we use the following estimation: ξ m = P m − P m−1 2 , m = 2, 3, 4, 5, where m = 1 corresponds to the initial polarization. Figure 3 illustrates the estimation of the perturbation of grid functions. These observations suggest an appearance of slight perturbations of the polarization in the variation of the electric field. The linear growth of the estimation indicates the numerical stability of the algorithm.

Simulation of ferroelectric polarization switching in BaTiO 3 thin films
Let us consider the results of computational experiments on the example of modeling of the characteristics of the polarization switching in c-domains of barium titanate thin film, that is induced by sinusoidal electric field E(t). Figure 4 shows the simulation results of the dynamics of the averaged polarization over the crystal thickness P(t) during the polarization switching. The ferroelectric hysteresis loops as polarizationelectric field functional dependencies P(E) are demonstrated in Fig. 5. The simulation results are visualized for different values of order of fractional derivative α  [60] are also plotted in Fig. 6. In a series of computational experiments, the value of the fractional  [60] derivative order α was varied from 0.75 to 1 with a step of 0.05. From the results obtained, it can be concluded that a decrease in the order of the time-fractional derivative corresponds to a narrowing of the dielectric hysteresis loop while maintaining its shape.
In order to perform more detailed analysis, we estimate statistical characteristics of the data obtained by the fractal model to experimental data under variation of the fractional derivative order α. Concretely, estimations of standard deviation, the coefficient of determination, and average absolute percentage error are listed in Table 4. The fractional derivative order α is varied from 0.8 to 0.9 with a step of 0.025. Our findings suggest that the best agreement with the experimental data corresponds to the order of the fractional derivative α=0.85 with the value of the coefficient of determination equal to 0.985. The value of the average absolute percentage error of 14% gives an overall measure of the goodness of fit between simulated polarization and the corresponding experimental data and indicates here a good forecast accuracy.
In this instance, the observations indicate that the main characteristics of the hysteresis loop, such as coercive field E c = 0.95 × 10 5 V/m and residual polarization P r = 0.24 C/m 2 , coincide with the data described in the literature (e.g., [22,59,61,64,65]).
It should be pointed out that the physical meaning of the application of the fractional differential approach in our case is presented by the formalization of timedependent changes in polarization related to the ferroelectric domain structure restructuring mechanism.
Here the space-time distribution of polarization is specified by not a classical diffusion process, but in the language of the theory of anomalous diffusion, the corresponding regime can be associated with "slow" wandering. The latter can be interpreted as a deceleration of the speed scale (the function changes faster with a smaller time value than with a larger one). Also note that the primary designation of the Landau-Ginzburg-Devonshire-Khalatnikov model is a prediction of the state of a nonlinear system based on the underlying thermodynamic principles, but not an approximation of a polarization hysteresis loop directly as it takes place in application of the well-known Preisach model (as an example). In our case, the model parameters are the parameters of the physical experiment (the frequency, the amplitude of the electric field, the initial value of the polarization, etc.) and the thermodynamic characteristics of the crystal (specified for each ferroelectric material). The only strictly approximate parameter included in this Landau model is the scaling factor ν for the field, which has to be adapted for simulation of polarization hysteresis ferroelectrics in high frequency fields. An increase in this coefficient results in a narrowing of the hysteresis loop and its extension along the ordinate axis. Thus, the order of fractional differentiation α as a numerical characteristic of the memory effect is an additional control parameter of the fractional model of polarization switching, which can be tuned to provide better agreement with the experimental data.
In the aspect of studying the influence of the values of control parameters on the behavior of the simulated characteristics, further we consider the so-called size  effect, which is presented here by the dependence of the characteristics of ferroelectric hysteresis on the sample thickness. Figure 6 shows the hysteresis loops calculated by means of the generalized Landau-Khalatnikov model with varying the thickness of BaTiO 3 .
Numerical experiments indicate that the shape of the hysteresis loop does not depend on the crystal thickness L for bulk materials. As known, in this case the generalized Landau-Khalatnikov model can be reduced to the case of the classical model described by the Cauchy problem for an ordinary differential equation (the gradient term can be neglected for bulk objects). At the same time, for low-dimensional structures (nanowires, thin films, island coatings, etc.), taking into account the Ginzburg term and introducing Robin boundary conditions are relevant for simulations.
We can observe a significant distortion of the hysteresis loop with a decrease in the film thickness, such as deformation, narrowing of the loop, formation of the double loop. In our computational experiments, the film thickness of 3 nm is estimated to be critical, since the hysteresis loop disappears at this thickness. These results agree with the data obtained for nanowires and thin films of barium titanate by independent authors in [66,67].

Conclusion
In this paper, we performed the detailed numerical study of the fractional modification of Landau-Ginzburg-Devonshire-Khalatnikov model proposed to describe the polarization switching in ferroelectrics. To solve time-fractional cubic-quintic PDE, we constructed the computational scheme based on the combination of a finite-difference approximation of Caputo derivative and an iterative procedure. The computational algorithm was implemented in Matlab to perform numerical simulations of ferroelectric hysteresis dependencies of polarization on applied field. In order to verify simulation results, the numerical convergence study was conducted. As a result, this allowed us to estimate the optimal values of computational parameters for computational experiments.
We demonstrated simulation results on the example of barium titanate thin film. The polarizationelectric field dependencies for BaTiO 3 thin films were calculated by means of fractional Landau-Ginzburg-Devonshire-Khalatnikov model. The results of computational experiments were demonstrated at varying the order of fractional derivative. The observations suggest that the possibility of varying the order of fractional differentiation provides the results of the implementation of the differential model, which better reproduce the experimentally observed regularities in comparison with the integer analogues. As a numerical characteristic of the memory effect, the order of fractional derivative can be used as an additional control parameter of the dynamic model that can be tuned to ensure optimal agreement with the experimental data.
Thus, we can conclude that the fractional generalization of the Landau-Khalatnikov model formalized the dynamic responses of ferroelectrics under external exposure that significantly expands the range of functional capabilities of the numerical modeling methodology as applied to the study of this class of fractal physical systems with memory.