Ferromagnetic hysteresis by Heisenberg partition function and its processing methods in nanoparticle magnetization modeling

The Heisenberg ab initio theory of magnetization is developed to apply for multilayer 1 nanoparticles. The theory is based on distribution and partition functions modification with 2 account the difference between exchange integral and closest neighbour numbers, that change 3 the system of resulting transcendental equation for magnetization and its reversal to form either 4 a paramagnetic type curve or hysteresis loops patterns. The equations are obtained within the 5 Heisenberg partition function construction by Heitler diagonalization of energy matrix via irreducible 6 representations of permutation symmetry group. A combination with the Gauss distribution gives 7 the explicit expression for the partition function in the asymptotic limit at large spin range in terms 8 of transcendent function. The exchange integral, as a parameter of the equation of state (material 9 equation) is evaluated from Curie temperature value by means of a formula derived within the 10 presented theory. Methods of data processing from the simultaneous solution of the material equation 11 system are proposed. The multi-valued function of hysteresis loop is found by combination of 12 graphical approach and special procedure for elimination of mistaken peaks and prolapses of the 13 patterns. The theory and computation methods are applied to spherical particles with separate 14 surface layers consideration. The contribution of the surface layers, that are specified by number of 15 closest neighbors and exchange integrals into overall magnetization, is studied for two-layer and 16 three-layer models, that are discussed and compared graphically. 17


Mathematical and physical background. Energy partitioning through the distribution function
whereĝ is a Gibbs canonical distribution operator, used for description of equilibrium states of physical systems. In (1), we noteĤ as the Hamiltonian operator, k -Boltzmann constant, T -absolute temperature. The Gibbs operator depends on external parameter a of the system under consideration.
We left the only such parameter in the magnetism phenomena description, its choice is done below. The normalization condition for (1) reads: The thermodynamic internal energy of a system U in the statistical approach is represented as the mean value of Hamiltonian: where F = U − TS (S -entropy), F is a thermodynamic potential, Helmholtz free energy. Being the function of state, its differential is exact. Other thermodynamic variables are defined similarly to (3), as the mean values of corresponding operators. Proceeding the previous point, the connection between the free energy and the partition function is settled as: F = kTln 1 Z . Plugging it in (2) The identity (4) is differentiated, taking the tiny change of temperature dT and external parameter δa: The work δA is defined by where b is the internal parameter. As the result, the first law of thermodynamics is obtained: Then, as the Pfaff's form of free energy F is exact, it could be written as follows: The following expression is substituted in the derived first law of thermodynamics (7), accounting the rule for linear combinations of independent differentials da and dT, by which the pairs of variables can be further correspondingly equated as: On this ground, the state equation is derived: The pair of parameters, necessary for current magnetization patterns visualization, are M = b, H = a. M -is the magnetization, H -is the external magnetic field force, As the result, the connection between prominent quantities -the magnetization and the external 120 magnetic field force -is set. Quadratic surface lattice 6 Simple cubic lattice 8 Cubic space-centered lattice 12 Cubic face-centered lattice # Correspondent condensed matter structure for listed z values substitution.

Heisenberg partition function and exchange interactions properties
122 W. Heisenberg builds foundations of the ferromagnetic theory, basing on the general quantum 123 statistics fundamentals, described in the preceded section with the effective use of universal symmetry 124 of a multi-electron system with respect to the permutation group and the Pauli principle, which makes 125 it generally applicable for various materials [11,29]. He incorporates the  results for 126 exchange interaction description and, as the main tool for energy spectrum evaluation, the Heitler's 127 one [15], that express the multi-electron Hamiltonian in terms of exchange integrals and characters of 128 irreducible representations of the permutation group. 129 The exchange interaction is described by the number of closest neighbors z -it is settled at the table 130 1. Note also, that the exchange integral J -decreases exponentially with the increasing distance between 131 atoms. The number of closest neighbors for Iron is z Fe = 8, for Nickel and Cobalt z Ni = z Co = 12, 132 basing on the knowledge about particular materials crystal structure; it is visualized on Figures 23 133 (Iron) and 24 (Nickel and Cobalt), placed in Appendix. 134 Heisenberg derives an expression for energy with the help of the characters of irreducible representations, listed in [15] and corrected in [11]. Its distribution within the σ−subspace is given by Gauss function. The partition function Z by such construction is: m -magnetic quantum number, which denotes spin projections, ∆E 2 -mean quadratic energy deviation, f σ -dimension of the irreducible representation, s -maximum value of m, n -number of electrons, k -the Boltzmann constant, T -Kelvin temperature. Two assignations of dimensionless variables are made: and here β is a combination of internal material parameters, α -the dimensionless external magnetic field (H), µ B -the Bohr magneton. After the summation of (12), we arrive at The most probable magnetic quantum number m 0 is determined by that is a form of the equation of state m 0 (α) in terms of the dimensionless pair of quantities. Next, note, that the symbol K is a function, which logarithm is small comparatively with the second term at large n and m 0 /n. It does not depend on m 0 , so it is omitted. As in computations we mostly worked with Z, one more notation is introduced for the argument of hyperbolic function: between paramagnetic and ferromagnetic ranges.

139
The terms β m 0 n and β 2 m 0 nz define the smoothness of hyperbolic tangent functions with argument ξ.

140
In turn, the smoothness of hyperbolic tangent functions determines whether m 0 is single-valued at 141 each solution or if exists a set of m 0 values - Figure 1.

142
To define the parameter β values at borders of switching from a paramagnetic state into a ferromagnetic one, we work with the condition: As it is seen from equations (16,17), the tangent of the curve (17) for m 0 n = 0 goes under a smaller angle with the x-axis than the tangent of the curve (16), that define number of intersections of the curves, hence -the magnetization kind. Below we give a lot of examples of the case. The condition (18) is equated in the critical temperature point T = θ and β here gets the assignation β θ . It brings us to the quadratic equation: The previously mentioned expression with square root 1 − 8 z figures further in computations by the Heisenberg theory (inclusively for the exchange integral J evaluation) as the consequence of a quadratic equation solution way via the discriminant. The solution consists of a pair of roots: General expression β = zJ kT in critical point takes the form: and, with the roots (20) application for the exchange integral values we get: We, as W. Heisenberg, take one of two possible roots as physical one. For defining the most convenient values of β θ , we plug the numbers of closest neighbors z values in (20), which are: If z = 8, then If z = 10, then β If z = 12, then β According to In a primary source the Taylor expansion is done for the expression (16), taking small values of ξ. It brings us to the approximate description of a high temperature region by Heisenberg theory (namely, higher than Curie point temperatures): The cubic and higher terms are discarded due to its minor significance, i.e.: Further the m 0 is equated with the expression for x. Referring to the previous point of exchange integral J reevaluation, β θ is expressed here through Then, by the binomial factoring rule, substitution of β θ from (28), putting 1 T out of the brackets, tear minus out of the brackets to exchange the positions of T and θ in θ − T, the resultant expression of approximation is For Iron (z = 8) it contracts to On the other hand, the ξ expression (17) is to be prepared for high temperature magnetization without any approximations. Similarly, we substitute β θ from (28): and include the last expression in the system of two equations (the visualization of computation results is given in Figure 26 and Figure 27) of Appendix: The system is solved with the temperature change, accounting that α parameter should be small; for 156 example, α = 0.01. The mean magnetization value is taken per one electron, so n = 1, as we compute    The expression (17) for mean magnetization per electron acquire the next method application: substituting the partition function with n = 1 (Z 1 ):

Mean per electron magnetization
The most probable magnetic quantum number m 0 per one electron is determined by (16) is as it is a derivative of a compound function. Next, we compose the system of two equations: where m 0 is the main variable. The system (35) is considered to be the most convenient for software to find out the solutions in case of paramagnetic curves creation, but for the ferromagnetic case we recommend to imply the second way of gaining the intersection points with the usage of the system with inverted hyperbolic function: The solutions of (35) and (36) are the same, but still, the processing and its speed might be different 171 due to different solutions density and location relatively to the axes - Figure 3. Before building the behavior -in order to find out the consequence of functions sets behavior and properties.

174
According to Figure 3 we conclude, that with applying the first system (35) the functions set 175 thickens in one point, while with applying the second system (36), the functions set thickens along 176 the straight line, parallel to y 1,2 axis. It makes the difference in data processing speed and quality. We        For nanoscale magnetism studies we regarded working with spherical nanoparticles, as their shape is the most convenient for primary estimations due to its uniform central symmetry - Figure 8. In defining the number of active magnetic electrons we adapted a common formula for a spherical  As it is shown, the grids represent the order of layer division. The yellow grid shows the inner sphere, which is accounted homogeneous for its parameters values. The space inside the thick orange grid is segregated, when we define the second layer. The space inside the rare orange grid is additionally segregated, when the third layer is defined. Without any layer segregation, phenomena from space outside yellow grid are omitted.
shell volume. It is the representation of an annulus in a three-dimensional case [9], defined as the subtraction of the smaller nested sphere with radius r 2 from the bigger total sphere with radius r 1 , for its volume we write: We considered two main configurations of nanoparticle models useful in nanoscale physics: a two-layer 243 spherical nanoparticle - Figure 9, and a three-layer spherical nanoparticle - Figure 10. Nonetheless, 244 more number of layers can be distinguished using the same principles, illustrated throughout this 245 section. In preparing the number of electrons parameter n, we come to calculation of each atomic layer 246 consistency.

247
Each additional layer has the thickness, equal to the thickness of a mono-atomic layer. The 248 approximate diameter of each element atom is taken as an elementary layer unit by d = 2r, where r -249 atom radius. Each element value is represented in Table 3.   We have obtained needed visual information about particles structure and their size parameters from Figure 9 and Figure 10. R is a total nanoparticle radius. Moreover, we can establish, that in implying new layers we just exhale them without additional transformations of previously modeled ones. There are three number of electrons values, separate for each layer: n B -number of electrons in a bulk part, located in the center of a particle, n I -number of electrons in an intermediate layer with the thickness of mono-atomic layer, n S -number of electrons in a surface layer with the thickness of mono-atomic layer. They are derived transforming (37) for our case: Then, the most suitable nanoparticle size is established numerically - 1 · 10 −9 3 · 10 −9 5 · 10 −9 1 · 10 −8 1.5 · 10 −8 5 · 10 −8 1 · 10 −7 n S 93 1045 3030 1.25 · 10 4 2.84 · 10 4 3.2 · 10 5 1.28 · 10 6 n B 44 2678 14209 1.25 · 10 5 4.37 · 10 5 1.69 · 10 7 1.36 · 10 8 Figure 11. Establishing the smallest border of nanoparticle multilayer model workability. Bluemonoatomic layer atoms number, red -central part atoms number. X axis -atoms number, Y axisnanoparticle radii. Left and right plots differ in the atoms number range.
size for our modeling is with 5 nm ≥ R ≥ 3 nm. With decreasing the radius up to 1 nanometer, the 253 surface layer becomes thicker than the bulk layer -in this case film models are to be applied instead 254 of bulk. As for the bigger radius -it lesser causes the impact, that is why up to 100 nanometers it 255 gradually becomes inconspicuous - Figure 11.

256
The next, microscopic, parameters for nanoparticle division: exchange integral J, number of 257 closest neighbors z, and the combined material parameter, composed from J, z and thermal energy kT.

258
The last constitutes general material characteristics, as it was mentioned in the Curie magnetization 259 section. The evaluation of these group of parameters is done sequentially on the ground of physical 260 reasons.

261
Analyzing Figure 12, we observe the patterns of exchange interactions parameters. Where the lines are thicker, the interaction is stronger. As the interaction patterns of a bulk layer (marked red) fully cover the whole nanoparticle, the exchange integral remains the same as in previous model of mean magnetization per electron. So we just use the expression, which was primary derived in Curie region section: It is evaluated through the concrete values; z B is the same as z from mean per electron model. However, both monoatomic layer exchange patterns (marked blue) and corresponding exchange integrals require changes. We supposed that the values of exchange integral increases with approximate percentage ratio. For intermediate layer its value is bigger in 0.15, in surface layer -in 0.25.
Analogically it is made for number of closest neighbors for each part, implying the same command of discarding non-integer numbers parts as in Figure 34 of Appendix. The difference is in other percentage relation: The percentage in this case is selected to obtain homogeneous decrease of each with accounting the crystal structure of each material, see Figure 23 and Figure 24 at Appendix. The parameter β is just separated for each layerβ B , β I , β S , and z B , z I , z S evaluation results are inserted in its evaluation: The results of z, J and β particular evaluation are represented in Figure 35 for Cobalt mono-elemental 262 nanoparticle (Appendix).  Let the interaction between layers is neglected, except the contributions in neighbours numbers and exchange integrals value, supposing that the total Hamiltonian of the whole body is just the sum of Hamiltonians of each layer. Under such assumption the partition function can be separated for both models (Z I I and Z I I I -index I I for a two-layer nanoparticle, index I I I for a three-layer nanoparticle), using the property of exponent e a+b = e a · e b : whereĤ II is the Hamiltonian of a two-layer nanoparticle. Hence and, after implying both of substitutions into Gibbs operator general definition (1), we obtain It leads to the Gibbs operator expression for two layersĝ wI I (during this section index "w" is used for denoting energy levels instead of canonical "n", in order to avoid mixing up with denotation of electrons number "n"):ĝ defining energy spectra for each layer, when Hamilton operator is already applied Here E wB -energy spectrum for a bulk layer, E wS -energy spectrum for a surface layer.

268
For defining three layers operators we do alike transformations, transforming the general definition (1) into 1 As the result, the Gibbs distribution expression for three layersĝ wI I I iŝ The connection of the energy spectra with partition function for each layer is given by E wB -energy spectrum for a bulk layer, E wI -energy spectrum for an intermediate layer, E wS -energy spectrum for a surface layer. Let us proceed to the connection of obtained transformations. All the values are included in ξ B , ξ I and ξ S ,that looks identically, for example, Hence, the partition function for the corresponding layer reads On this ground, for magnetization variable m 0 we have: As the result, the expression for m [BS] 0 is derived analogically as in Section 4, for mean magnetization per electron model. The difference is in number of electrons impact. For three-layer spherical nanoparticle the expression (16) changes the form with substituting the partition function Z 3 : 0 . In building the resultant magnetization patterns, the system of equations from mean magnetization per electron model (35) is changed into the pairs of systems like this: The alternative variant, like (36), which is easier for data processing, writes, e.g. as: where m The resulting sums m  The slope angle of each layer auxiliary function is different, but they show that the saturation 295 region of every layer matches for the whole nanoparticle.
296 Figure 14 gives the explanation, where ferromagnetic pattern origin is placed at high temperatures   Then, the comparison of two-layer and three-layer nanoparticle models is done. As it is described 301 in Figure 15, the intermediate mono-atomic layer gives a significant impact in total nanoparticle 302 magnetization.
303 Figure 16 proves, that paramagnetic behavior remains the same in temperature slightly upper 304 than the Curie point. On this base we consider this range applicable in practical usage.

305
If to visualize nanoparticle size studies (Figure 17),a multi-layer model gives a possibility to define 306 mono-atomic layers impact for nanoparticles with various sizes from Table 4, accounting the border 307 of film magnetization emergence from Figure 11. An Iron element is determined applicable in such 308 situation modeling with z B = 8.

309
The represented graph shows the border value of a nanoparticle bulk-significant size, as the bulk       Multi-layer hysteresis of a spherical nanoparticle is to be compared as well. For this aim Figure   330 21 is algorithmically built, Figure 22 is made with manual checking to emphasize each separate 331 loop pattern changes. Every layer gives not only the impact in a bulk hysteresis, but the total loop 332 deformation, which can be interpreted as the one of rectangular loop deformation methods.

392
The overall algorithmic sequence allows applications in general field of hysteresis loops 393 construction and analysis. Namely, the issue of hysteresis occurrence is solved numerically, and, 394 not less important, the prolapse and peaks generated by discrepancies are controlled and eliminated.

395
The example we trace is the hysteresis loops problem, appears in many areas of physical knowledge,    to intersection points lists creation, we produce the next loop operations, which are similar to the 413 programming loops, created in Figure 28.    Float(round(Float(x))) commands are put to use in discarding non-integer parts of a number, 439 as every n denotes a concrete set of elements nodes, which has to be integer in any cases. eval f (x) 440 command is primary intersected in the previous command, forming a compound request.

441
The loop occurrence checking repeats Figure 32 procedure, which founds out multiple   Figure 30 (B). The first pair shows the single solution for each curve. The second part shows the non-standard case for Iron (analogically with Figure 4). The third pair shows the homogeneous hysteresis, which gives rectangular loops. All the resulting patterns are illustrated here -Figure 33

On substruction 448
The first subtraction stands in defining the module of subtractions of two closest points - Figure   449 39. This operation adverts to the magnetization curve pattern and works with its correction. As the  The loop is to be separated into five areas - Figure 43 (B). The first area contains two first  In central part rebuilding for a single-curve pattern, the smoothness procedure is applied as well.        (or vice versa). So, we create two corresponding lists of peak indexes and discard these two correct 484 peaks from these lists. For maximums: the first peak is always correct, so we discard the first peak   (red), grey peaks show the minimums behavior. The first black peak and the last grey peak are correct, that is why they are not compared with the black point value (peak coefficient). Central peaks are discarded. The first sample is mean per electron model loop, the second -bulk nanoparticle loop, correction results of both are presented in Figure 42. First is for mean per electron loop correction, second is for bulk nanoparticle loop correction. Blue circles are loop division coefficients. The prolapse, which module is smaller than these coefficients, is discarded. In (C) different percentage multipliers are implied in a loop division coefficient evaluation, for illustrating the flexibility of subtraction checking procedure. The last of (C) is the summation for a non-mistaken graph; it shows the prospective of hysteresis deformation; so the summation checking can be used not only for corrections, but as a deformation procedure of primary correct curves.
When we defined borders of loop division (four border indexes in total), we obtained, at which 497 indexes maximum and minimum patterns will be subdivided (into three parts each pattern and it is duplicated with minus sign for negative loop values -it is done for saving symmetry too.
As it was mentioned before, prolapses and peaks originate randomly, that is why in some cases 505 the summations require recheck. For example, when there is only one mistaken peak for the whole 506 loop pattern -all the borders match at one value then. The lacking borders can be defined manually, 507 but there exists one more computational solution. All the obtained barders are repeatedly processed.

508
If all the borders are already defined, algorithm displays the message on a screen, reporting: which 509 concrete prolapse part of the whole magnetization pattern is properly processed. If the software lacks 510 at least one index -it is immediately created at this step. Each border corresponds to one of two the computation scale is very large. It is hard to process data with huge magnitudes (10 29 order, for 548 example), especially if the hardware is not primary purposed for this aim.

549
It is easy to prove, that the observed peaks and prolapses are built mistakenly. The first evidence 550 is -exists the full observable central symmetry of all possible solutions. When the mistaken omission 551 happens -the symmetry is contravened. On this ground all the intersections points recording is always done in exclusively one sequence -from α1 to α2. The second evidence is -visual comparison of an 553 automatic list of i-points and their manual analogue. Before the development of automatizing, all the 554 points were written out by hand, so the first loops were built manually. The last is the reason, why the 555 following modeling is recommended to be used in its simplest configurations with small magnitudes.

556
It was found out in the result, the matching curves with various scales are identical, so the difference is 557 always in their axes scale, not in a shape of graph increments.

558
It is important to mention, professional and narrowly targeted for transcendental tasks packages 559 should solve the system without highlighted prolapses. If such difficulty still arises, the detailed 560 succession how to correct the graphs is outlined here and some recommendations are included as well. The following correction stage is considered to be the most complicated, as it implies complex 564 smoothing of the sharpest peaks. Though, all prolapses fixing is recommended to start from here. prolapsing patterns coefficients generalization.

571
Before data converting we are to create the peak axis. It numerates each peak by its index number.

572
As it is a draft step of algorithm, we can leave the α axis by the way of subtraction/summation 573 argument. The following alternative just makes visualization more precise - Figure 38.

574
When correct peaks indexes are rewritten in a new list, our task is to define the borders of curve 575 cutting. Each function pattern is cut into several areas. Namely, into four areas, but the fourth area the last peaks from this list and nominates the next border role for two future borders: index before 581 the first peak index and index after the last peak index. Though, both borders have to be centrally 582 symmetrical for an ordinary single-pointed curve, so there are some additional operations for border 583 reflections, described in the next paragraphs of this section. If the curve draft is primary correct, the 584 subtraction checking step can be omitted or both the borders could be set equal to 0 at this step. It 585 leads us to "cutting" the correct curve in its center, so the central correction region will contain 0 points.

586
The second subtraction and the third procedure (summation) are developed for corrections of 587 hysteresis loops. They work analogically to the curve pattern correction. The main difference is in 588 separate checking of loop maximums (marked red) and loop minimums (marked blue). The corrections 589 can be done for one loop part only, and then it can be reflected, producing the same resulting loop -590 Figure 44.

591
However, the aim of this step of work was to develop this checking with its future prospective in 592 hysteresis biases modeling as well. Or, even, to forecast such subtraction/summation applications in 593 other functions building, without any limitations due to their symmetry.