## Abstract

This chapter presents the prevalent mathematical models of irreversible processes in polycrystalline ferroelectric materials when they are subjected to intense electrical and mechanical influences. The main purpose of such models is to describe the dielectric hysteresis loops, with which the models of Rayleigh and Preisach coped well, though they were developed almost a 100 years ago. Nevertheless, in order to describe the whole gamut of material properties in irreversible polarization-depolarization processes, it was required in the last three decades to develop new approaches and methods that take into account the material structure and the physics of the process. In this chapter, we attempted to collect the most common one-dimensional models, with a view to give a brief description of the basics and approaches with the application of working formulas, algorithms and graphs of numerical calculations. On one-dimensional models, the basics of three-dimensional models are worked out, such as evolutionary laws, domains switching criteria, generalizations from “hysteron” to the polarization surface, and so on, so they are a necessary step in modeling. However, some of them proved to be so effective that they obtained the right to independent existence, as happened with the Preisach model, which found application in dynamic systems. This research is based on published articles, monographs, proceedings of conferences, and scientific reports of individual collectives published over the past 20–25 years.

### Keywords

- mathematical model
- hysteresis
- domains
- polarization
- strain
- ferroelectric materials
- ceramics
- constitutive relations

## 1. Introduction

Polycrystalline ferroelectric materials or ceramics are active materials that, by virtue of their internal structure, have the ability to convert mechanical energy into electrical energy, and vice versa. This means that by acting on the sample an electric field or mechanical stress, we observe a response in the form of an electric displacement and the strain [1]. When the intensity of external electric fields or a mechanical stresses are small, the deformations and electrical displacements caused by them are also small. Such processes are called reversible; they say that there is a linear response. The modeling of such response leads to the construction of constitutive relations in the form of linear algebraic equations, which, like the generalized Hooke law, connect external and internal parameters. We can say that the mathematical model of the behavior of such materials is described by the linear algebraic operators in which the elements of elastic, piezoelectric and dielectric constant tensors are found experimentally. The overwhelming majority of problems concerning the calculation of the physical characteristics of transducers with polarized before saturation piezoceramic elements are solved within the linear response of active materials. In this case, the complete formulation of the problem includes the equations of motion, the equations of electrostatics, geometric relationships, and the constitutive relations. The general solutions obey the corresponding initial and boundary conditions. Such problems in the mathematical plan are linear. In simple cases, they can be solved analytically, and in more complex cases, we have to use numerical methods, for example, the finite element method.

The situation changes dramatically as soon as external loads reach thresholds, and their intensity continues to increase. In this case, irreversible processes begin, expressed in the fact that the response of the material will already be nonlinear. The consequence of this is that under increasing loads we have some nonlinear equations, while for decreasing ones we have other nonlinear equations. The constitutive relations become nonlinear and ambiguous. Mathematically, they can no longer be described by simple algebraic relations, but it is necessary to use operator relations of hysteresis type.

Due to the small volume of the article, we confine ourselves of modeling irreversible polarization-depolarization processes by an electric field and mechanical stresses under isothermal processes. Irreversible processes associated with relaxation properties, with the influence of temperature, with the features of the influence of size the ferroelectric granules, the dynamics of processes, and some others will not be considered here. The main circle of questions will be connected with the analysis of existing mathematical models describing the response of the material to external influences of high intensity for the isothermal process: we will consider the principles of constructing the constitutive relations, analyze them, and formulate some conclusions.

First, we note an interesting regularity: many irreversible processes have a similar response in the sense that the relationships between external and internal parameters are mathematically described by similar relationships. For example, in plasticity media, stresses cause elastic and residual deformations; in ferromagnets, the magnetic field leads to induced and remnant magnetization; in ferroelectrics, the electric field generates induced and residual polarization, etc. In the case of cyclic processes, the responses are described by hysteresis dependences, as shown in Figures 1, 2, 3. Therefore, it is not surprising that the mathematical methods being developed for the study of certain processes are often used to describe others. Mathematical modeling of nonlinear responses of polycrystalline ferroelectric materials plays a significant role [2], therefore, the creation and use of such models is based on experimental data.

## 2. Some experimental data

The criterion for the correctness and adequacy of the work of any model is the acceptable coincidence of the predicted data with experimental data. It should be noted that a qualitative experiment is a very complex study, so most of them reflect only certain properties with simple effects. We note only those works that reflect the electric and elastic response due to the action of the electric field and mechanical stresses. Basically, these are the works where the properties of ferroelectric ceramics of the perovskite type are investigated: for example, BaTiO_{3}, or a ceramics containing lead: PZT, PLZT 8/65/35. Interesting results [3] on the response of PLZT 8/65/35 on the effect of electric and mechanical fields and similar results [4, 5, 6, 7, 8, 9, 10] for complex loads show that the loops of electrical and deformation hysteresis depend significantly on the intensity of the operating fields. Uniaxial mechanical compressive stresses along the electric field axis affect the ability of domains to rotate: the more intense the mechanical stresses, the fewer domains are able switching along the electric field as it is shown in Figures 4 and 5.

The types of considered ceramic materials are full ferroelectrics-ferroelastics, so the response of such materials to the effect of purely mechanical loads is of considerable interest. Some results of such tests can be seen in Figures 6 and 7. From this it is easy to see that under purely mechanical impacts, the “solid body”—“solid body” phase transition takes place, the material from isotropic (anisotropic) becomes anisotropic, the elastic modules of the material (tangents to the curves) change, residual strains appear that satisfy the property of incompressibility of the material (compare the values of longitudinal and transverse strains).

We note one more interesting property in the response of the material, when small loops of dielectric and strain hysteresis are investigated [9]. Small loops of dielectric and strain hysteresis due to the action of an electric field are shown in Figures 8 and 9. Small loops of strain hysteresis due to the action of mechanical stresses are shown in Figures 10, 11. Figure 10 reflects the stretching-compressing process, and Figure 11 reflects only pure compression followed by an increase in intensity.

By behavior of small hysteresis loops, one can judge the change in the elastic, dielectric, and piezoelectric modules of the material. In addition, one can see the property of incompressibility of a material, by estimating the numerical values of strains at zero stresses. Summing up, we can say that the main task of mathematical modeling of irreversible processes is the construction of hysteresis operators taking into account the changing anisotropy of material properties.

## 3. Simple models

To compare the efficiency of a particular model, it is needed to have graphic images in the form of hysteresis loops. All the graphs given later were obtained in [2], where the algorithms are also described.

### 3.1. The Rayleigh model

One of the first models describing hysteresis relationships was the Rayleigh model [11], proposed in 1887 for magnetization processes of iron. Applying it to the polarization processes of polycrystalline ferroelectric continuum, it is necessary to replace the magnetic field by electric field and magnetization by polarization. Mathematically, the branches of the dielectric hysteresis are described by parabolic relationships:

where

Conclusion: a quadratic dependence between the electric field and polarization approximates the dielectric hysteresis for small and medium intensity values of the electric field.

### 3.2. Evolutionary models

Investigating the nonlinear behavior of ferroelectric materials, the authors [12] assumed that the dielectric displacement

Under the assumption of linear functions entering into the previous relations and the linear dependence of the dipole moment on the electric field

a formula for the electric displacement was obtained

The creep functions

Conclusion: irreversible parameters were introduced in implicit form, for the find of which evolutionary laws were applied. The constitutive relations were obtained in the form of creep integrals. The creep functions are constructed under the assumption of linear dependences for functions entering into evolutionary dependencies. For the formulation of increasing and falling branches of hysteresis, inequalities for functions describing the load and response of the material were formulated.

### 3.3. Models of the theory of plasticity

In view of the similarity of the plasticity and polarization phenomena noted earlier, for qualitative description of the polarization effects, use the rheological models of the theory of plasticity. At the basis of our subsequent actions lie analogues of mechanical and electrical quantities: the generalized coordinate—the electric charge; generalized speed—current; coefficient of elastic compliance—capacity; the generalized force is the electromotive force. In the transition to continuous media, forces are replaced by mechanical stresses, displacements by strains, etc. As a result, one can write the following correspondence:

where the indices “e” and “0” indicate the induced and residual components, respectively, and

By connecting these elements in series (Figure 14a), according to the conditions

If we collect a battery of one of the chains in parallel, we see that the resulting polarization for such a connection is determined by summing the individual polarization components of each chain, and the electric field is the same for any of them. Let the chains be distributed continuously [2, 21, 22], then there is a distribution density function

Taking a battery of continuously distributed chains for a representative volume, we can determine it polarization for any value of the electric field. It is possible to perform calculations using the generalized Prager model if, in addition to the probability function of the distribution density, three functions

We will assume that the generalized Prager operator is generated by a battery of parallel connected chains with functions of the following form:

The result of the calculations is shown in Figure 16.

Conclusion: varying the parameters of the Gaussian distribution, as well as the functions included in the differential inclusions, it is possible to substantially change the shape of the hysteresis loop. This circumstance gives confidence that by the selection of these parameters a good coincidence of the modeled and experimental hysteresis loop in the one-dimensional case can be achieved.

### 3.4. The Preisach model

This model was proposed in 1935 by Preisach [23]. The model is based on the concept of hysteron [24, 25], which approximates the switching of 180° domain. Simulation of the polarization of ceramics is associated with its representation as a set of a very large number of 180° hysterons with the probability density function

defined on the plane

The distribution function in the vicinity of the coercive fields had to have a pronounced peak, which allows approximating its using known distributions, for example Gauss, with subsequent determination of the parameters entering into it. There are also distribution functions in the form of polynomials degree of intense and coercive fields [26, 27, 28, 29]. For a convenient mathematical representation, we introduce the concept of an elementary dipole hysteresis operator (Figure 21) or the “relay” operator [22] using the functions of sets. Suppose that for any fixed pair

and

In view of the finiteness of the oscillations of any continuous function between

The parameters

where

where

Various forms of loops of dielectric hysteresis can be found in [2]. For example, if we neglect the correlation coefficient in the Gaussian distribution and use the symmetry condition for the function with respect to the variable

To estimate the influence of the parameter

A Preisach model has become widely used not only in the description of magnetic and ferroelectric hysteresis, but it is intensively used in calculating the damping coefficients of many dynamical systems, including the dynamic of dipole switching [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. Two-dimensional models were also developed [49, 50, 51, 52, 53, 54, 55, 56].

Conclusion: Preisach model in the simplest case “includes” the domain structure, but operates only 180° domains. In the mathematical plan, the approximation of the real loop is carried out by elementary rectangular hysteresis loops. A model is intended only for finding the residual polarization. This model does not include mechanical stresses, which significantly reduces its practical application.

### 3.5. Model of orientation switching

A representative volume is considered as a set of

Denote the plane perpendicular to the vector

Let

For 90° domains switching, we have

A schematic representation of all domains axes

As an example [2], a hysteresis loop (Figure 25) was calculated for the next set of parameters:

Conclusion: in physical terms, this model is closer to the physics of the phenomenon of polarization of polycrystalline ferroelectrics, but does not take into account the effect of neighboring domains on each other during the polarization process. Because of this, the loop acquires angular shapes and is not suitable for describing the differential properties of the material. On the other hand, it allows us to find a value of the angle in which the directions of all the spontaneous polarization vectors are located after reaching the saturation polarization. Another drawback of the model is that it does not include mechanical stresses, which does not allow us to investigate ferroelastic phenomena.

### 3.6. Energy model of switching

In this model, a representative volume is also considered, including a set of domains oriented in space in an arbitrary manner [2] and a residual polarization vector is determined by averaging. The main difference is in describing the orientation of domains and determining the criteria for switching domains. For example, for a ferroelectric of the perovskite type, the axes of the local system are introduced

where

Conclusion: the energy model of switching at physical essence takes into account the physical phenomena of polarization of polycrystalline ferroelectrics but also does not take into account the effect of neighboring domains on each other during the switching process. The dielectric hysteresis loop has angular shapes and is not suitable for describing the differential properties of a material. In addition, this model operates only with the residual parameters of polarization and strain and does not include induced components. However, this model has an unquestionable advantage because it already includes a mechanical stresses, which allows us to describe not only ferroelectric, but also ferroelastic phenomena.

### 3.7. The Giles-Atherton model

The model is based on the so-called “limiting” (or anhysteretic) curve, derived analytically on the basis of Weiss theory and Boltzmann statistics [60, 61, 62, 63, 64, 65, 66, 67]. If in the polycrystalline ferroelectric material were no mechanisms for locking the walls of the domains, then after the removal of the electric field, the polarization of the representative volume would be zero. With great assumptions, one could say that such a material would behave like polar liquid, where the electric field tries to align dipoles along one direction, and the thermal fields impede it: impacts try to destroy the alignment pattern, and after removing the electric field, the thermal fields neutralize the polarization. Of course, a different domain switching mechanism is observed in polycrystalline ferroelectrics, but the switching model is borrowed from polar liquids.

To determine the dependence connecting the polarization of a representative volume with the strength of the effective electric field [60] let the representative volume of the ferroelectric material contain

Using the positions of the Weiss theory, we assume that the field of forces acting on the ferroelectric domain is a sum of the electric field

The potential energy of a domain that is similar to a dipole depends on the direction of its polarization vector and is expressed by the formula:

To determine the distribution of the domain axes in the presence of an effective field orienting them, it is necessary to use the Boltzmann theorem of statistical mechanics. Under the conditions of thermodynamic equilibrium, the distribution law for domains in the presence of a conservative electrostatic field differs from the law of their distribution in the absence of this field by a factor

In the one-dimensional theory, the direction of the electric field coincides with the direction of the axis

In this model, a reversible (induced) part

The reversible part is the state parameter. It can be defined as:

where

where

To solve the equation, we need use numerical methods, for example, the fourth-order Runge-Kutta method. When the electric field varies according to the harmonic law, we obtain dielectric hysteresis loops. The model includes five parameters:

Some examples of numerical calculations in which the influence of the coefficient k on the form of the hysteresis curves is studied are shown in Figures 28, 29, 30.

Conclusion: Giles-Atherton model deals with a huge number of distributed domains in a representative volume. Using of Boltzmann statistics made it possible, first, to construct an extremely possible polarization in the form of a Langevin distribution; second, to obtain the density of distribution of domains as a function of exponential type; third, to use the switching criterion. In this case, the switching is reduced to the energy rotations of the domains in dependency of the intensity of the electric field. The model contains both a reversible and a residual part of the polarization. The model contains a set of five parameters, which are chosen from the condition of coincidence of the calculated large hysteresis loop and experimental data. Due to shortcomings of the model, one can be attributed include a very rough approximation in the description of domain rotations in ferroelectrics, where this process is replaced by the process of rotation of dipoles in polar liquids in the construction of limiting polarization. This immediately affects when trying to build small hysteresis loops, where the results of calculations lead to large discrepancies.

### 3.8. Modeling of ferroelastics

It is interesting to note that, despite its shortcomings, the Giles-Atherton model can be used to describe deformation of ferroelastics [68]. Repeating the main contents of the previous model, we can introduce the limiting strain, which in the one-dimensional case of stretching-compression reduces to calculating the probability integral of the imaginary argument:

The total deformation consists of the reversible

The large and small loops of strain hysteresis, calculated from this model, can be seen in Figures 31, 32.

In [68, 69], two modifications were made for these models in order to obtain small hysteresis loops. First, to the induced component of polarization and strain was added a summand proportional to the loading field. Second, restrictions on the process of domain switching during the construction of small cycles are introduced. The results of this modernization are shown in Figures 33, 34. The obtained results are in good agreement with the experimental data shown in Figures 8 and 11, not only qualitatively, but also quantitatively.

Conclusion: the dignities and drawbacks of the model for ferroelastics are the same as those of the Giles-Atherton model.

## 4. Discussion

Each of the presented models performs the basic function related to the description of hysteresis dependencies. However, each of them is based on different prerequisites. Therefore, the results of the work of a particular model differ from the corresponding experimental data. The Rayleigh model gives a coincidence with the experimental data for small and medium fields. The Preisach model is widely used in areas where the integral properties of hysteresis are important. The model of orientation switching and model of the energy switching allow us to determine the cone of the angles all domains for which the directions of the vectors of the spontaneous polarization lined up inside this cone after the polarization to saturation, but the hysteresis loops have angular shapes. The Giles-Atherton model, taking into account the modernization, allows describing not only large but also small hysteresis loops.

Advantages and disadvantages of each of the considered models were presented earlier. Each of the models has a number of parameters that influence on the shape of the hysteresis loops. To determine these parameters, it is necessary to compare the calculated and experimental data. However, in the mathematical plan, this model belongs to the class of inverse problems and has a number of specific difficulties.

To sum up, it can be noted that the one-dimensional models allow describing hysteresis dependences of irreversible polarization processes. The more accurate the model is based on the physics of the phenomenon, the more accurate the results of its work. For example, those models that include lot of domains are more accurate in describing quantitative dependencies. However, analysis shows that at the present moment a universal mathematical model of the polarization and deformation of polycrystalline ferroelectric materials has not yet been developed, which would accurately reflect the data of the experimental dependences. Many problems that are importance remain unresolved. Unresolved problems include the following:

small loops of dielectric and strain hysteresis have been poorly investigated;

existing models do not reflect functional dependencies on the effects of any loading components with simultaneous exposure to electric fields and mechanical loads;

the question of the validity of the application of the developed models to thin films remains open;

there are no studies of how the anisotropy of the material changes with the simultaneous action of an electric field and mechanical stresses;

very poorly represented models, which take into account the effect of neighboring domains on current switching;

there are practically no studies related to the analysis of the similarity and differences in the existing models of Preisach, Rayleigh, Giles-Atherton, plasticity materials, except for [70, 71];

in fact there is no mathematical analysis of the effect of model parameters on the final result, there is no formulation of a minimum set of parameters;

very poorly represented models, which use the functions of the density of the distribution of domains.

## 5. Conclusions

One-dimensional models that reflect the hysteretic properties of polycrystalline ferroelectric media on the external effects of high-intensity electric and mechanical fields are considered. The bases of construction of each model are disassembled. The fundamentals of the construction of each of the well-known Rayleigh models, evolution models, models of plasticity theory, Preisach models, models of orientation switching, energy switching models, the Giles-Atherton model are analyzed and the results of their work in the form of hysteresis loops are presented. The main advantages and disadvantages of each model are noted; in particular, the differential and integral properties of hysteresis dependences are marked. It is noted that each of the models has a number of parameters, the choice of which affects the shape of the hysteresis loop. However, in the mathematical sense, such a choice generates additional difficulties, which are connected with the solution of inverse problems. Unresolved problems in the general problem of modeling the polarization of polycrystalline ferroelectric materials are noted. A list of works is given in which the main ideas and results of each of the earlier models are reflected.

## Acknowledgments

This work was supported by the Russian Foundation for Basic Research (Grant 17-08-00860-a).