Estimating the demagnetization factors for regular permanent magnet pieces

In this paper, methods for finding the demagnetization factor, specifically for permanent magnets, are considered and compared. Widely applied and cited literature expressions are compared to high resolution numerical modeling in order to establish the applicability of the expressions. Overall, the expressions are found to be very closely correlated with the modeling results. In the second part, a simple geo-metrically based method to find the demagnetization factor is proposed and compared to the numerical modeling results. Fairly good correspondence is found, indicating the applicability of this simple method


I. INTRODUCTION
Demagnetization correction of magnetometry data is important but is often overlooked or ignored. When measuring magnetic properties using open flux loop systems, such as a vibrating sample magnetometer (VSM), it is important to correct the recorded data for magnetostatic demagnetization in order to remove the influence of sample shape. For closed flux loop systems, such as hysteresigraph, this is not the case due to the closed flux path. However, for these systems, the sample must typically be quite large and regularly shaped. For VSM measurements, the sample is typically smaller and can be less regularly shaped.
It is well known that when magnetizing a sample, the dipole field from it will generate an opposing magnetic field, known as the demagnetization field, H d . The magnitude of this field increases toward saturation as it is proportional to the magnetization of the sample, M. For ellipsoids and other idealized shapes, such as infinite plates or rods, in an applied magnetic field, H appl , the internal field in the sample, H, can be written as where the demagnetization factor, N, is a constant value. In other cases, N and thus H will be position dependent within the sample, making the situation more complex. In principle, the demagnetization factor is a tensor as it relates two vectors, but in general, it is a good approximation for soft magnets to assume that M is aligned with H appl and thus just look at the magnitudes of the fields, with N as a position dependent factor. For hard magnets in the absence of strong applied fields, M is aligned with H, albeit in the opposite direction. When correcting magnetometry data for the effect of demagnetization, it is useful to reduce N to a single value, allowing the conversion of an M vs H appl dataset into an M vs H dataset. In order to do this, an average value of N that will be representative of the sample should be determined. In general, such an average will represent a certain spread of N within the sample, but the bulk of the sample will be close to the average value. This is especially true if an orientation is chosen where the average N is small. 1 For a number of simple geometrical shapes, such average values have been estimated by expressions given in the literature. These include prisms, 2-6 ellipsoids, 7-9 and cylinders, 10-14 as well as a number of more complex shapes. [15][16][17][18] Although the demagnetization factor is constant within ellipsoids, the exact value is not simple to derive. The methods, expressions, and approximations mentioned in the literature vary in complexity and thus have direct applicability to correct experimental data. In addition, it is not always clear what the quality of the expressions and approximations is and thus how accurate a correction of experimental data would be. Notably, Chen and co-workers developed elaborate but very accurate methods, as shown in Refs. 5, 11, and 13. Here, some of the widely used expressions for correcting magnetization data will be evaluated. For prisms, the expression by Aharoni 3 is very widely used, as indicated by the 1214 citations to the paper (July, 2021). For ellipsoids, the classical expression from 1945 by Osborn 7 is still in wide use. The paper has 2780 citations. For cylinders, many expressions have been suggested. Here, the one by Joseph 10 has been chosen due to its relative simplicity. The paper has 188 citations.
The expressions for ellipsoids and cylinders require the use of elliptical integrals. Previously, this was a complication and has led to papers being published with tabulation and approximations for the expressions, e.g., Ref. 8. However, today, the original expressions can relatively simply be implemented in software packages, such as MATLAB.
Here, we will compare the demagnetization constants derived from the equations mentioned in the literature with the ones found by detailed numerical modeling. Specifically, this will be carried out for parameters relevant to hard permanent magnets. The aim is to evaluate how accurate the expressions are and determine whether they apply for permanent magnets.
In addition, a simple method for estimating the values of N for any of the three shapes is given. This allows the estimation of the demagnetization field while carrying out measurements in the laboratory and thus the possibility of adjusting the measurements to cover the desired range of internal fields. In addition, it allows determination of which orientation to choose for the sample during measurement if this is not a priori obvious.

II. MODEL
In order to analyze the effect of demagnetization, a numerical model was created in the software package COMSOL Multiphysics. For the initial work of modeling a prism, a prism was created with side lengths a, b, and c, as shown in Fig. 1. Around this, a volume with three times the side lengths was defined for high resolution meshing. Finally, a volume with ten times the side lengths was defined around the prism, with a lower resolution mesh.
A permanent magnet is characterized by a strong magnetocrystalline anisotropy, seeking alignment of the magnetic moments. Such a magnet can be modeled using the state function B = Br + μ r μ 0 H, where Br is the remanence of the magnet, μ r is the relative (or recoil) permeability, and μ 0 is the vacuum permeability. With the applied field along the x direction, the relative permeability is anisotropically defined as a tensor, Here, the value of μ ∥ = 1.05 is used along the direction of magnetization, in accordance with the stated value from magnet manufacturers. It represents minor imperfections and misalignments, which are removed when applying the magnetic field, increasing the magnetization toward saturation. The value of μ = 1.17 in the perpendicular directions is the accepted value found in Ref. 19. This off-axis term represents the strong magneto-crystalline anisotropy. Using the method described in Ref. 20 and assuming a saturation of M 2 0 = 1.6 T, the value of μ = 1.17 is equivalent to an anisotropy of which is close to the conventional value for type NdFeB magnets 21 of about 5 MJ/m 3 . In the remaining volumes in the model, the state function is that of vacuum B = μ 0 H. The average demagnetization constant of the prism is found by setting a remanence value, typically 1.4 T, and then analyzing the resulting magnetization and internal field in the direction of the applied field. When the applied field is zero, Eq. (1) becomes N = − H M . For calculating N in a certain direction, the volume averages of the field and magnetization components in that direction are used. Thus, Nx = − H x M x . The stability of the calculated value of demagnetization was tested as a function of the level of mesh refinement. It was found that above ∼10 6 mesh elements, the solution value of N becomes stable within less than 0.1%.
The side lengths a, b, and c were varied between 0.1 and 1.0 in steps of 0.1, giving a total number of model runs of 1000. In addition to the prism, cylindrical rods and ellipsoids were also modeled. The ellipsoids were defined with three semi-major axes j, k, and l along the cardinal axes, as shown in Fig. 1. These were varied as in the prism case, also giving 1000 model runs. For the cylindrical rods, the height of the cylinder, d, was varied in the range 0.02-20. Keeping the radius, r, fixed at 0.3, this gave 86 model runs. Here, all model runs were conducted with the remanence along both the cylinder axis and perpendicular directions.
In the same way as for the prisms, regions with a high mesh density were defined around the cylinders and ellipsoids at three times the diameters or lengths, and regions with a lower density were defined at ten times the diameters or lengths.

III. RESULTS
To probe the dependence of the remanence and of the relative permeability, two initial tests were conducted. In both, a prism was modeled with b = c = 1 and a = [0.3 0.6 1 2 3], with the field along the a direction. In the first run, the permeability is kept as mentioned above, and the remanence is varied in the range 0.2-2.0 T. In the second run, the remanence is kept constant, and the relative permeability is set to be isotropic and varied in the range 1-100.
The results from the variation of the remanence show, as expected, that within machine precision, the demagnetization is exactly the same, independent of the remanence.   Figure 2 shows the demagnetization factor vs the isotropic relative permeability. It is observed that the demagnetization factor decreases as the permeability increases. The reason for this is less alignment with the remanence. It is an indication that although the demagnetization factor is often considered to be only a function of geometry, in reality, it is also a function of the permeability, albeit fairly weak. The expected result of N = 1 3 for the cube is found at the lowest value of the permeability. The dashed lines show the corresponding values calculated from the Aharoni expression in Ref. 3. It should be noted that the values relevant for permanent magnets are very close to 1, as stated above.

A. Prisms
With the initial results in mind, the full 1000 model runs for the prisms were conducted. For each one, the average demagnetization was calculated as described above.
Using the geometrical input, the corresponding demagnetization constants according to the method presented by Aharoni 3 were also calculated. The initial observation of Fig. 3(a) is that the expression by Aharoni is a good approximation to the modeled demagnetization factors as the red dots closely follow the dashed line, indicating a unity slope. In Fig. 3(b), the model values have been subtracted from the Aharoni values. We see that the two sets of values agree to within around 0.01, with the deviation being largest in the most extreme cases of very large or very small demagnetization. In these cases, the geometries are close to those of infinite rods or sheets, where the approximations of N = 0 or N = 1 are often employed. From Fig. 2, we know that the correlation would have been much worse if we had modeled using higher values of the relative permeability. Thus, it has been shown that for permanent magnets or materials with properties similar to these, the expression given by Aharoni is a very good approximation to find the demagnetization factor.

B. Ellipsoids
Figure 4(a) shows the correlation between the classical expression by Osborn 7 and the numerical model. Due to the way the expression is defined by Osborn, where j ≥ k ≥ l is required for the ellipsoid, the data plotted are limited to this subset. The three colors indicate the demagnetization of these ellipsoids along the three cardinal axes. Looking at the difference in Fig. 4(b), it is observed that this correlation is indeed remarkably good for permanent magnets as the maximum deviation is below 0.01.

C. Cylinders
For the cylinders, again, the initial observation of Fig. 5(a) is that the expressions are good. Here, both the axial and perpendicular demagnetization factors have been modeled, calculated, and plotted against each other. Looking at the difference in Fig. 5(b), we see that in both directions, the deviation is less than 0.025, and again mainly in the most extreme cases. This shows that the expression derived by Joseph 10 is good for permanent magnets.

IV. SIMPLE METHOD FOR ESTIMATING THE DEMAGNETIZATION
As a way of estimating the demagnetization factor, a straightforward method based on projected areas is proposed. This method can readily be used in the lab while preparing and initiating experiments as it is simple to apply.
Put simply, the projected area of the face of the sample orthogonal to the applied magnetic field is divided by the sum of the projected areas in all three Cartesian directions. Thus, for a cube of side length a with the field normal to one of the sides, the projected area normal to the field is a 2 , while the total projected area in all three directions is 3 ⋅ a 2 , so N = 1/3. Similarly, for a sphere of radius r, the projected area is πr 2 , and in all three directions, it is 3 ⋅ πr 2 , again giving N = 1/3, as expected.
In the following, it is demonstrated that this expression gives a fairly good estimate of N for any variation in prisms, cylinders, or ellipsoids.

A. Prisms
For each one of the modeled prisms discussed in Sec. II, the demagnetization factor has been calculated by the projection method according to with the field along the direction of a. The letters are permuted for other field directions. The classical textbook case of infinite sheets or rods magnetized normal to the plane can be approximated by letting a → 0 or a → ∞, resulting in N proj → 1 or N proj → 0, as expected.
The values found by this method have been compared to those found in the numerical modeling. In Fig. 6, it observed that the simple expression actually gives a fairly good approximation to the value from numerical modeling. The difference between the two shows a maximum deviation of just above 0.02, but again primarily at low and high values of N.

B. Ellipsoids
For an ellipsoid characterized by three semi-major axes j, k, and l, the projection method gives with the field along the l direction. The letters are permuted for other field directions. Figure 7 shows that also for the ellipse situation, the projection method is a fairly good approximation. From the difference plot, the maximum deviation is seen to be below 0.04 for any of the directions of magnetization.

C. Cylinders
For a cylinder with radius r, length d, and the field along the axis of the cylinder, the projection method gives Considering the textbook case of an infinite cylinder magnetized along the axis, which can be approximated by letting d → ∞, the result is N proj → 0, as expected. With the field perpendicular to the axis of the cylinder, the projection method gives Here, the result for the textbook infinite cylinder magnetized perpendicular to the axis becomes N proj → 1/2, as expected from symmetry. Again, by comparing the values from the simple expression to those from the numerical modeling (Fig. 8), we see that the correlations are good but not quite as good as the one for prisms. This can clearly be seen from the difference plotted in Fig. 8(b). The difference reaches values just above 0.03 in the worst cases. It is smallest when closest to the extreme cases, opposite to what was observed for prisms and ellipsoids.

V. DISCUSSION
Throughout Secs. III and IV, the absolute deviations have been considered. Conventionally, deviations would be given in relative terms. In Fig. 9(a), the results from Fig. 6(b) normalized by the demagnetization factor from the model have been replotted. This shows that the relative deviations become quite large especially for the lowest values of N. However, when correcting magnetization data using Eq. (1), it is actually the absolute deviation that is important as this multiplied by M is the error in the calculation of the internal field. An example of this is shown in  Here, an experimental dataset from a 3 × 3 × 3 mm 3 permanent cube has been corrected using the value of N = 0.33. Deviations of ±0.03 from this value have also been plotted. These deviations represent the worst of the absolute deviations from the projection method, but the effect of these is still quite small. For N = 0.33, the deviation is equivalent to about 10%. Had the value of N been, e.g., 0.15, the deviations would represent 20%, but the actual change in the corrected data would have been the same. Thus, it is the absolute deviation that is important for practical data correction.
A concern about the applicability of results presented here could be the validity of the numerical model compared to experimental measurements. A thorough comparison will not be conducted here but could be considered in a future publication. Instead, it is noted that numerical modeling of magnetic structures using COMSOL or other software packages has proven to be extremely accurate, compared to analytically calculated expressions 16,17 or experimental measurements of constructed assemblies. [22][23][24] This is AIP Advances 11, 075028 (2021); doi: 10.1063/5.0060897 11, 075028-6 ARTICLE scitation.org/journal/adv a testament to the ability of numerical models to calculate the local magnetic fields in all parts of the assemblies, based on the shapes and properties of the magnets. Numerical modeling is considered a reliable tool to use for the design of even very accurate magnet assemblies. [25][26][27][28] In Secs. III and IV, the presented comparisons have throughout been related to the average values of the demagnetization factors. In Fig. 10, the local demagnetization factor has been spatially resolved within versions of the three geometries. It is noted that, as expected, for the ellipsoid, the demagnetization factor is uniform. This is due to the uniform nature of the magnetization and internal magnetic field in this geometry. For the prism and cylinder, a significant distribution is observed. In all three cases, the parameters have been chosen to give relatively high average demagnetization factors of 0.538, 0.585, and 0.542, for the prism, ellipsoid, and cylinder, respectively. These high values increase the spread of N within the prism and cylinder. 1 However, even though the distribution is quite broad, the deviation is most significant at the edges, and an average value can be usefully applied. The applicability of the average value is validated by the correspondence to the single values calculated by FIG. 10. Spatial variation in the demagnetization factor within the three considered geometries. In all, the field is applied along the x-direction. expressions given in the literature and also found by the proposed simple expression.

VI. CONCLUSION
The aim of this study was to evaluate the quality of the expressions widely used in the literature for estimating the demagnetization factors of different shapes. This was carried out by comparing the values found by the expressions to detailed numerical modeling. Specifically, for permanent magnet materials with relative permeabilities close to one, it has been shown that expressions mentioned in the literature for prisms, ellipsoids, and cylinders are very good. For prisms and ellipsoids, the deviation between expression and model values is less than 0.01, often much less. For cylinders, it is a little more, up to 0.025. This shows that using expressions mentioned in the literature for correcting data measured on permanent magnets or other low permeability materials is valid.
In addition, a simple and readily applicable way of finding the demagnetization factors has been presented. It is based on calculating the ratio of projected surface areas of the sample, which can be carried out in your head. The values found by this expression are less accurate that those found by expressions mentioned in the literature but still fairly accurate. The deviation between the model and projection based values is always less than 0.04 for all shapes but often much less than this. Thus, this method is applicable for quickly calculating approximate demagnetization factors.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.