Exact formulation for hysteresis loops and energy loss in Stoner–Wohlfarth systems Open Access

Magnetic materials exhibit a broad variety of compositions and microstructures and a wide spectrum of responses to a varying magnetizing field $H⃗$⁠.^1–3 

Although such behavior is the result of three-dimensional (3D) magnetization process, a two-dimensional theoretical approach is often regarded as appropriate, in addition to being perfectly adequate for a large class of materials, such as steel sheets, amorphous and nanocrystalline alloys, and thin films. It is also possible, in many cases, to consider a given material as an ensemble of mesoscopic regions (“particles”), either interacting or independent, each characterized by the local polarization $J⃗s$ and a magnetic anisotropy {uniaxial [see Fig. 1(a)], biaxial, cubic, etc.}, and even single-domain, if sufficiently small.

The Stoner–Wohlfarth (SW) theory,^4–8 in its native formulation (which considers single-domain particles as independent, endowed with uniaxial anisotropy, and at T = 0 K), fits the simplifications mentioned above and is, thus, widely exploited in modeling. The main fields where the SW approach (often with some modifications) plays an important role are (1) magnetic nanoparticles,^9–12 notably for magnetic hyperthermia;^13,14 (2) ferromagnetic thin films^15–17 and magnetic random access memories;^18–20 (3) amorphous and nanocrystalline materials;^21–25 (4) soft magnetic composites;^26,27 (5) permanent magnets;²⁸ and (6) study of the effect of spin-polarized currents.^29,30

In the SW model, the response to an applied field $H⃗$ of a particle characterized by the local uniaxial anisotropy $K↔$ is described by a vector hysteresis transducer providing the constitutive law for the reversible coherent rotations and the irreversible jumps of the local polarization $Js⃗$⁠. In other words, one gets the equilibrium orientation γ of each $J⃗s$ as a function of the anisotropy constant K, its orientation φ, and $H⃗$ [see Fig. 1(a)]. With $J⃗$ being the polarization of the whole particle assembly, the $J⃗(H⃗)$ constitutive equation of the system comes out after integration over the probability density functions (pdfs) of K and φ. The design of magnetic components needs to implement $J⃗(H⃗)$ in numeric modeling codes: a step with high computational cost, scaling up with the device size. Indeed, an explicit γ = γ_H(K, φ) solution for the $J⃗s$ orientation is not available so far, being the equilibrium angle obtained either following numeric procedures^25,31 or using analytic approximations.

In the first case, the problem is commonly faced by adopting the astroid representation for the SW transducer [see Fig. 1(b)] and exploiting the so-called “tangent method,”⁷ as well as searching for the minimum of the Gibbs free energy of each particle utilizing, e.g., the Newton–Raphson method. These are time-consuming approaches, not efficient when implemented in codes used to design magnetic devices.

Researchers following the second way tried, for example, to recover the SW loop branches with ad hoc analytic functions,³² or proposed an H = H(J) “inverse” constitutive law,¹² or limited the analytic solution to the low-field regime.³³ Alternatively, an SW-like hysteron with different kinds of approximations for the anisotropy energy term,^34,35 or with an equivalent field playing the role of the anisotropy,³⁶ has been also adopted.

This work is not focused on the search of the γ = γ_H(K, φ) single particle solution [see Fig. 1(a)], afterward integrated over the whole system, but it deals at the start from the entire particle assembly by handling each quantity as a random variable (rv) with associated pdf. Statistical methods^37,38 are then utilized to retrieve p_H(γ): the equilibrium angle pdf, ending with an exact expression for the hysteresis loops (and consequently for the energy losses) of the whole system. The single particle behavior comes out as a particular case. Apart from its importance from the basic viewpoint, this result turns out to be useful for applications (chiefly in the magnetic loss prediction³⁹), shortening the cumbersome and slow procedure needed for the design of devices.

In this paper, after a review of the SW model features essential for our purposes (Sec. II A), and some hints about the numerical approach and the adopted statistics (Sec. II B), the core of the work is discussed in Sec. III. Here, the statistical approach is outlined, and the pdfs controlling the equilibrium orientation of the local polarization vs a quasi-static $H⃗$ are worked out, separately accounting for the irreversible (Sec. III A) and reversible (Sec. III B) phenomena. In both these subsections, the pdf evolution is described starting from the demagnetized state and following the system magnetization along the branches of the hysteresis loop driven by an alternating field [Fig. 1(c)], separately for the four Q quadrants shown in Fig. 1(d). Eventually, the obtained constitutive equation (and thus the hysteresis loop and loss) is successfully tested against the numerical procedure, both for the entire assembly (Sec. III B 5) and for the case of a single particle (Sec. III B 6). An approximate analytic expression for the energy loss, not involving the knowledge of the hysteresis loop, is worked out and effectively checked (Sec. IV), and finally, in Sec. V, some possible advancements of the proposed approach are listed.

The introduction reports the strategies followed in the literature to find the γ equilibrium angle for each particle. In the following, the statistics of the system is described by the ψ(K, φ) joint probability density function, with K > 0 and φ limited to the [0;π/2] range for symmetry reasons.

In this work, we have abandoned the problem of finding the exact analytical solution for the local $J⃗s$ equilibrium orientation γ of a single SW particle in favor of a statistical approach that considers the whole assembly of SW particles from the beginning. Accordingly, the quantities K, φ, and γ (−π to +π) are treated as random variables, along with their ψ(K, φ) and $pH(b,Q)(γ)$ pdfs (the first normalized to 1 for each φ angle). In the latter, “b” indicates the loop branch considered [see Fig. 1(c)], “Q” is the quadrant to which γ pertains [see Fig. 1(d)], and H plays the role of a parameter. All along the paper, the evolution of the γ values, which are bound to increase from −π to π [see Fig. 1(d)], dictates the sequence adopted to discuss the Q quadrants (III → IV → I → II).

Along each loop branch, the applied field triggers the irreversible rotations (IRs) and drives any local $J⃗s$ of the whole particle ensemble to its equilibrium angle (reversible rotations) at the same time. Despite the ensuing complex interplay between reversible and irreversible processes, the technique developed here treats them separately, and the γ pdf is consequently worked out in two steps: before finding $pH,IRR(b,Q)(K,φ)$⁠, i.e., the pdf determined by the IRs only, which bring the local $J⃗s$ to their absolute minimum (Sec. III A), and subsequently accounting for the coherent reversible rotations (by means of the $J⃗s$ relaxation to the equilibrium orientation), ending with $pH(b,Q)(γ)$ (Sec. III B). This approach is sketched in Fig. 2.

Eventually, after the integration of $pH(b,Q)(γ)$cos γ over the [−π; +π] γ domain, one gets the constitutive law of the system [Sec. III B 5, Eq. (28)].

At this stage, we assume $H⃗$ to drive the irreversible processes only, again leaving the local $J⃗s$ to point along one of the two sides of the φ easy axis directions of the particle to which it pertains. For this reason, the magnetization distribution among the quadrants, rearranged by the occurrence of the IRs, is again described by a pdf stated in terms of K and φ instead of γ: $pH,IRR(b,Q)(K,φ)$⁠, for the three loop branches of Fig. 1(c), as listed in Sections III A 1 -III A 4, and illustrated in Fig. 3 for a given φ = φ₁.

For H increasing from DEM [Fig. 1(c)], Eqs. (1b) and (1c) state that the $J⃗s$ in quadrant III smoothly (i.e., without IRs) crosses the γ = −π/2 threshold and move to quadrant IV when φ > π/4 and H > H_III→IV = (K/J_s)sin 2φ [a value always lower or equal to H_c,K(φ): the threshold for IRs]. As stated above, we are not interested at this stage in the equilibrium position of $J⃗s$⁠, so we account for these $J⃗s$ reversible rotations by assigning the same pdfs to these quadrants: $pH,IRR(FMC,IV)=pH,IRR(FMC,III)$⁠. Note that this “doubling” of the $pH,IRR(FMC,Q)$ does not affect the normalization because each $J⃗s$ does not pertain to Q = III and Q = IV for the same field, the Step 2 (Sec. III B) taking charge of finding the actual $J⃗s$ orientation vs H. As H increases from zero [where K_c,0 = 0, see Eq. (4)] to the peak value H_p, the $J⃗s$ in quadrant III or IV switch to quadrant I if K < K_c,H(φ), and the corresponding $pH,IRR(FMC,Q)(K,φ)$ is reported in Fig. 3(b). Note that $pH,IRR(FMC,II)$ = 0 because the system geometry [Fig. 1(a)] does not allow irreversible rotations to the second quadrant, a fact in agreement with the detailed discussion reported in the  Appendix, which shows that no energy minima of u_H [Eq. (1a)] can appear in Q = II.

When the field is reduced from H_p to zero, the polarization decreases from its peak value J_p = J(H_p) to the remanence J_r. Along this branch, no IRs occur, and thus, $pH,IRR(REC,Q)(K,φ)$ is frozen to the value reached along FMC at H = H_p: $pHp,IRR(FMC,Q)(K,φ)$ [see Fig. 3(c)]. Here, again, we have $pH,IRR(REC,IV)$ = $pH,IRR(REC,III)$ because the “doubling” of the pdfs, pointed out in FMC persists, but, as above, it does not affect the normalization.

Eventually, the system is brought from the remanence to the negative tip of the hysteresis loop, lowering H from 0 to −H_p, with the magnetization evolving as described by $pH,IRR(NDB,Q)(K,φ)$ [Fig. 3(d)]. For what concern the IRs and the smooth $J⃗s$ transitions discussed in Sec. III A 1, a “symmetrical” reasoning applies here, with the quadrants I, II, and IV correspondingly playing the role of Q = III, IV, and I in FMC, and the γ = π/2 crossing occurring at the threshold H_I→II = −H_III→IV [≤H_c,K(φ)], when φ > π/4.

The quantities K_c,H(φ₁) and K_c,Hp(φ₁) of Fig. 3, signaling a variation of the $pH,IRR(b,Q)$ values for a given φ = φ₁, become the K_c,H(φ) and K_c,Hp(φ) threshold functions when considering the entire 0 ≤ φ ≤ π/2 range. As H, starting from the demagnetized state (DEM), oscillates between the ±H_p peak values of the hysteresis loop, the K_c,H(φ) threshold follows it, swinging between K_c,0 = 0 and K_c,Hp(φ), with the $pH,IRR(b,Q)(K,φ)$ function that accordingly changes. An example for Q = I across NDB is depicted in the (φ, K) plane of Fig. 4.

Given H and γ, and considering the scheme of Fig. 2, for the three branches of the hysteresis loop indicated in Fig. 1(c), the $pH(b,Q)(γ)$ pdf, which accounts for irreversible and reversible processes, is obtained through standard statistical methods.^37,38

According to Sec. III A 4 and Fig. 3(b), along this branch, the K_c,Hp(φ) curves do not play any role (and thus are not depicted in Fig. 5), whereas the K_c,H(φ) ones not only signal the K thresholds where the $pH,IRR(FMC,Q)$ pdfs [retrieved from Fig. 3(b)] change but also represent a portion of the D_H integration domains (magenta regions in Fig. 5, for Q = III and Q = IV) border.

a) $−π≤γ≤−π/2(Q=III)̄$

b) $−π/2≤γ≤0(Q=IV)$

c) $0≤γ≤π/2(Q=I)$

d) $π/2≤γ≤π(Q=II)$

Figure 7 shows, for the whole −π ≤ γ ≤ π range, an example of $pH(b,Q)(γ)$ vs H evolution, calculated with Eq. (26), and the corresponding hysteresis loop [Eq. (28)] with J_s = 1.61 T. We have assumed ψ(K, φ) =f(K)g(φ), with the marginal pdfs supplied by Eqs. (5a) and (5b) and $K$ = 3000 J/m³.

The effectiveness of the statistical method is apparent in Fig. 8, where the quasi-static hysteresis loops, computed at different peak polarization J_p using Eq. (28) and the corresponding loss figures (inset), perfectly recover the ones obtained via the numerical procedure presented in Sec. II B [again it is assumed ψ(K, φ) = f(K)g(φ), with the marginal pdfs furnished by Eqs. (5a) and (5b)].

The ability of the analytic procedure to reproduce the elemental loops obtained when H oscillates with peak values |H_p| = H_c,K(φ₁) [Eq. (3)] and J = J_s cos γ measured along the $H⃗$ direction is apparent in Fig. 9, from the comparison with the cycles numerically computed.

The developed mathematical tool analytically predicts the response of an assembly of non-interacting SW particles driven by a quasi-static alternating magnetic field. The exact integral relationships for the hysteresis loops and losses allow one to bypass the computationally slow numeric procedure necessary to implement the SW model in software used to predict the behavior of devices containing magnetic components.

The versatility of this approach lies in its skills listed below.

1. It works with any anisotropy distribution, permitting one to easily predict the role played by the sample texture (e.g., the effect of a macroscopic easy axis induced either by annealing or by applied stress).

2. The same analytic-statistical procedure, although more complicated, can be exploited with any field history (e.g., rotating fields,^31,34 loops with bias, and asymmetrical minor loops).

3. Marginal modifications of the approach allow it to account for the presence of domain walls inside the particles, following the idea suggested in Ref. 22.

4. The whole mathematical structure is prone to a generalization for systems of particles endowed with biaxial anisotropy.

5. Under whatever field history and accounting for interactions, in Ref. 40, a mathematical technique capable to track the irreversible switches in a particle assembly described by the SW model has been worked out (so preserving the non-local memory of the system), together with a graphic representation that turns out to be an extension of the one operating in the Preisach model. One can then envisage a coupling between such a tool and the one developed here, with the aim to assess a more general analytical approach to the particle systems description.

This research work was carried out in the framework of the Grant No. 19ENG06 HEFMAG project, which was funded by the EMPIR program, and co-financed by the Participating States and the European Union’s Horizon 2020 research and innovation program.

The author has no conflicts to disclose.

C. Appino: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.

For a general γ′ equilibrium orientation of $J⃗s$⁠, with $γlow′≤γ′≤γ$ (where γ depends on the Q quadrant considered, and $γlow′$ is its lower limit), the $Δφ*(γ′)$ range, corresponding to a minimum of the particle Gibbs energy u_H(K, φ; γ′) [Eq. (1a)], is determined by the overlap between two φ regions: the first one, where $KH*(φ;γ′)$ [Eq. (8)] is positive, and the second one satisfying Eq. (9), $uH″*(φ;γ′)>$ 0, as illustrated in Fig. 11 (remember that φ is limited between 0 and π/2 for symmetry reasons).