A thermodynamic approach to rate-type models in deformable ferroelectrics

The purpose of the paper is to establish vector-valued rate-type models for the hysteretic properties in deformable ferroelectrics within the framework of continuum thermodynamics. Unlike electroelasticity and piezoelectricity, in ferroelectricity both the polarization and the electric field are simultaneously independent variables so that the constitutive functions depend on both. This viewpoint is naturally related to the fact that an hysteresis loop is a closed curve in the polarization–electric field plane. For the sake of generality, the deformation of the material and the dependence on the temperature are allowed to occur. The constitutive functions are required to be consistent with the principle of objectivity and the second law of thermodynamics. Objectivity implies that the constitutive equations are form invariant within the set of Euclidean frames. Among other results, the second law requires a general property on the relation between the polarization and the electric field via a differential equation. This equation shows a dependence fully characterized by two quantities: the free energy and a function which is related to the dissipative character of the hysteresis. As a consequence, different hysteresis models may have the same free energy. Models compatible with thermodynamics are then determined by appropriate selections of the free energy and of the dissipative part. Correspondingly, major and minor hysteretic loops are plotted.


Introduction
Ferroelectrics belong to a large class of materials that exhibit a spontaneous alignment of dipoles which is ascribed to their mutual interaction. All polar crystal classes show a spontaneous polarization without applying a mechanical load and in the absence of an external electric field. Such a behaviour is due to a non-vanishing electric dipole moment associated with their unit cell and depends on temperature. Polar crystals such as poled piezoelectric ceramics exhibit both piezoelectricity and pyroelectricity. The direct piezoelectric effect is the change of polarization when applying a mechanical stress, whereas in the converse piezoelectric effect the application of an electrical field creates mechanical deformation in the crystal. Pyroelectricity can be described as the ability to generate a temporary electric potential in response to a temperature change. A recent review of mathematical theories appropriate for the description of the nonlinear coupling of electrical forces and mechanical deformation is presented in [1].
Ferroelectric materials may be strongly polarized even in the absence of any applied electric field provided the temperature lies below a characteristic temperature θ C , called the Curie temperature. Such materials are required by symmetry considerations to be also piezoelectric and pyroelectric at temperatures higher than θ C . They are thus a challenging subject for mathematical modelling in that they exhibit nonlinear behaviour and hysteresis [2,3]. Moreover, the combined properties of memory, piezoelectricity, and pyroelectricity make ferroelectrics very useful, e.g. for sensor applications [4].
However, unlike piezoelectricity and pyroelectricity, in ferroelectricity both the polarization and the electric field are simultaneously independent variables so that the constitutive functions, including free energy, depend on both.
On the basis of these ideas, the physical literature developed constitutive properties of deformable ferroelectrics by requiring that the equilibrium configuration is determined by finding the minimum free energy with respect to polarization. In turn, according to the Ginzburg-Landau-Devonshire model [5,6], the free energy is evaluated as the sum of the Landau potential, the exchange energy due to gradient terms, the dipolar term, and the electrostatic term. The Landau-Devonshire potential is a polynomial of even-order powers of the polarization P, usually limited to the sixth order. To account for the transition, the lowest order term is parameterized by the Curie temperature θ C so that the material shows hysteretic properties when the temperature is below θ C (see [7, § 69]). Many phase-field models for piezo-ferroelectric materials are based on this approach (see, for instance, [8,9] and reference therein).
Though mainly in connection with ferromagnetism and plasticity, an alternative dissipation-based approach was considered on the basis of an evolution function that should generate the hysteretic curve of the material. Borrowing from the ferromagnetic context, various one-dimensional models are investigated [10][11][12][13] where the time derivativeṖ of P is given byṖ =P(P, E,Ė), E being the electric field. A special subclass is represented by the unidimensional Duhem model P = g 1 (P, E)Ė + g 2 (P, E)|Ė|.
Many models of this kind have been proposed within the ferromagnetic context depending on the choice of the characteristic functions g 1 and g 2 (see, for instance, [11, pp.130-150]). In some special cases, the corresponding free energies (or storage functions) have been determined [14,15]. As pointed out in [14], however, for a fixed Duhem model a large class (really, a convex set) of free energy potentials is allowed. On the other hand, as we prove in this paper, for a fixed free energy potential many different Duhem models are thermodynamically compatible.
In this paper, we follow a dissipation-based approach. To our mind, the present modelling of hysteresis deserves significant improvements and the continuum thermodynamics is likely to be the natural framework where models can be established. For the sake of generality, the deformation of the material and the dependence on the temperature are allowed to occur.

Plan of the paper
The paper develops two main topics. First we set up a joint model of electroelasticity and piezoelectricity, in deformable media, by a systematic application of the second law instead of minimum properties. In this regard, we select the pertinent variables of electromagnetism in matter such as the electric field E and the electric polarization Π in the reference configuration, the electric polarization being given by a constitutive equation that involves, besides deformation and electric field, also spontaneous polarization P as a function of temperature only. Next hysteresis modelling is developed for deformable ferroelectrics, where both E and P are taken as independent variables.
In a one-dimensional setting, where the deformation and the temperature are viewed as parameters, the contemporary dependence on both E and P is related to the fact that an hysteresis loop is a closed curve surrounding a region in the (E, Π)-plane or that Π(E) would not be a single-valued function. Accordingly, along any hysteretic curve,Ṗ andĖ are not independent of each other and a connection follows from the entropy inequality. Remarkably, the corresponding relation betweenṖ andĖ consists of a twofold dependence: the former is fully characterized by the free energy and the latter is related to the dissipative character of the hysteresis. One-dimensional models compatible with thermodynamics can then be determined by appropriate selections of the free energy and of the dissipative part. For definiteness, four thermodynamically consistent models are established and, correspondingly, hysteretic loops are plotted. Each model is characterized by an anhysteretic function A (E), entering the free energy, and a hysteretic function ξ(E, P) closely related to the form of the hysteresis loops.

Balance equations and second law inequality
We consider a body occupying the time dependent region Ω ⊂ E 3 . The motion is described by means of the function χ (X, t) providing the position vector x ∈ Ω in terms of the position vector X in a reference configuration R, and the time t, so that Ω = χ(R, t) and u = x − X denotes the displacement vector. The deformation is described by means of the deformation gradient F = ∂ X χ, or in suffix notation, The symbol ∇ denotes the gradient in the current configuration Ω, whereas a superposed dot denotes time differentiation, namely · = ∂ t + (v · ∇).
Throughout, E and H denote the effective electric and magnetic fields, whereas P is the polarization vector, M is the magnetization, q and J are the free electric charge and current, respectively. Since all these fields are dependent on the frame of reference, for definiteness we let E, H and J (as well as P and M) be the fields in the moving body as measured by an observer following each material particle that is in the local frame of reference at rest with the material; all points of the frame move with the velocity v of the body at the point under consideration. Owing to the polar character of ferroelectrics, the stress T need not be symmetric.
We restrict our attention to non-conducting ferroelectrics. Inside the body, we assume the magnetization and the electric current are zero, M = J = 0, whereas the magnetic induction B and the displacement vector D are given by Their time evolution is ruled by Maxwell's equations within matter. We let ε be the internal energy density (per unit mass), T the Cauchy stress, b the mechanical body force (per unit mass), q the heat flux vector, r the external heat supply (per unit mass), ρ the mass density, and q the electric charge density. For later convenience, let π = P/ρ denote the polarization density per unit mass. The balance equations of linear momentum, angular momentum, and energy can be written in the form 1 where P denotes the convective time derivative of P [16, eq.(3.6)], also referred to as Truesdell vector rate (see, for instance, [20]). Consider the entropy inequality where θ indicates the absolute temperature, η denotes the entropy density (per unit mass) and q represents the heat flux vector. Substitution of ∇ · q − ρr from the energy equation (2) 3 gives Letting ψ = ε − θη (Helmholtz free energy), we may rewrite this inequality as follows Introducing the Gibbs free energy upon some rearrangements we can write (3) in the form We state the second law of thermodynamics by saying that either inequality (3) or (4) has to hold for any process consistent with the balance equations.
Since ferroelectrics exhibit an hysteretic behaviour at temperature below some critical value, we allow constitutive equations to involve both the electric field E and the polarization density π as independent variables. For generality, the temperature and the deformation are taken to affect the response of the material, too. Accordingly, thermodynamic restrictions may be scrutinized in the Eulerian description, by assuming as the set of independent variables (see, for instance, [17]). In this case, all thermodynamic potentials and constitutive relations must fulfil some universal principles, such as the principle of objectivity and the second law of thermodynamics.

Objective fields and rates
By the principle of objectivity (sometimes referred to as principle of material frame indifference [18,19]), the material properties of a body should not depend on the observer, no matter how he moves. More precisely, as the statement of the principle of objectivity, we assume that the constitutive equations are form invariant within the set of Euclidean frames [20]. Accordingly, they have to be invariant under any change of frame with rotation matrix Q, In particular, vectors E and π transform according to Since ∂ x x * = Q, the deformation gradient F transforms as a vector, too, It is easy to check that the Cauchy-Green tensor C = F T F is an invariant tensor, whereas F T E and F T π are invariant vectors. Other invariant vectors are F −1 E and F −1 π. In particular, we let

Remark 1
It is worth noting that E and Π are strictly related to the Lagrangian electric and polarization fields in the reference configuration, E R and P R . Indeed, we have E R := F T E = E and P R := J F −1 P = ρ R Π, respectively [7]. That is why, E and Π are considered as electromagnetic fields, in matter, rather than E and P.
Both fields E and Π are Lagrangian (or material) so that their time derivatives, are invariant vectors. On the contrary, the standard time derivative of a vector field in the Eulerian (or spatial) description is not objective, since the velocity field depends on the choice of the Euclidean frame.

Remark 2
In order to establish a constitutive rate equation for the polarization, it is preferable to use Lagrangian (or material) rather than Eulerian variables. Indeed, the former approach automatically satisfies the principle of objectivity. However, the two points of view are somewhat related as it will be shown later when rate-type constitutive equations are scrutinized. We briefly recall here some well-known objective time derivatives of the vector field w(x, t) that will be useful in the following (see, for instance, [20]): Although different invariant formulations can be used to describe constitutive relations of ferroelectrics (see, for instance, [21,22]), hereafter we choose as the set of independent field variables instead of (5). Henceforth, abusing the notation, θ indicates the absolute temperature in the material description, θ(X, t).
In order to express inequalities (3) and (4) in terms of the invariant fields E and Π, we first compute their time derivativesĖ In addition, taking into account (6) we find for any vector field w. Letting for simplicity and using (7) 1 , the second law inequality (3) takes the form After multiplying this inequality by J/θ , using (1), the Nanson's formula and the identity ρ R = Jρ, we obtain where Since E · Π = E · π and then φ = ψ − E · Π, using (7) 2 inequality (4) may be rewritten as Thermodynamic restrictions follow from these inequalities, but strongly depend on the functional dependence of the free energy. In the next sections, we will deal with two different scenarios: electroelastic and ferroelectric deformable materials.

Electroelasticity and piezoelectricity
Restricting our attention to electroelastic materials, the free energy is assumed to be independent of Π, as well asĖ andΠ, namely On the other hand, Π is determined through a constitutive function (see [7, § 25]) Moreover, ψ and Π are required to be differentiable functions of their arguments. Accordingly, and substitution into the entropy inequality (10) results in The identity allows us to write the entropy inequality in the form The arbitrariness and linearity ofθ , ∇ Rθ ,Ḟ, andĖ implies the relations and the remaining inequality reduces to Using the Gibbs free energy φ = ψ − E · Π, relations (11) become 2 Remark 3 From (12), we can say that η, Π, and T follow from the potential φ via the differentiation with respect to θ , E, and C, respectively. Indeed, (12) 4 proves that T is a symmetric tensor and (9) allows T to be written as As a consequence, as it must be in view of (2) 2 .
As a special case, let us consider the following constitutive relation for Π, where B, Γ are third-order and second-order constant tensors, respectively, whereas P is a temperaturedependent vector. Since C = 1 in the reference configuration, we regard C − 1 as the convenient measure of mechanical deformations. Hence, as E and C − 1 approach zero, the polarization vector Π reduces to P (often called spontaneous polarization). Assuming that Γ is symmetric, from (12) 2 we infer (12) and (13), we have

Piezoelectricity
To model piezoelectricity, we find it natural to refer to the reference configuration. First, observe that in the linear approximation (small deformations) Then, from (14) and (15) we have So, letting we finally obtain the small strain constitutive equations of piezoelectricity 3 which describe the direct and converse piezoelectric effects, respectively. When the applied electric field vanishes, E = 0, these relations reduce to thus showing that some polarization is produced by deformation in piezoelectric materials. On the other hand, the application of an electric field produces some stress even if S = 0, Small strain constitutive equations for dielectrics are recovered from (16) provided that the piezoelectric tensor Λ vanishes and the spontaneous polarization P is neglected, where the dyadic term is referred to as electrostrictive stress. When C and σ are positive definite, by reversing these relations we obtain where K = C −1 and Q = C −1 σ −1 (in suffix notation Q i jkl = C −1 i jkm σ −1 ml ) is called the electrostriction coefficient.

Pyroelectric effect
We observe that the vector P provides the polarization P when the applied electric field E and the stress T are zero. If this is the case, the variation in P on heating is reversed on cooling and this is consistent witḣ Since the fields P and P are collinear, this relation can be rewritten in scalar formṖ = P (θ )θ , where P can be viewed as the pyroelectric factor, that gives the change of polarization per unit of variation in temperature. Typically, P may be identified with the spontaneous polarization, and therefore, it turns out to be a decreasing function of θ which vanishes above the Curie temperature [4]. A negative value of the slope P is also in agreement with the observation that an increase in temperature increases the disorder of the dipoles and therefore tends to reduce the spontaneous polarization of the medium.

Electrocaloric effect
It is a phenomenon in which a material shows a reversible temperature change under an applied electric field. It is often considered to be the physical inverse of the pyroelectric effect. During isentropic processes, Assuming (14) and (15), then and the isentropic condition takes the form Letting C := ρ R θ∂ θ η > 0 denote the heat capacity, we finally obtaiṅ In the special case of isotropic materials, the fields P and E are collinear, and then, this relation can be rewritten in scalar formθ where P and E are the components of P and E in the common direction and P is the pyroelectric factor. The conventional electrocaloric effect is related to a negative value of P . Hence, a decrease (increase) in the strength of an externally applied electric field yields a decrease (increase) in temperature. Nevertheless, a positive slope P is also possible if P and E are antiparallel, a situation that would correspond to an inverse electrocaloric effect [23].

Deformable ferroelectric materials
Unlike electroelasticity and piezoelectricity, in ferroelectricity both the polarization and the electric field are simultaneously independent variables so that the constitutive functions depend on both of them. Hence, in order to formulate consistent models for deformable ferroelectric materials, it is necessary to assume that the free energy has a more general functional dependence than in Sect. 3: in the temperature range θ ∈ (0, θ C ] (hereinafter referred to as ferroelectric regime) ψ is also required to depend on Π, and the same dependence also applies to φ. On the other hand, it is well known that all the ferroelectric materials exhibit piezoelectric effect due to lack of symmetry [2].
With the aim of recovering the results of the previous section as a special case, we also suppose that in the ferroelectric regime Π is partially determined through a constitutive function, namely where P is an independent variable that takes into account the residual polarization density. The adjective "residual" is justified by assuming that Ξ vanishes when C = 1 (no deformation) and E = 0. Under stationary conditions, Π reduces to the so-called spontaneous polarization, and therefore, it is expected to disappear in the paraelectric regime in which θ > θ C . Hereafter, abusing the notation, we write Since ψ and Ξ are required to be differentiable functions of all their arguments in the ferroelectric regime, it followsψ Substitution into the entropy inequality (10) results in Multiplying by θ/ρ R > 0 and using the identity we eventually can write the entropy inequality in the form As in Sect. 2, the arbitrariness and the linearity ofθ, ∇ Rθ , andḞ imply the relations The remaining inequality reads By virtue of (19) 2 , two separate inequalities follow, the customary heat conduction inequality q R · ∇ R θ = J q · ∇θ ≤ 0 and the reduced inequality involvingĖ andṖ, which is proper to the ferroelectric regime. Using the Gibbs free energy φ = ψ − E · Π, relations (19) and (20) Following the same reasoning as in Remark 3, even in this case the Cauchy stress tensor can be represented by means of the formula (13).

Remark 4
Although the values E and P are taken as independent variables within the ferroelectric regime, the same need not be true for their objective rates, as a consequence of (21). Indeed, under the assumption thaṫ E andṖ are independent of each other and may be arbitrarily chosen, (21) holds if and only if When φ is a C 2 function, these equalities are contradictory in that we may infer This suggests thatĖ andṖ have to be linearly related to each other when the free energy φ depends on both E and P.
Remark 5 Since π = FΠ, in the spatial description (18) must be replaced by where P P P = FP. Then, exploiting the identities (8), the reduced inequality (21) may be rewritten by involving the Cotter-Rivlin rate of E and the Oldroyd rate of P P P, namely ∂ Eφ + π · E + ∂ P P Pφ , · P P P ≤ 0, whereφ(F, θ, E, P P P) = φ(F T F, θ, F T E, F −1 P P P).

Remark 6
The occurrence of the objective derivatives P in (2.2) and E, P P P in (22) is not a subjective choice.
Rather, P follows from the balance equation (see [16]). Instead, E and P P P are a consequence of the selection of the vectors E = F T E and Π = F −1 P/ρ as appropriate variables to describe the properties of the material.

Rate-type models
Many different vectorial schemes are compatible with inequality (21), most of which can be gathered in the class of rate-type models. Indeed, such inequality may be satisfied assuming that the l.h.s. is balanced by a nonnegative-valued dissipation function d, namely The dissipation function d may depend not only on all the independent variables, but also on their time derivatives. A similar equation can be obtained by involving the free energy ψ instead of φ, A particular case follows by choosing d = ξ Ė n , where n ∈ N, ξ is a nonnegative-valued function (independent of the rates) and · denotes the vector norm 5 . Accordingly, may be regarded as a rate-type constitutive relation for the residual polarization of a ferroelectric material. If n = 1, this differential equation is invariant under the time transformation t → c t, c > 0, and hence, any associated model is rate independent. Relation (23) becomes unidimensional in character when P and E have the same common direction, and leads to the well-known class of Duhem models (see, for instance, [11,17]) provided that n = 1 and ∂ P φ does not vanish. Quite similar models follow by choosing a different dissipation function, for instance, by letting where n ∈ N and ν is a nonnegative-valued function. At any time whenṖ does not vanish, this equality may be rewritten involving the free energy ψ as In particular, it may be simply satisfied by letting Remark 7 In the case n = 1, similar relationships (where H and M take the place of E and P) provide the basis for a vectorial energy-based model successfully applied to the incremental formulation of ferromagnetic hysteresis [24]. When n = 2, it readsṖ = 1 ν (E − ∂ P ψ) := − 1 ν ∂ P φ, looking like the kinetic equation of some phase-field theory [8,9,25,26]. Indeed, after properly choosing φ, this simplified approach is close to the Ginzburg-Landau-Devonshire theory in the special case where the free energy is independent of the gradient of polarization [5].
We finally stress that similar relations can be obtained in the spatial description starting from (22) rather than (21).

A general scheme for the hysteretic behaviour
More complexṖ -Ė dependencies (even those involving memory) are compatible with inequality (21) (see, for instance, [13]). Nevertheless, in order to obtain a simple scheme for an exhaustive description of hysteretic processes in the ferroelectric regime, here we restrict attention to the evolution equation (23) with n = 1, namely Since the hysteretic behaviour cannot take place when ξ = 0 and the dissipation function vanishes, ξ is referred to as hysteretic function. Although the major hysteresis cycle can even be described by a phase-field approach (see, for instance, [27,28]), the scheme devised here has the advantage of being able to model even minor loops, as will be highlighted by means of examples. Some assumptions are now made in order to ease the comparison with electrolasticity. First, we assume the following splitting Ξ (C, θ, E) = Ξ 1 (C, θ) + Ξ 2 (θ, E). (25) As a consequence, we can write and then so that (24) takes the form As a special case, consider the linear constitutive relation for Π in the form As already mentioned with regard to (18), this ensures that Π reduces to P when C = 1 and E = 0. On the other hand, the dependence of ψ and Π on P is lost in the paraelectric regime, when θ > θ C . In that regime, all results of Sect. 3 continue to hold true, while specifying that Ξ 1 , Ξ 2 are not necessarily continuous across θ C . Under the assumption of small deformations, a linearization procedure such as that performed in Sect. 3.1 provides the stress charge equations where S := 1 2 (∇u +(∇u) T ). Unlike (16), however, here P = ρ R P is not a function of θ , but is an independent variable whose evolution equation is obtained from (24) by taking its small strain approximation. This can be achieved by multiplying it by ρ R and then replacing E, P and Π with E, P and P/ρ R , respectively, that is, where Φ = ρ R φ(θ, S, E, P). We report that constitutive assumptions (27) are similar to eqns. (13) and (25) in [9] where P is referred to as ferroelectric polarization vector, a vector-valued order parameter whose evolution is ruled a priori by a Ginzburg-Landau equation.
Finally, a spontaneous deformation may be introduced as a function of P (see, for instance [9]), so that the linearized strain tensor splits into the sum In this case, both polarization P and stress T depend on the elastic part S el , only. Then, after replacing S with S − S t in (27), we obtain A full analysis of this model is, however, beyond the scope of this article.

Unidimensional models of ferroelectric hysteresis
The fields E and P (as well as Π) are assumed here to be collinear with a given fixed unitary vector a. Let E and P be the components along the common direction. Hence φ turns out to be a function of E, P, parameterized by θ and C, on which the dependence is understood and not written. From (26) and (25) it follows Ξ 1 and Ξ 2 being given constitutive functions with Ξ 2 (0) = 0. As a result, Eq. (24) becomes unidimensional in character and reads where By integrating (30) along any closed curve in the (E, Π)-plane, we have This inequality implies that the curve is run in the counterclockwise sense. Looking at time intervals whereĖ = 0, upon division byĖ eq. (30) becomeṡ Except at stationary points whereĖ = 0,Ṗ/Ė equals the derivative of P with respect to E, dP/dE. Let The residual polarization P is then given by the differential equation We are now looking for explicit representations of χ 1 and χ 2 separately. While χ 1 is fully determined by the potential ψ(E, P), χ 2 depends also on ξ , and hence, many different models can be constructed by using the same energy potential. Although more complex dependencies of ξ are compatible with inequality (30), in order to obtain a simple but fairly general description here we let ξ = ξ(E, P) ≥ 0.
In relation to the wide literature on ferroelectrics, we observe that (30) involves merely a potential function φ and a hysteretic function ξ . Here we have no recourse to switching surfaces analogous to yield surfaces in plasticity (see e.g. [8] and refs therein) or regions of reversible behaviour (see e.g. [3,29] and refs therein).

The anhysteretic curve
Before looking for the expressions of χ 1 and χ 2 and then determining the model of ferroelectric materials, we start from assuming the existence of a characteristic non-decreasing odd function A such that the dependence on the parameter (temperature) θ being understood and not written. Taking into account that the hysteretic function ξ describes the dissipative and hysteretic properties of the material, we refer to A as anhysteretic function [17] and to the corresponding graph A in the (E, Π)-plane as anhysteretic curve. Indeed, by (30) the thermodynamic condition corresponding to (34) is Hence, on the anhysteretic curve A is likely to represent the non-hysteretic behaviour, where χ 1 can be viewed as the anhysteretic residual susceptibility. This in turn suggests that we let χ 1 be a function of E only (possibly parameterized by C and θ ), namely Finally, the characteristic property (34) and the required counterclockwise property (31) of ferroelectric cycles suggest that (36) Remark 8 By virtue of (29) and (33), it follows dΠ If ξ = 0, and hence χ 2 ≡ 0, (37) simplifies and describes a family of models with paraelectric behaviour where χ represents the anhysteretic differential susceptibility. Accordingly, we require that it must be positive as well as the slope of each branch of the hysteresis curve generated by (37), namely In particular, the former condition implies Unlike the counterclockwise property (36), inequalities (38) are not a consequence of the second law; they are assumed as physically sound restrictions so that the polarization on each branch of the hysteresis curve be a monotone function of the electric field, as observed.

The free energy potential
The expression of χ 1 is given by (35), but at the same time it is related to ψ by (32) 1 . Therefore, it is first necessary to establish the form of ψ compatible with these requirements. To this end, we now select L, G, and F being differentiable functions to be determined. Upon substitution of ∂ P ψ and ∂ E ψ into (32), we obtain which implies In order to determine L, we observe that condition (36) reads Since ξ ≥ 0, it follows which is easily satisfied for any β > 0 by letting Summarizing all these results, we conclude that 1. the free energy ψ takes the form where β is a positive constant possibly dependent on θ and

the anhysteretic function depends on χ as follows
According to (29), we get

Remark 9
In order that all hysteresis loops shrink to the anhysteretic curve as θ approaches θ C from below, we require that 1/β vanishes as θ → θ − C , for instance, by letting where h is the positive part of 1 − θ/θ C , namely Taking advantage of this remark, from results (1) and (2) we obtain so that both implicitly depend on θ . For any E, χ(E) tends to A (E) as θ → θ C . Moreover, the hysteretic effects expressed by χ 2 become smaller and smaller as θ approaches θ C from below and completely vanish in the paraelectric regime. As a final comment, the anhysteretic function A , the positive parameter α, and the hysteretic function ξ completely characterize the differential model (37) of the material.

Remark 10
Taking into account (38), A and ξ must verify the following monotonicity conditions

The hysteretic function ξ
To determine an explicit representation of χ 2 , we first observe that the only thermodynamic requirement on ξ is the non-negative valuedness. In addition, we require that it fulfils monotonicity conditions (42). With this in mind, we now examine the hysteretic properties associated with a given φ and four different functions ξ .

Model
This model is rather rough but nevertheless describes both bilinear cycles and the saturation of Π . It is able to reveal that a suitable choice of ξ can introduce a saturation threshold analogous to the yielding threshold in plasticity. We confine our attention to the strip |Π | ≤ γ where ξ takes nonnegative values and satisfies (42) provided that ξ 0 ≤ χα/ h. From (37) and (41), it follows that Hence, the hysteretic behaviour is modelled via A and the positive parameters α, ξ 0 , γ . The saturation property of the model can be characterized by the vanishing of χ − |χ 2 | (see also [14]), namely If this is the case, the anhysteretic curve A is given by all pairs (E, P) such that E ∈ R and The saturation property is fulfilled for any value of θ in (0, θ C ) provided that ξ 0 depends on E and h as follows = ω E cos ωṫ where h ∈ (0, 1) depends on the fixed value of θ ∈ (0, θ C ). In Fig. 1, cyclic processes with large amplitudes (major loops reaching saturation) and small amplitudes (minor loops) are depicted assuming A (E) = k E, k = γ = 3 , α = 5/12, h = 5/6, ξ 0 = 1/2.

Model II
From (37) and (41), it follows that This is a very simple model: it is characterized by the positive parameters ξ 0 , α, and the anhysteretic function A , only. Condition (42) is fulfilled within the region Assuming a fixed temperature θ < θ C , it is of interest that (43) has some analogies with a differential equation widely investigated in the literature for the model of hysteresis in ferromagnetic materials [30], The correspondence is made formal via the identifications The advantage of the present scheme, however, is to give evidence to the compatibility of the differential equation (44) with thermodynamics; the compatibility is guaranteed if the functions f and the anhysteretic differential susceptibility g are subject to the relation where h/α is an arbitrary positive parameter that depends only on θ . As a particular case, compatibility holds if g is piecewise constant and f is correspondingly piecewise linear and increasing, as devised in [10,31]. Unfortunately, even if g → 0 as |E| → ∞, this relation prevents the model to exhibit the saturation property unless τ → 0 as |E| → ∞. This remark suggests a generalization of (44) by letting τ be replaced by τ g itself (see Model IV).
Hysteresis cycles in the (E, Π)-plane are obtained by solving the system Fig. 2, cyclic processes with different amplitudes are depicted assuming and h/α = 0.2, ξ 0 = 2.8. In this case, g = χ(E) = 1.3[1 − tanh 2 (1.3 E)] + 0.1. The solution starting from E 0 = Π 0 = 0 with increasing values of E(t) is referred to as initial polarization curve. Its first and second derivatives at zero are given by The latter represents the curvature at zero of the initial polarization curve and is positive (as experience indicates) provided that A (0) + ξ 0 h 2 /α 2 > 0, a condition that certainly occurred if A (0) = 0 as in example (46) (see the details in Fig. 2).

Model III
We address this model to describe the asymptotic saturation of Π for large values of |E|. Monotonicity property (42) is met provided that ξ 0 ≤ α/ h. From (37) and (41), it follows that Hence, the hysteretic behaviour is modelled via A and the positive parameters α, ξ 0 , only. Hysteresis cycles in the (E, Π)-plane are generated by the system Ė = ω E cos ωṫ where h ∈ (0, 1) depends on the fixed value of θ ∈ (0, θ C ). In Fig. 3, cyclic processes with large and small amplitudes are depicted assuming and h/α = 0.4, ξ 0 = 0.2. In this case, g = χ(E) = 0.8[1 − tanh 2 (0.9 E)] and the saturation property is fullfilled.

Model IV
As in Model II, also here the vanishing of χ as |E| approaches infinity is a way of modelling the saturation property of hysteresis. Unlike that case, however, here no threshold is set for polarization, but saturation is reached only asymptotically. Monotonicity property (42) is met provided that Hysteresis cycles are obtained by solving the system starting from different initial states (E 0 , Π 0 ). For instance, in Fig. 4 cyclic processes with large and small amplitudes are depicted assuming and h/α = 0.4, ξ 0 = 2.5 (τ h = 1). In this case, g = χ(E) = 0.9[1 − tanh 2 (0.9 E)] and the saturation property is fullfilled. The first and second derivatives at zero of the initial polarization curve are given by In summary, all models describe rate-independent hysteresis in soft ferroelectric materials. Models I and III involve hysteretic functions ξ which are proportional to |Π −A (E)|. In particular, Model I shows 'perfectly bilinear' hysteresis loops in connection with a piecewise linear anhysteretic function A and constant values of χ and χ 2 . When χ = |χ 2 |, it accounts for the appearance of the saturation regime. In Model III, the anhysteretic function involves a hyperbolic tangent. Monotonicity condition and saturation property are satisfied throughout the (E, Π)-plane. In Models II and IV, the hysteretic function ξ is proportional to [Π − A (E)] 2 and the anhysteretic function involves a hyperbolic tangent. Both satisfy the condition of monotonicity only in a strip around the anhysteretic curve, but, unlike Model II, Model IV also enjoys the saturation condition.

Conclusions
The purpose of this paper was to establish a systematic constitutive scheme, for deformable ferroelectric materials, compatible with thermodynamics. In the first part of the paper, the well-known model of linear piezoelectricity is rewritten by using invariant form of the electric and polarization fields, E and Π. As is customary, Π is determined by a constitutive equation that depends on deformation, electric field, and temperature (see (14)). To describe deformable ferroelectrics, we assume a similar relation, (18), where the residual (or spontaneous) polarization, P, is an independent variable. Indeed, by the natural view that an hysteresis loop in the unidimensional case is a closed curve surrounding a region in the (E, Π)-plane, the novel starting point of this paper is a set of constitutive equations where the electric field E, the residual polarization P, and their time derivativesĖ,Ṗ are simultaneously independent variables. These constitutive functions are consistent with the second law of thermodynamics in that satisfy the entropy inequality. Among other results, the standard balance of torques is a consequence of thermodynamic restrictions. In addition, the second law implies a general inequality involvingṖ andĖ so that, along any hysteretic curve, a connection betweenṖ andĖ follows.
In the unidimensional case, two features have to emerge from this paper. First, the thermodynamic restriction, which is expressed by an inequality, makes the relation betweenṖ andĖ to consist of two parts. One is fully characterized by the free energy. The other one is related to the dissipative character of the hysteresis. Models compatible with thermodynamics can then be determined by appropriate selections of the free energy and of the dissipative part. Second, once physical reasons indicate a dependence on E of the reversible differential susceptibility as in (35), the present scheme provides the corresponding free energy φ as shown in Sect. 5.2. In particular, see (39), φ turns out to be represented as a quadratic function of P. It is of interest that while the anhysteretic function is fully determined by the free energy potential, the dissipative part of the model, characterized by χ 2 , involves both the free energy influence, via the denominator ∂ P φ, and the non-negative purely thermodynamic term ξ . For any free energy φ then, there are a family of dissipative parts characterized by the chosen function ξ .
Relative to known models of hysteresis, Model II looks very similar to the Coleman-Hodgdon model that describes rate-independent hysteresis in soft ferromagnetic materials [10,31]. Despite of this similarity, we stress that the inequality f − g ≥ 0 (that is, the requirement that the hysteresis loss associated with a primitive hysteresis loop be not negative, see [10, eqn (2.18)]) is not enough to ensure thermodynamic compatibility. Indeed, we proved here that a stronger condition is required in order to build up a compatible free energy potential ψ, namely f − g = for some = 1/α > 0 (see Eq.(45)). In particular, this thermodynamic condition implies the inequality +∞ |E| [ f (y) − g(y)]e −τ y dy > 0 foranyE ∈ R that is stronger than (1.5) in [10], where the vanishing of the integral is allowed.