An Alternative Formulation of the Harrison Model

: The approach advanced by Harrison puts in the spotlight the fundamental role of bistability in hysteresis modeling. The description is based on physical premises concerning irreversible thermodynamics. In the original model, the upscaling of irreversible phenomena acting on the micromagnetic level is carried out by the introduction of a phenomenological parameter β . In the present paper, an alternative approach is proposed. The outputs of individual outputs of elementary hysteresis units (hysterons) are considered like in the stop model. A veriﬁcation of the proposed model is carried out using measurement data for a praseodymium–dysprosium ribbon sample and a cylinder core made of cobalt-based amorphous material.


Introduction
Hysteresis modeling is important for practitioners dealing with descriptions of magnetic properties of materials readily used in Nondestructive Testing and Evaluation.This topic is also in the spotlight of MDPI's Applied Sciences journal, just to mention recent publications [1,2].As was pointed out in [1,3], the hysteresis phenomenon is prevalent in different areas of science, yet its occurrence in ferromagnetic materials remains a paradigm [4].
The present contribution is inspired by an excellent review work published recently [5].It extends several concepts which have been mentioned in the previous publications [6 -9].The fundamental concept of the paper is to modify the Harrison model of ferromagnetism to make it more useful for engineering purposes.The modification relies on the summation of responses of elementary bistable units operating on the quantum level using an approach like the stop models.
The original description advanced twenty years ago [10] may be perceived as a "toy model" because of its spurious simplicity.However, it is based on a profound physical interpretation; moreover, it is consistent with the laws of irreversible thermodynamics.It distinguishes the reversible and irreversible magnetization processes and offers a simple interpretation of the hysteresis loop onset as related to quantum-mechanical effects.
The irreversible magnetization processes are expressed in dimensionless units with the following relationship as: in which m stands for reduced magnetization, and h is reduced applied field strength, whereas t denotes temperature, referring to the Curie temperature.This self-consistent relationship (since reduced magnetization m appears on both sides of Equation ( 1)) puts in the central perspective the concept of the "effective" field, which is the true internal field acting on magnetic moments within the ferromagnetic substance.It is different from the externally applied field.In physical units, this concept may be written as: where the dimensionless coefficient α controlling the strength of interaction forces (referred to as the molecular field constant) is due to Weiss [11].For a friendly introduction to the Weiss concept, the Readers are referred, e.g., to the textbooks [12,13], or, e.g., to [14].The theory of effective field has been used to explain spontaneous magnetic ordering and phase transitions in ferromagnets [3], and it is widely used in solid-state physics; some widely used hysteresis models such as the Jiles-Atherton formalism [15] heavily explore and rely on the concept.It is, however, interesting to mention that the idea of mutual interactions between individual entities has influenced seemingly different research studies, e.g., those devoted to collective behavior [16,17] or focused on applied mathematics and computer science [18][19][20].In the context of research on Nondestructive Testing and Evaluation, it should be stressed that the expression for the effective field may include additional terms related, e.g., to eddy currents, thermal viscosity, mechanical stresses or demagnetization effects [21], which opens an exciting perspective to develop new types of sensors of physical quantities based on monitoring of magnetic properties.

A Brief Introduction to the Harrison Model
In the original Harrison model, the Weiss' mean field constant α is assumed to be positive, just like in, e.g., the Jiles-Atherton approach-this corresponds to the positive feedback in the system.The argument of hyperbolic tangent in Equation ( 1) is written as the weighed sum h + m.In that case, the curve given with Equation ( 1) is S shaped (if one considers m = m(h) dependence) or N shaped (if one considers h = h(m) dependence).It passes through the second and the fourth quadrants of the m-h plane; the extreme points are referred to as bifurcation points.
Equation (1) may be easily inverted analytically; the description having magnetization as its input and field strength as its output is referred to as an inverse hysteresis model.The appropriate expression becomes: where the reduced temperature t plays the role of a weighting parameter.The curve given with Equation (3) for a positive t, 0 < t < 1, has two extreme points whose coordinates are computed from equating dh/dm to zero.The fragment of the curve between these extreme points is decreasing; in ferroelectric devices, this corresponds to the "negative capacitance" effect, which recently has been the subject of intensive research [22,23] because of its potential applications, e.g., in low-power nanoscale devices.In ferromagnetic materials, in most cases, aside from in precise experiments carried out for single iron crystals [24], the intermediate region is hard to observe experimentally.
As field strength reaches the threshold value equal to h BIF , an abrupt jump to another hysteresis loop branch occurs.The magnetization value m EXIT at which both branches coincide cannot be given in an analytical form; however, it is straightforward to compute it numerically.
Two extreme cases are: • t = 0 when the curvature of the h = h(m) function disappears.Then, the h = h(m) dependence is represented by a zigzag resembling the letter "N".The value of field strength at the bifurcation point (considered equivalent to coercivity) reaches unity; • t = 1 when the hysteresis loop disappears.
Appl.Sci.2023, 13, 12009 3 of 14 In the original model, Harrison used a phenomenological parameter β, which was introduced to scale the irreversible hysteresis curves operating on the micromagnetic level up to the domain size (mesoscopic) level, on which the anhysteretic processes occurred.The afore-given statement is clear if one considers the expression (1) written not in relative units but physical ones.It reads [10]: where µ 0 = 4π × 10 −7 H/m is free space permeability, m B = 9.27 ×10 −24 Am 2 is the Bohr magneton, k = 1.38 × 10 −23 J/K is Boltzmann's constant and dimensionless α is Weiss' mean field parameter (cf.( 2)), whereas T measured in K is the absolute temperature.M s stands for saturation magnetization.From the inspection of the above-given relationship, all quantities have a clear physical interpretation, whereas β is merely a numerical means to scale the switching phenomenon, occurring at quantum scale up to the scale of magnetic domains.Harrison explained his approach in this way: " . . . the new vector quantity βm B can be called a "supermagneton".In this model the "supermagneton" is a new quantum particle, that is acted upon by the effective field" [6].
In other words, the original approach corresponds to the multiplication effect for a single hysteresis operator; thus, it means that the model is subject to averaging and homogenization, not just because of the presence of the averaged quantity α, but also because all operators are treated as the same.According to the model developer, the Weiss coefficient accounts for the exchange field between atomic moments, while the domain coefficient, i.e., β, describes the atomic moment summation process [6].
The outcome of hysteresis loops with realistic shapes were obtained by summing field strength contributions corresponding to the irreversible and reversible re-magnetization processes, cf.[10] or Figure 3 in [8].The concept of the present paper is to introduce a weighted summation from several elementary irreversible hysteresis operators and, thus, to treat the elementary Harrison unit as a generalized stop operator [5,25,26].In this context, the idea behind the present paper corresponds to the one outlined in recent publication [27], in which, however, the hysteretic characteristics of elementary units are assumed to be composed of broken line segments.The proposed approach considered in the present paper allows one to obtain better flexibility and fitting accuracy.
It should be emphasized that there exists a profound connection between the Harrison approach and theory of phase transitions advanced by Landau [28].A refinement of the Landau theory, known as the Landau-Ginsburg-Devonshire approach, plays a crucial role in the understanding of the basic physics of ferroelectricity [14,23].

An Alternative Elementary Unit ("Hysteron")
As pointed out earlier, the reduced field strength h modeled with the original Harrison description rarely reaches values close to zero; for example, for Fe-Si electrical steels, the maximal value h crit is around 0.498 [8] since the reduced temperature t ∼ = 293/1043 = 0.281 [10].Our goal was to find out a functional dependence which would behave qualitatively like the one given with (3) yet able to reach h = 1 within the bounds m ∈ −1; 1 .For such a functional relationship, it might be easier to assign weights to individual units and interpret the results.Inspired by the work [29], which stressed the flexibility of logistic mapping and its usefulness in phenomenological descriptions of ferromagnetic hysteresis, we considered the following relationship for the description of the double-well energy landscape in an elementary hysteretic unit as: where a is a constant.Following the monograph [3], we considered a connection between hysteresis phenomenon and metastability.The proper framework for examination of this relationship is thermodynamics of irreversible phenomena.If we consider a system with local memory, just three variables, H, X (Bertotti's notation, cf.chapter 2.1.4,here identified as M) and T, provide a complete characterization of states.If H and X are conjugate work variables, then G(H, T) = F (X, T) − H X is the corresponding Gibbs energy (F denotes the Helmholtz free energy of the system).To consider an energy profile with at least two minima symmetrically set with respect to the X axis, we need a polynomial of the order four.
Imposing the condition ∂F(m)/∂m = h for the equilibrium of competing thermodynamic forces [3], we can write: Thus, the dependence of the applied field versus magnetization becomes a polynomial of one order lower, i.e., of the order three.Let us notice that, due to symmetry, this polynomial has just odd order terms.Of course, it is possible to consider higher-order polynomials, which may possess local global minima, yet, in that case, the polynomials will also possess only odd terms.
Figure 1 depicts the dependencies F(m) and h(m) for an arbitrary value a = −2.
Appl.Sci.2023, 13, x FOR PEER REVIEW 4 of 14 just three variables, H, X (Bertotti's notation, cf.chapter 2.1.4,here identified as M) and T, provide a complete characterization of states.If H and X are conjugate work variables, then G(H, T) = F (X, T) − H X is the corresponding Gibbs energy (F denotes the Helmholtz free energy of the system).To consider an energy profile with at least two minima symmetrically set with respect to the X axis, we need a polynomial of the order four.
Imposing the condition ()/ = ℎ for the equilibrium of competing thermodynamic forces [3], we can write: Thus, the dependence of the applied field versus magnetization becomes a polynomial of one order lower, i.e., of the order three.Let us notice that, due to symmetry, this polynomial has just odd order terms.Of course, it is possible to consider higher-order polynomials, which may possess local global minima, yet, in that case, the polynomials will also possess only odd terms.
Figure 1 depicts the dependencies () and ℎ() for an arbitrary value  = −2.It is interesting to remark that the specific choice of (4) imposes several invariable values for some loop parameters such as remanence  = ∓1/√2 or  = ∓1/√6 ∓0.408 regardless of the  value.It can be proven that the choice  = −0.75√6−1.837 allows one to obtain ℎ = 1.A fixed remanence value might be considered a limitation since ( 4) is most appropriate for those materials for which the squareness of hysteresis loops (the ratio of remanence magnetization to saturation magnetization) [30] is close to 0.7.However, we would like to remark that alternative formulations for the dependence ℎ = ℎ() may also be applied.For example, a polynomial with odd powers is: where the subscript denotes that powers with non-zero values of coefficients also behave qualitatively in a similar way.The relationship (6) is depicted in Figure 2. The remanence values for this polynomial are  ∓0.849, and the bifurcation magnetization  ∓0.602, whereas the corresponding value of field strength (coercive field strength) ℎ 1.22.The value of magnetization when two branches of hysteresis loop coincide is approximately  0.9267.It is interesting to remark that the specific choice of (4) imposes several invariable values for some loop parameters such as remanence m rem = ∓1/ √ 2 or m crit = ∓1/ √ 6 ≈ ∓0.408 regardless of the a value.It can be proven that the choice a = −0.75√ 6 ≈ −1.837 allows one to obtain h crit = 1.A fixed remanence value might be considered a limitation since ( 4) is most appropriate for those materials for which the squareness of hysteresis loops (the ratio of remanence magnetization to saturation magnetization) [30] is close to 0.7.However, we would like to remark that alternative formulations for the dependence h = h(m) may also be applied.For example, a polynomial with odd powers is: where the subscript denotes that powers with non-zero values of coefficients also behave qualitatively in a similar way.The relationship (6) is depicted in Figure 2. The remanence values for this polynomial are m rem ≈ ∓0.849, and the bifurcation magnetization m crit ≈ ∓0.602, whereas the corresponding value of field strength (coercive field strength) h crit ≈ ±1.22.The value of magnetization when two branches of hysteresis loop coincide is approximately m out ≈ ±0.9267.A thorough examination in which combinations of coefficients for polynomials with odd powers yield qualitatively similar loops is beyond the scope of the present paper and is postponed to another publication.However, we would like to point out that the possibility to use different elementary units as building blocks for the generalized stop operator together with appropriate weighting allows one to obtain a certain level of flexibility for the description, which shall be shown in the subsequent part of this section.

The Effect from the Use of Different Elementary Units and Different Weights
During preliminary tests, we simulated the effect of just two elementary units, one given with (5), where the  value is set to  = −3, and the other given with (7).For the first unit, the loop coordinates at bifurcation points are  ∓0.408 and ℎ 1.633.Figure 3a depicts the component hysteresis loops, whereas Figure 3b shows the outcome loop for the case when no weighting was assumed, and thus the ratio of contributions was 50:50 (solid black line), and for the case when the weights were assigned following the ratio 80:20 (red dashed line).
(a) A thorough examination in which combinations of coefficients for polynomials with odd powers yield qualitatively similar loops is beyond the scope of the present paper and is postponed to another publication.However, we would like to point out that the possibility to use different elementary units as building blocks for the generalized stop operator together with appropriate weighting allows one to obtain a certain level of flexibility for the description, which shall be shown in the subsequent part of this section.

The Effect from the Use of Different Elementary Units and Different Weights
During preliminary tests, we simulated the effect of just two elementary units, one given with (5), where the a value is set to a = −3, and the other given with (7).For the first unit, the loop coordinates at bifurcation points are m crit ≈ ∓0.408 and h crit ≈ 1.633.Figure 3a depicts the component hysteresis loops, whereas Figure 3b shows the outcome loop for the case when no weighting was assumed, and thus the ratio of contributions was 50:50 (solid black line), and for the case when the weights were assigned following the ratio 80:20 (red dashed line).A thorough examination in which combinations of coefficients for polynomials with odd powers yield qualitatively similar loops is beyond the scope of the present paper and is postponed to another publication.However, we would like to point out that the possibility to use different elementary units as building blocks for the generalized stop operator together with appropriate weighting allows one to obtain a certain level of flexibility for the description, which shall be shown in the subsequent part of this section.

The Effect from the Use of Different Elementary Units and Different Weights
During preliminary tests, we simulated the effect of just two elementary units, one given with (5), where the  value is set to  = −3, and the other given with (7).For the first unit, the loop coordinates at bifurcation points are  ∓0.408 and ℎ 1.633.Figure 3a depicts the component hysteresis loops, whereas Figure 3b shows the outcome loop for the case when no weighting was assumed, and thus the ratio of contributions was 50:50 (solid black line), and for the case when the weights were assigned following the ratio 80:20 (red dashed line). (a)

The Effect from Treating the "a" Parameter as a Random Variable
By analogy to the Preisach formalism in which different "hysterons" are assigned different coercive field strengths, and thus it is natural to work with their distributions, we assume that it is possible to treat the "" parameter in (6) as a random variable and to monitor the effect of the varying number of hysterons and their variance on the shape of the outcome hysteresis loops.In this context, our analysis is like the one presented in [31].
There is a direct correspondence between the values of parameter "" and coercive field strength ℎ for hysterons whose ℎ = ℎ() dependence is given with (6), namely, Therefore, producing a distribution of the "  " parameter with the mean value  = −0.75√6 is equivalent to the introduction of hysteron distribution with coercive fields around unity. Figure 4 depicts the effect of increasing the number of hysterons.The mean value of coercive field strength is set to unity, the variance is set to 0.05.The seed of pseudo-random is in each case set to the same initial state (the appropriate MATLAB command is rng('default')); thus, the sequences of pseudo-random numbers drawn from the normal distribution have the same initial components.We consider the cases when there are ten hysterons; next, we increase their number to one hundred.It can be noticed that an increase in number of hysterons makes their distribution look more like the Gaussian distribution, which can be appreciated from a visual comparison of obtained histograms.On the other hand, the computations of the outcome characteristic become more cumbersome as the number of hysterons is excessively high.It is beyond the scope of the present paper to present an identification method which will enable one to resolve how many individual hysterons should be taken into account and how to assign weights to their outputs; however, we think that artificial neural networks (ANNs) might be useful for the purpose since the very principle of superposition of stop operators resembles to a large extent the architecture of classical feed-forward ANNs.

The Effect from Treating the "a" Parameter as a Random Variable
By analogy to the Preisach formalism in which different "hysterons" are assigned different coercive field strengths, and thus it is natural to work with their distributions, we assume that it is possible to treat the "a" parameter in (6) as a random variable and to monitor the effect of the varying number of hysterons and their variance on the shape of the outcome hysteresis loops.In this context, our analysis is like the one presented in [31].
There is a direct correspondence between the values of parameter "a" and coercive field strength h crit for hysterons whose h = h(m) dependence is given with (6), namely, 3 Therefore, producing a distribution of the "a" parameter with the mean value a mean = −0.75√ 6 is equivalent to the introduction of hysteron distribution with coercive fields around unity.
Figure 4 depicts the effect of increasing the number of hysterons.The mean value of coercive field strength is set to unity, the variance is set to 0.05.The seed of pseudorandom is in each case set to the same initial state (the appropriate MATLAB command is rng('default')); thus, the sequences of pseudo-random numbers drawn from the normal distribution have the same initial components.We consider the cases when there are ten hysterons; next, we increase their number to one hundred.It can be noticed that an increase in number of hysterons makes their distribution look more like the Gaussian distribution, which can be appreciated from a visual comparison of obtained histograms.On the other hand, the computations of the outcome characteristic become more cumbersome as the number of hysterons is excessively high.It is beyond the scope of the present paper to present an identification method which will enable one to resolve how many individual hysterons should be taken into account and how to assign weights to their outputs; however, we think that artificial neural networks (ANNs) might be useful for the purpose since the very principle of superposition of stop operators resembles to a large extent the architecture of classical feed-forward ANNs.An additional refinement to the presented approach might be based on the use of the Gaussian Mixture model for several hysteron distributions differing in mean values and variances or, e.g., some other type of kernel (Epanechnikov, triangular).In the first case mentioned, we think that the MATLAB toolbox Cupid [32] might be useful.Figure 5    An additional refinement to the presented approach might be based on the use of the Gaussian Mixture model for several hysteron distributions differing in mean values and variances or, e.g., some other type of kernel (Epanechnikov, triangular).In the first case mentioned, we think that the MATLAB toolbox Cupid [32] might be useful.An additional refinement to the presented approach might be based on the use of the Gaussian Mixture model for several hysteron distributions differing in mean values and variances or, e.g., some other type of kernel (Epanechnikov, triangular).In the first case mentioned, we think that the MATLAB toolbox Cupid [32] might be useful.Figure 5    The basic statistics on the mixed distribution are obtained by typing: [myDist.Mean, myDist.SD], which, in our example, yields ans = 1.0025 0.1012.

The Reversibility Effect
An advantage of the original Harrison model is that this description considers not just purely irreversible effects but also it explicitly speaks of reversibility, expressed with the anhysteretic curve.As already discussed in [9], there exists a discrepancy between this approach and the Jiles-Atherton (JA) formalism as far as the physical meaning of irreversible and reversible magnetization processes is concerned.Jiles and Atherton used the so-called effective field, being the sum of the applied field and the bulk magnetization multiplied by a positive constant α (by analogy to the Weiss mean field), as the argument of the function, which was supposed to describe the anhysteretic curve.Zirka et al. criticized the physical assumptions underlying the Jiles-Atherton description, indicating that the summation of irreversible and reversible magnetization components postulated by the model developers leads to the non-physical behavior of modeled curves after field reversals, which is to some extent patched using pseudo-parameter δ M , whose role is to suppress the negative dynamic susceptibility (dM/dH < 0) regions [33].However, the presence of δ M cannot be justified from the derivation of model equations.The source of problems with the physical interpretation of the JA description is the fundamental assumption that the hysteresis loop is obtained by introducing offsets (shifts) from the "anhysteretic" curve along the M (and not along the H) direction.A counter-example of the correct approach in the irreversible thermodynamics framework is the GRUCAD model [34], also studied by some of the authors of the present contribution [35].The Harrison model bears a striking resemblance to the GRUCAD model in the sense that it distinguishes irreversible and reversible contributions as field strength and not magnetization components.
If one follows the Weiss definition of the effective field, Equation ( 2), then, in the language of control science, it can be interpreted as positive feedback in the system of magnetic moments.Harrison has explicitly stated that positive feedback immediately leads to the occurrence of bistability, and the m = m(h) dependence becomes S shaped (Figures 1 and 2 resemble rather the letter "N", since we consider the mapping h = h(m)).Even Jiles and Atherton have noticed that their equation for the "anhysteretic" curve, in which the effective field is the argument of the Langevin function, may give rise to the occurrence of a hysteresis loop (what contradicts its physical meaning, cf. Figure 2 in [15]).At this point it seems crucial to remark that, in several recent publications, Silveyra and Conde Garrido [36,37] argued that the effective field should rather be considered as a demagnetization effect if it were to be considered as the argument to the anhysteretic curve; thus, the feedback would be negative.
Regardless of the academic debate on the sign of the mean field parameter, a question arises as to whether one really needs an equation for the anhysteretic curve.A reliable method, which in most cases is sufficient for practical applications, is to recover the anhysteretic curve as the middle curve from the major loop branches [8,9,38].If one accepts that hysteresis loop branches may be described with shifted hyperbolic tangent function, then an equation for the inverse anhysteretic function may be derived in a closed form [39]; thus, there is no need to introduce some complicated functions [40] or their combinations [41].
In the Harrison model, the role of the anhysteretic curve is merely to introduce a skewing of the h − m coordinate system, which means a virtual contraction of the coordinate related to field strength.This can be appreciated in the example of a point marked with a black dot in Figure 6, which originally has non-zero field value but becomes the remanence point in the irreversible field vs. magnetization coordinate system.dynamic susceptibility at coercivity (tan Θ = dℎ/d), and we are concerned with a transformation of the axis ℎ clockwise by the complementary angle /2 − Θ so that it overlaps the  axis.The figure may be compared to Figure 5 in [10] or Figure 10 in [39].
The transformation matrix which may be used for computation of coordinates of points belonging to the irreversible loop marked as (b) may be written as: and its inverse as: The forward transformation may also be interpreted as an attempt to force the straight line with non-zero slope equal to the reciprocal of differential anhysteretic susceptibility to be aligned with the  axis.
The afore-given transformation can be also used in the analysis of magnetic circuits with open air gaps, which introduce "skewing" of M-H dependencies.

Results
Exemplary modeling was carried out for a Pr8Dy1Fe60Co7Ni6B14Zr1Ti3 ribbon (a hard magnet, examined previously in [9]).The details on sample preparation are provided in [42].Figure 7 depicts the decomposition of the measured hysteresis loop into the irreversible hysteresis loop related to the bistability phenomenon and the anhysteretic curve in accordance with the approach outlined in the previous section.The figure may be compared to Figure 5 in [10] or Figure 10 in [39].
A similar affine transformation was considered recently in [39].If addition of irreversible and reversible field strength components holds, then the reverse processing of experimental hysteresis loops is also possible.Therefore, the irreversible hysteresis loops may also be recovered from the measured ones using an affine transformation where the angle marked as Θ in Figure 6 is related to a measurable quantity, i.e., the reciprocal of dynamic susceptibility at coercivity ( tan Θ = dh/dm), and we are concerned with a transformation of the axis h an clockwise by the complementary angle π/2 − Θ so that it overlaps the m axis.
The transformation matrix which may be used for computation of coordinates of points belonging to the irreversible loop marked as (b) may be written as: and its inverse as: The forward transformation may also be interpreted as an attempt to force the straight line with non-zero slope equal to the reciprocal of differential anhysteretic susceptibility to be aligned with the m axis.
The afore-given transformation can be also used in the analysis of magnetic circuits with open air gaps, which introduce "skewing" of M-H dependencies.

Results
Exemplary modeling was carried out for a Pr 8 Dy 1 Fe 60 Co 7 Ni 6 B 14 Zr 1 Ti 3 ribbon (a hard magnet, examined previously in [9]).The details on sample preparation are provided in [42].Figure 7 depicts the decomposition of the measured hysteresis loop into the irreversible hysteresis loop related to the bistability phenomenon and the anhysteretic curve in accordance with the approach outlined in the previous section.
In modeling, we considered just two hysterons without weighting to make the structure really simple; the dependence h = h(m) generating bistability for the first hysteron was h 51 = m 5 − 0.5m.The other unit was of type (5), where the value of the a parameter was taken as a = −3.3.The results of modeling are presented in Figure 8.It can be stated that, despite its simplicity, the model is able to capture the shape of the measured hysteresis curves quite well, which may be inferred from a visual comparison with previous modeling results [9] (drawn with red dotted lines).It is noticeable that the shape of the measured hysteresis curve is better reproduced with the new approach.The J-H dependencies modeled with the original model from [9] generally overestimated the values which may be considered "figures of merit", i.e., coercivity and remanence polarization, and the shapes of modeled hysteresis curves remained in qualitative agreement with the experiment.In terms of quantitative assessment, the discrepancy between coercive field strength modeled with the present model and its experimental counterpart is merely 1.05%; this value should be compared to 11.4% discrepancy for the model from [9] and the experimental value for this specific point.For remanence polarization, the respective discrepancy values are 11.4% (the present model) and 17.2% (model from [9]).In modeling, we considered just two hysterons without weighting to make the structure really simple; the dependence ℎ = ℎ() generating bistability for the first hysteron was ℎ =  − 0.5.The other unit was of type (5), where the value of the  parameter was taken as  = −3.3.The results of modeling are presented in Figure 7.It can be stated that, despite its simplicity, the model is able to capture the shape of the measured hysteresis curves quite well, which may be inferred from a visual comparison with previous modeling results [9] (drawn with red dotted lines).It is noticeable that the shape of the measured hysteresis curve is better reproduced with the new approach.The J-H dependencies modeled with the original model from [9] generally overestimated the values which may be considered "figures of merit", i.e., coercivity and remanence polarization, and the shapes of modeled hysteresis curves remained in qualitative agreement with the experiment.In terms of quantitative assessment, the discrepancy between coercive field strength modeled with the present model and its experimental counterpart is merely 1.05%; this value should be compared to 11.4% discrepancy for the model from [9] and the experimental value for this specific point.For remanence polarization, the respective discrepancy values are 11.4% (the present model) and 17.2% (model from [9]).
The present model works much better, particularly for lower values of magnetic flux density, and the value of coercive field strength is very well reproduced.The departures of modeled hysteresis curves from the experimental data points around the "knee" and in the region approaching saturation may be due to the arbitrary choice of hysterons with specified characteristics and considering just two of them.We believe that much better accuracy may be achieved if more units are taken into account.It is evident that use of many hysterons in the analysis results in better capability of reproducing experimental results (due to increased number of degrees of freedom of the description-this explains the successful story of the Preisach formalism); on the other hand, we wanted to keep the  of just two hysterons.The J-H dependency of the first one was given with expression (7); for the other one, it was of similar type, yet with different values of coefficients, i.e., For the latter hysteron,  ∓0.896, and the bifurcation magnetization  ∓0.678, whereas the corresponding value of field strength (coercive field strength) ℎ 0.894.The value of magnetization when two branches of a hysteresis loop coincide is approximately  0.966.The results of modeling are depicted in Figure 9.During modeling, no weighting of hysteron outputs was assumed.The present model works much better, particularly for lower values of magnetic flux density, and the value of coercive field strength is very well reproduced.The departures of modeled hysteresis curves from the experimental data points around the "knee" and in the region approaching saturation may be due to the arbitrary choice of hysterons with specified characteristics and considering just two of them.We believe that much better accuracy may be achieved if more units are taken into account.It is evident that use of many hysterons in the analysis results in better capability of reproducing experimental results (due to increased number of degrees of freedom of the description-this explains the successful story of the Preisach formalism); on the other hand, we wanted to keep the model structure as simple as possible in our preliminary tests.
As a second attempt to apply the developed description in practice, we carried out modeling for a cylinder core wound of Co-based amorphous alloy AMMET212.As pointed out by Harrison, his approach should be valid both for soft and hard magnetic materials.The considered material is extremely soft since the measured value of coercive field strength is 5.1 A/m (this is a value five order of magnitude lower than for the praseodymiumdysprosium sample shown in Figure 8).Again, we considered a combination of just two hysterons.The J-H dependency of the first one was given with expression (7); for the other one, it was of similar type, yet with different values of coefficients, i.e., h prim 531 = 7m 5 − 5m 3 − 0.5m (11) For the latter hysteron, m rem ≈ ∓0.896, and the bifurcation magnetization m crit ≈ ∓0.678, whereas the corresponding value of field strength (coercive field strength) h crit ≈ ±0.894.
The value of magnetization when two branches of a hysteresis loop coincide is approximately m out ≈ ±0.966.
The results of modeling are depicted in Figure 9.During modeling, no weighting of hysteron outputs was assumed.
of just two hysterons.The J-H dependency of the first one was given with expression (7); for the other one, it was of similar type, yet with different values of coefficients, i.e., ℎ prim = 7 − 5 − 0.5 For the latter hysteron,  ∓0.896, and the bifurcation magnetization  ∓0.678, whereas the corresponding value of field strength (coercive field strength) ℎ 0.894.The value of magnetization when two branches of a hysteresis loop coincide is approximately  0.966.The results of modeling are depicted in Figure 9.During modeling, no weighting of hysteron outputs was assumed.

Conclusions
Summing up to this point, in this contribution we have presented the following concepts: - The phenomenological parameter β introduced by Harrison is abandoned since it implies averaging and oversimplifies the model, which affects its usefulness;

Figure 3 .
Figure 3.The component loops (a) and corresponding outcome loops (b).

Figure 4 .
Figure 4.The histograms obtained (a) for ten and (b) for one hundred hysterons with mean value of coercive field set to unity and variance set to 0.05.
depicts the probability density function (PDF) and the cumulative density function (CDF) for a mixture of two normal distributions with different mean values, variances and weights.The plot was obtained using the commands: myDist = Mixture(0.65,Normal(0.95,0.05), 0.35, Normal(1.1,0.1)) myDist.PlotDens

Figure 5 .
Figure 5.The PDF and the CDF for a mixture of two distributions.The basic statistics on the mixed distribution are obtained by typing: [myDist.Mean, myDist.SD], which, in our example, yields ans = 1.0025 0.1012.

Figure 4 .
Figure 4.The histograms obtained (a) for ten and (b) for one hundred hysterons with mean value of coercive field set to unity and variance set to 0.05.

Figure 4 .
Figure 4.The histograms obtained (a) for ten and (b) for one hundred hysterons with mean value of coercive field set to unity and variance set to 0.05.
depicts the probability density function (PDF) and the cumulative density function (CDF) for a mixture of two normal distributions with different mean values, variances and weights.The plot was obtained using the commands: myDist = Mixture(0.65,Normal(0.95,0.05), 0.35, Normal(1.1,0.1)) myDist.PlotDens

Figure 5 .
Figure 5.The PDF and the CDF for a mixture of two distributions.The basic statistics on the mixed distribution are obtained by typing: [myDist.Mean, myDist.SD], which, in our example, yields ans = 1.0025 0.1012.

Figure 5 .
Figure 5.The PDF and the CDF for a mixture of two distributions.

14 Figure 7 .
Figure 7.The decomposition of measured hysteresis loop into irreversible and reversible components.

Figure 7 .
Figure 7.The decomposition of measured hysteresis loop into irreversible and reversible components.

Figure 8 .
Figure 8.The measured (dots) and the modeled (solid lines) hysteresis curves.Additionally, the modeling results for the same sample from Ref. [9] (Przybył et al., 2022) are shown as red dotted lines.

Figure 8 .
Figure 8.The measured (dots) and the modeled (solid lines) hysteresis curves.Additionally, the modeling results for the same sample from Ref. [9] (Przybył et al., 2022) are shown as red dotted lines.

Figure 8 .
Figure 8.The measured (dots) and the modeled (solid lines) hysteresis curves.Additionally, the modeling results for the same sample from Ref. [9] (Przybył et al., 2022) are shown as red dotted lines.

Figure 9 .
Figure 9.The measured (dots) and the modeled (solid lines) hysteresis curves for a cobalt-based amorphous alloy.

Figure 9 .
Figure 9.The measured (dots) and the modeled (solid lines) hysteresis curves for a cobalt-based amorphous alloy.