Thermodynamically-Consistent Modeling of Ferromagnetic Hysteresis

Models of ferromagnetic hysteresis are established by following a thermodynamic approach. The class of constitutive properties is required to obey the second law, expressed by the Clausius–Duhem inequality, and the Euclidean invariance. While the second law states that the entropy production is non-negative for every admissible thermodynamic process, here the entropy production is viewed as a non-negative constitutive function. In a three-dimensional setting, the magnetic field and the magnetization are represented by invariant vectors. Next, hysteretic properties are shown to require that the entropy production is needed in an appropriate form merely to account for different behavior in the loading and the unloading portions of the loops. In the special case of a one-dimensional setting, a detailed model is determined for the magnetization function, in terms of a given susceptibility function. Starting from different initial magnetized states, hysteresis cycles are obtained by solving a nonlinear ordinary differential system. Cyclic processes with large and small amplitudes are established for materials such as soft iron.


Introduction
Hysteresis is a phenomenon relevant to various areas of science and means that the non-linear relation between two physical quantities, say input and output, changes depending on the increasing or decreasing phase of the input. In particular, ferro-magnetic hysteresis phenomena, along with the variation in magnetic susceptibility, affect the positional accuracy in magnetic resonance imaging systems [1] and occur during a typical charge-and-discharge process of a high-temperature superconducting magnet for NMR applications [2]. Hence, much effort has been devoted to the reduction and correction of magnetic hysteresis in magnetic-resonance imaging devices.
The first detailed model of hysteresis traces back to Duhem [3]. If u is a piecewise monotone input, then the output x is given bẏ where a superposed dot denotes the time derivative. Duhem-like models have been developed and investigated in several contexts, such as circuit theory [4,5] and ferromagnetic materials [6,7]. Next Preisach [8] modeled hysteresis by introducing two thresholds characteristic of the material [9,10]. Lately, further models have been developed by means of the Langevin function [11,12] and potential functions [13]. A generalization of the Preisach model was investigated in [14,15] through hysteresis operators, and a connection with thermodynamics was developed through hysteretic (clockwise and counterclockwise) potentials and a dissipation operator. Duhem-like rate equations seem to be the most convenient schemes for describing any type of hysteresis. Moreover, to our mind, once the balancing (dynamic) laws of a continuum are established, the second law of thermodynamics has to be the key point to characterizing admissible constitutive properties. Following that, in this paper, we develop a thermodynamic approach to ferromagnetic hysteresis by requiring the consistency of the constitutive functions with the second law expressed by the Clausius-Duhem inequality. Indeed, while the second law states that the entropy production, say γ, is non-negative for every admissible thermodynamic process, we follow the assumption that γ itself has to be considered as given by a constitutive function, the entropy η and the other constitutive functions. This view in essence traces back to Green and Naghdi [16], though thereafter no significant application has been developed in the literature. Lately, we have made recourse to this scheme in connection with hysteresis in plasticity [17] and ferroelectrics [18].
The purpose of this paper is to establish a model of hysteresis for ferromagnetic materials. First, general thermodynamic relations are expressed in a three-dimensional setting. The (ferromagnetic) body is allowed to be deformable, and hence, balance equations and constitutive assumptions involve mechanical and electromagnetic properties. Since hysteresis is determined by a non-linear relation between the rate of magnetization and the magnetic field, it is non-trivial to comply with the objectivity principle, whereby the constitutive equations are required to be invariant relative to Euclidean transformations. It follows that both objectivity and the balance of angular momentum hold if the magnetization and magnetic field are expressed by Lagrangian fields.
Next, with the restriction to collinear fields, we establish explicit models of hysteresis suitable for describing soft iron materials. As a thermodynamic restriction, it follows that the hysteresis curve is run in the counterclockwise sense. Examples are given of cycles with different properties of the asymptotic regime (saturation).

Notation
We consider a body occupying a time-dependent region Ω ⊂ E 3 . The motion is described by means of the function χ(X, t), providing the position vector x ∈ Ω = χ(R, t). The symbols ∇ and ∇R denote the gradient operator with respect to x ∈ Ω and X ∈ R. The function χ is assumed to be differentiable; hence, we can define the deformation gradient as F = ∇R χ, or in suffix notation, F iK = ∂ X K χ i . The invertibility of X → x = χ(X, t) is guaranteed by letting J := det F > 0. For any tensor A, we define |A| as (A · A) 1/2 . Throughout (x, t) ∈ Ω × R. We let v(x, t) be the velocity field. For any function f (x, t), we letḟ be the total time derivative;ḟ = ∂ t f + (v · ∇) f . A prime denotes the derivative of a function with respect to the argument.

Balance Equations
We consider a ferromagnetic, deformable body where electric conduction and electric polarization are negligible. Let ρ(x, t) be the mass density. The balance of mass leads to the local continuity equationρ Let T be the mechanical Cauchy stress tensor and b be the mechanical body force. The equation of motion can be written in the form where f M is the force per unit volume of magnetic character. In stationary conditions or in the approximation of a negligible electric field, we have where H is the magnetic field, M the magnetization (per unit volume), and µ 0 the permeability of free space. The balances of angular momentum and energy are taken in the form where ε is the internal energy (per unit mass), m = M/ρ, L is the velocity gradient, L ij = ∂ x j v i , q is the heat-flux vector, and r is the energy supply (per unit mass). Let η be the entropy density and θ the absolute temperature. As the second law of thermodynamics, we take the following statement: the inequality holds for any process compatible with the balance equations. The non-negative scalar γ, namely, the (rate of) entropy production per unit mass, is assumed to be given by a constitutive function. Hence, the thermodynamic process consists of η, q, r, γ, and the other functions occurring in the balance equations.
In terms of the Helmholtz free energy ψ = ε − θη the entropy (or Clausius-Duhem) inequality (4) can be written as To simplify the description of the material properties, it is understood that H and M are the fields measured in the reference locally at rest with the body.

Euclidean Invariance and Power Representation
The internal energy ε, the entropy η, and the free energy ψ are invariant under a change of frame. Hence, they can depend only on invariant quantities. A change in frame F → F * given by a Euclidean transformation, such that x → x * , is expressed by Under the transformation (6), the deformation gradient F and the magnetic field H change as vectors: and hence they are not invariant. Yet invariant scalars, vectors, and tensors occur in connection with F and H.
We first look at invariants of mechanical character. The right Cauchy-Green tensor C and the Green-St. Venant strain tensor E, defined as and the like for E. Consequently, the scalar the power T · L is apparently non-invariant. Decompose L in the classical form where D is the stretching tensor and W is the spin; we have be the second Piola stress. We observe that sinceĖ = F T DF, so The referential heat flux and temperature gradient are invariant, and so is the power: It then follows that Hence, we obtain Incidentally, For later convenience we notice that, by (9) and (10), while JT · L = T RR ·Ė + JT · W.

Consistency with the Balance of Angular Momentum
While the fields H and M M M enjoy Euclidean invariance, we now look for specific requirements induced by (2). We go back to the form (5) of the Clausius-Duhem inequality and note that −ψ + µ 0 H ·ṁ = (−ψ + µ 0 H · m)˙− µ 0 m ·Ḣ.
Hence, we let φ = ψ − µ 0 H · m and write inequality (5) in the form To fix our ideas, let θ, F, H, ∇θ be the set of variables for the functions φ, η, T, q, and γ. Computation ofφ and substitution result in The arbitrariness of∇θ,θ and L,Ḣ implies The constraint (2) results in (13) and the requirement (13) holds if ∂ F φ is related to ∂ H φ. Any fieldH of the form f (J)H is objective. Hence, we let φ depend on F through E = (F T F − 1)/2 and jointly on F and H throughH, Since where Consequently, by (14) and (15), it follows that the requirement (13) holds identically for any magnetic fieldH Owing to the form (11) of the power, the pair M M M, H seems more convenient to describe the magnetic behavior in deformable bodies. That is why we then proceed with the choice of H, i.e., f = 1, for the referential magnetic field.

Thermodynamic Restrictions
The Euclidean invariance suggests that we investigate the Clausius-Duhem inequality (5) in the Lagrangian description. Hence, we consider J times inequality (5) and use the representations (7)-(9) of the powers T · L, q · ∇θ, and µ 0 ρH ·ṁ to obtain Hereafter, we use the referential fields For later developments, it is convenient to consider the free energy Moreover, to save writing, we let By (10) and the definition of T RR , we have Consequently, Equation (16) is then rewritten to read The purpose of modeling ferromagnetic hysteresis suggests that we take (θ, F, H, M, ∇θ,Ḟ,Ḣ) as the set of independent variables, or alternativelyṀ in place ofḢ. Yet invariance requirements demand that the dependence on the derivatives occurs in an objective way. Moreover, the Euclidean invariance of the free energy φ implies that the dependence of φ R is a function of Euclidean invariants. Now, θ, E, H, M M M are invariants, and hence we let and the like for η R , T RR , q R , and γ.
Compute the time derivative of φ R and substitute in (19) to obtain The (linearity and) arbitrariness of ∇Rθ,Ë,Ḧ,θ, W implies that The symmetry condition in (21) is just the balance relation (2) of angular momentum. Hence, (20) simplifies to In the following analysis of (22), we neglect cross-coupling effects. Specifically, we assume T is independent ofḢ and ∇R θ;Ṁ M M is independent ofĖ and ∇R θ; q R is independent ofĖ andḢ. Consequently, inequality (22) splits into three sub-inequalities: The three functions γ H , γ T , and γ q are non-negative as particular cases of γ; i.e., γ H is the value of γ asĖ = 0, ∇R θ = 0 and the like for γ T and γ q . Equation (23) is investigated in the next sections; the joint occurrence ofḢ andṀ M M result in hysteretic properties of the material. As for Equation (24), the stress T RR , and hence T RR , can depend onĖ. This dependence is allowed in the form where Ξ is a positive semi-definite fourth-order tensor such that Sym → Sym. In view of (17), we have in suffix notationΞ ijRS = J −1 F iP F jQ Ξ PQRS . Hence, as must be the case, Equation (27) shows that the stress T consists of the elastic term J −1 F∂ E φ R F T , the magnetic term −µ 0 H ⊗ M, and the viscous termΞĖ. Equation (25) is the heat equation in the reference configuration. Fourier's law q = −κ(θ)∇θ is allowed so that Rate-type constitutive equations for q R are obtained by lettingq R be given by a constitutive function while q R is one of the independent variables [19].

Cyclic Processes
We first go back to inequality (19) and investigate cyclic processes of inviscid materials, Ξ = 0, in uniform temperature fields; ∇R θ = 0. In a cyclic process in the time interval Of course if T RR satisfies (26), then the corresponding integral (28) is non-negative.

Hyper-Magnetoelastic Materials
If M M M is not among the independent variables, then the arbitrariness ofḢ in (23) implies in addition to γ H = 0. The dependence of φ R on θ, E, H allows us to say that Equation (30) represents the constitutive equation of the magnetization in a hyper-magnetoelastic material.

Linear Magnetoelastic Materials
For definiteness, we look for constitutive equations associated with a special class of free energies. Let χ possibly depend on θ. Let be the magnetic susceptibility, per unit volume, in the current configuration. We assume that the free energy ρφ is the sum of a thermoelastic part ρΨ(θ, C) and a quadratic isotropic part due to magnetization: Replacing H = F −T H and multiplying by J, we find the form (32) shows that φ R is a function of invariant quantities. By (30), we have Hence, it follows that which represents the magnetization function in a linear paramagnetic material. The associated free energy in terms of M M M is obtained by a Legendre transformation of (32), Correspondingly, in the current configuration the free energy is For later convenience, we show that φ and ψ can be given by a joint dependence on H and M M M. Owing to (33), we can write (34) as and then Correspondingly,

Nonlinear Magnetoelastic Materials
According to Landau's pioneering approach [20], nonlinear isotropic paramagnets are associated with a free-energy function with a fourth-degree polynomial in the form where χ is given by the Curie-Weiss law When the applied external field H vanishes, it follows that M = 0 is a solution at any temperature. In addition, if θ < θ C , there exist infinitely many pairs of non-vanishing solutions: where e is a generic unit vector. Hence, M s (θ) can be viewed as the spontaneous magnetization modulus at θ < θ C . In the (θ, M) plane, the curve consisting of the two branches M = ±M s (θ) describes a super-critical pitchfork bifurcation which is typical of secondorder phase transitions. In addition, by letting M * = M s (0) = √ θ C /Cκ, we infer that M s is a decreasing function on (0, θ C ] so that 0 ≤ M s (θ) < M * (see the dashed line in Figure 1).  When the applied external field H does not vanish, we assume that M and H have a common direction and consider the pertinent components, M and H. Then, (37) becomes unidimensional in character and gives The corresponding curve in the (θ, M)-plane is drawn in Figure 1 (solid line) for a given positive value of H. For all θ θ C , there exists only one solution, say M 0 (θ), which approaches zero, whereas for θ θ C there are three solutions very close to solutions 0, ±M s of the homogeneous case.
Since the differential magnetic susceptibility, χ d , can be computed as the derivative with respect to H of the constitutive function for M, we infer Hence, if θ θ C , then M = M 0 (θ) is negligible and we get the Curie-Weiss law: Otherwise, when 0 < θ θ C we have M ≈ M 2 s (θ) = (θ C − θ)/Cκ so that Summarizing, the plot of χ d (θ) is given in Figure 2. Figure 1. The (Landau) differential magnetic susceptibility χ d .

Hypo-Magnetoelastic Materials
We now look at more general non-dissipative models consistent with (23). We let γ H = 0 but allow M M M to be an independent variable, whence it follows from (23) that Let n be a unit vector, n · n = 1. Any vector w can be represented as the sum of the longitudinal part, along n, and the transverse part, (1 − n ⊗ n)w, If the transverse part is undetermined, then we can write w = (w · n)n + (1 − n ⊗ n)g.
The second-order tensor M in (40) depends non-linearly on the strain E and the temperature θ, beyond M M M and H. By analogy with hypo-elastic materials ( [21], §99), we say that Equation (40) characterizes hypo-magnetoelastic materials.
Otherwise, if g = ΓḢ, theṅ A particular case follows by taking Γ such that thereby implying the vanishing of the dyadic term. Indeed, inner multiplication ofṀ M M = ΓḢ by ∂ M M M φ R and the use of (41) yield (38).

A Simple Hypo-Magnetoelastic Model
A special but significant class of hypo-magnetoelastic models is obtained assuming the free energy ψ is independent of M M M. In this case ∂ M M M φ R = −µ 0 H and In the special case Γ = 0, it follows that Otherwise, if Γ = 0, we can choose Γ =Γ(θ, C, H) such that (42) For definiteness, we exhibit a simple example assuming a quadratic expression of the free energy:

Ferromagnetic Hysteresis
Starting from the dependence on the set of variables (23) is required to hold with γ H ≥ 0. In addition, in an isothermal cyclic process, the inequality (29) has to be true. For definiteness, we now investigate hysteresis properties by letting Hence, M M M and H are subject to To select appropriate free energy functions φ R we observe that, in the hysteretic regime, M and H are neither independent nor are related in the form M = χH, as they are in the paramagnetic regime. If we assume ∂ M M M φ R = 0 and let then the requirement (45) and the representation formula (39) yielḋ This relation shows a particular case that follows by taking Γ such that where U is the Heaviside step function and α(θ)> 0 describes the possible dependence on temperature. In the material description, we have For ease of writing, we now understand that θ ∈ (0, θ C ), and hence U (θ C − θ) is omitted. Observe that Consequently, the constraint (47) holds with Γ = χJC −1 , and the representation (46) can be written in the forṁ (49)

One-Dimensional Models of Hysteresis
Assume the spatial fields H and M are collinear and the body is isotropic, or otherwise H and M are in the direction of easy magnetization (easy axis of the transversely isotropic body). We then let H = He 1 , M = Me 1 and take (e 1 , e 2 , e 3 ) be an orthonormal basis. Hence, we represent the deformation gradient in the form Inequality (50)

|Ḣ|.
Now, we consider the one-dimensional version of (48): where α = α(θ) > 0. Correspondingly, Equation (49) simplifies tȯ where M(H, θ) = (χ + µ 0 /α)H. Except at inversion points (whereḢ = 0), we have If ζ depends onḢ at most through sgnḢ, then Equation (52) is invariant under the time change t → t * = ct, c > 0, and then we say that the equation is rate-independent. As a check of consistency, we consider the limit behavior of non-dissipative materials, as is the case in some magnetic nanoparticles (see [22]). This is made formal by letting γ H = 0 and then ζ = 0 so that (52) reduces to dM dH = χ, which represents the differential susceptibility of a paramagnetic material. Let so we can write Equation (52) as a differential equation, for the unknown function M(H). The second term χ 2 sgnḢ describes hysteretic effects in that the slope changes depending on the sign ofḢ. Since χ 2 is proportional to ζ, the vanishing of the entropy production γ H results in χ 2 = 0. Hence, χ 2 = 0 is said to represent (the limit case of) hysteretic non-dissipative materials, and χ 1 represents the slope of the curve M(H) of a paramagnetic substance; possibly, the slope is not constant and depends on the values of M and H. When χ 2 = 0, we can view (54) as the positive, differential, magnetic susceptibility. We then require that Since α, ζ > 0, χ 2 satisfies according to the counterclockwise sense required by M dH ≤ 0.
In summary, the model is characterized by the paramagnetic susceptibility χ 1 = χ, the hysteretic function ζ, and possibly the temperature-dependent function α. By definition, χ 1 is fully determined by the free energy φ R , whereas χ 2 depends also on ζ. Hence, different models are obtained by using the same function φ R . The function χ 2 is connected with the entropy production through ζ, and as we will see in a while, governs the hysteretic properties of the material.
It is of interest to consider the case Hence, regardless of the form of ζ, as θ → θ C the curve M = M(H, θ C )= χH is just the magnetization curve of a paramagnetic material. As shown by (54), the hysteretic function ζ, together with α and the sign ofḢ, affects the differential susceptibility dM/dH. Indeed, dM/dH = χ 1 + χ 2 sgnḢ is the effective slope of the magnetization curve in the H-M plane, and dM/dH = χ 1 simply represents the average value of the possible slopes at a fixed point (H, M) of this plane.

Soft Iron Models
Soft magnetic materials are of interest because they are easily magnetized and demagnetized. They have low permanent magnetization (magnetic remanence) and low intrinsic coercivity, but have a high level of saturation and a high Curie temperature. To this class belong soft iron and isoperms, e.g., Fe-Ni-Cu alloys and Mn-Zn ferrites. A model for soft materials is now established within the previous scheme: Hence, we have Equation (56) is consistent with the second law of thermodynamics for a given function f and g = f − µ 0 /α. It is of interest that the constitutive relation (1.1) of [6] is similar to (56). By analogy with [6], we first consider a function g to be piecewise constant, and correspondingly, f is piecewise-linear. For definiteness, let µ 0 = 1, α > 2, and In this case, hysteresis cycles are obtained by solving the system Ḣ = ωH cos ωṫ  As θ → θ C , the parameters µ 0 /α and τ h tend to vanish so that g approaches f and (56) reduce to dM dH = f (H).
Hence, M = f (H) can be viewed as the magnetization curve of a paramagnetic material.
Some properties, e.g., counterclokwise orientation, are established in [6] by assuming that f ≥ g. This condition, which here implies µ 0 /α ≥ 0, entails that the energy expended in a complete traversal of a simple loop is non-negative ( [7] Equations (1.6) and (3.18)). However, it is not enough to ensure thermodynamic consistency with the existence of a free energy. A stronger requirement that guarantees this consistency is the existence of a positive constant > 0 such that f − g > for all H and M. In the present model, this property is trivially satisfied as f − g = µ 0 /α > 0. Unfortunately, it implies that f − g cannot vanish, not even at the limit as |H| goes to infinity, and this prevents the model from exhibiting the saturation property.
A model allowing for the saturation property can be obtained as follows. Let ζ 0 > 0 and Hence, The vanishing of χ = g(H) as |H| approaches infinity is a way of modeling the saturation property. On the other hand, this entails that f approaches µ 0 /α as |H| tends to infinity. Hysteresis paths are obtained by solving the system Ḣ = ωH cos ωṫ  As θ → θ C , the parameters µ 0 /α and τ h tend to vanish so that g approaches f and (56) reduce to dM dH = f (H).
Since we assume that g(H) vanishes as |H| approaches infinity, the same does f (H) in the limit θ → θ C . Consequently, M = f (H) can be viewed as the magnetization curve of a paramagnetic material with the saturation property.

Hysteresis Loss
Some remarks are in order on the dissipation due to hysteresis in the general scheme (51). Owing to the counterclockwise sense, MḢdt is the area enclosed in a cycle and also 1/µ 0 times the dissipation of the sample (also called hysteresis loss). For a closed curve in the H-M plane, it follows from (51) that If ζ is parameterized by temperature and strain but independent of H and M, then we can regard ζ as constant in a H-M cycle so that This is the case for the model (52), where the dissipation is proportional to the variation ∆H of the magnetic field and twice the hysteretic function ζ.
Otherwise, if ζ is given as in the soft iron model (56), then Owing to the explicit calculation of the loading and unloading curves that make up the cycle, in ( [7], Equation (3.14)) the following result is proved here, for simplicity, we assumeH = 0 for the center of the loop. Accordingly, the area of a loop of small amplitude ∆H is of order (∆H) 3 , whereas the area of the major loop is given by Since all models considered are rate-independent, the hysteresis loss is independent of the frequency at which the alternating magnetic field varies.

A Rate-Dependent Generalization
In order to jointly investigate hysteresis and frequency-dependent dissipation properties, we let where M(H, θ) = (χ + µ 0 /α)H and χ and α possibly depend on θ. Hence, (51) becomes Considering once again the one-dimensional version of (48), which represents a generalization of (56) with f = M and g = χ.
It is easy to check that this rate-type equation is rate-dependent in that the response to an AC magnetic field (i.e., a magnetic field that varies sinusoidally) depends on its frequency. For definiteness, let H(t) = H sin ωt, ω > 0. After introducing t * = ωt, we put Then,Ḣ = ωḢ * ,Ṁ(t) = ωṀ * , and assuming all parameters are constant, (57) becomes In the limit of small frequencies, ω → 0, the material behaves as a reversible paramagnet, with M = M(H, θ). Hence, χ(θ) + µ 0 /α(θ) may be considered as the static magnetic susceptibility. Otherwise, in the limit of high frequencies, ω → +∞, the ferroelectric material exhibits a frequency-independent hysteresis described bẏ Let > 0. If ζ 0 is small enough to satisfy the inequality This relation implies that M and H are not in phase under AC magnetic processes and then are related in a complex form. In addition, a dependence of ζ 0 on the derivativeḢ (and not only on its sign) would provide the same effect in the general case.

Generalization to Materials within Non-Uniform Fields
Within a quantum mechanical description, the interaction between magnetic moments is modeled by exchange integrals of the probabilistic densities ( [23], Ch. 15). The classical analogue of the interaction in non-uniform fields may be modeled by allowing dependence of the energy on the gradient of the magnetization or of the magnetic field ( [20], § 44).
For definiteness, we look for a model involving ∇R H. To account for a dependence on ∇R H, we consider the Clausius-Duhem inequality in the more general form with a possibly non-zero extra-entropy flux k R [24]. Hence, we express the Clausius-Duhem inequality as denotes the classical property that the hysteresis curve in the H − M plane is run in the counterclockwise sense. The free energy in the form (48) has the feature that, through the factor α(θ), the ferromagnetic behavior approaches the paramagnetic one as θ → θ C . Hysteresis is shown to be modeled by Equation (54), where χ 1 is the paramagnetic susceptibility and χ 2 affects the slope changes depending on the sign ofḢ. Hence, in general, hysteresis is governed by free energy φ R and a hysteretic function ζ. Two definite models have been established for the soft iron. In the first one, the free energy φ R and the hysteretic function ζ, are quadratic functions; the resulting constitutive equation is similar to Equation (1.1) of [6]. As shown by Figure 3, the saturation does not occur. In the second model, φ R and ζ are not in polynomial forms, and the saturation property was obtained (see Figure 4).
After discussing the dependence of the hysteresis loss on the quantities involved in an alternating magnetic field, some generalizations of these models were introduced: First, a rate-dependent generalization where hysteresis and frequency-dependent dissipation occur jointly. Second, a generalization to materials within non-uniform fields by allowing a dependence of the energy on the gradient of the magnetization or of the magnetic field.
A future improvement to the theory would be modeling materials where the mechanical hysteresis occurs in connection with the ferromagnetic hysteresis.