Accuracy of the analytical demagnetization tensor for various geometries

We investigate the accuracy of the analytical expressions for the magnetostatic demagnetization tensor for three different geometries by comparing these to the dipole field at far-away distances. We consider a prism, tetrahedron and cylindrical tile geometry. We show that for the prism and tetrahedron tiles the median relative error compared to the dipole field reaches below 10 − 7 at 200 tile radii away from the generating tile, while the cylindrical tile has a median relative error of about 10 − 5 at this distance. Even at a distance of 10 4 tile radii the relative error is below 10 − 2 for all tiles, also when single precision is used to calculate the magnetic field. This shows that the demagnetization tensor is an accurate way of calculating the stray magnetic field also for far-away distances.


Introduction
Calculation of magnetic field is important in a large number of fields, from MRI scanners to electric motors.Besides their common usages in motor and hard drives [17], magnetic fields are important in robotics for navigation [21,38], in magnetohydrodynamics for computing charged particle trajectories [22,3] as well as in geo-and astrophysics considering the Earth's and the Sun's magnetic fields [23,24,39].
In all applications, the calculation of the magnetic field generated by a collection of magnetic moments is needed.Unfortunately, calculating this field, which is known as the demagnetization field or the stray field, is numerically a time consuming task [1].While fast multipole methods [41,32,19] and fast k-space convolution methods [33] exists, and tensor-grid methods [15,1] are in development, a computationally efficient way of calculating the demagnetization field is using an analytical expression for the demagnetization tensor.This is the tensor relating the magnetization of the magnetic object to the produced field at a given spatial location, and it depends on the shape-dependent of the magnetic object [25].For a number of geometries with uniform magnetization, which we here term tiles, such as prisms [37], cylinder segments [28,36], tetrahedrons [29], ellipsoids internally [31] and externally [40], hollow spheres [34], cylinders of finite length [9], rings [35], Halbach cylinders [6] and a regular cuboid grid [26,16], to name a few, the demagnetization tensor can be calculated analytically.
The analytical expressions for some of the former geometries are implemented in the MagTense framework [7,8] and the Magpylib framework [30], while the regular cuboid geometry, using a Fast Fourier transform, is used in e.g. the MuMax framework [43,42] for micromagnetic calculations.
However, as shown by Ref. [10] the demagnetization tensor is inaccurate far from the generating source for a regular cuboid grid.In that work, an analytical expression for the demagnetization tensor for a regular cuboid grid, taken from Ref. [26], was studied.Using the scaling of the demagnetization field and the numerical precision of double floating numbers, the authors asserted that at tile separations greater than approximately 300 tiles, the analytical computation will contain no significant digits because of numerical cancellations.By defining the error relative to an exact (to machine precision) computation of the demagnetization tensor, they indeed found that the error reaches a value of around 1 at approximately 300 tile distance.For this reason micromagnetism simulations frameworks such as OOMMF [11] only consider the demagnetization field accurately close to the originating tile, while further away use an asymptotic expansion for the demagnetization field [10].Other micromagnetic codes, such as MaGICo [14,13,12] combine analytical and numerical computation of the demagnetization tensor at small and large distances, respectively.
While the demagnetization tensor for the regular cuboid grid has been shown to be inaccurate at large distances from the generating tile, it is not known what the accuracy is of the analytical expressions for the demagnetization tensor for all other geometries besides the regular cuboid grid investigated in Ref. [10].As these expressions are both used in micromagnetic [8,5] and magnetostatic calculations [27,30] it is crucial to investigate this in detail.

Theory
The analytical expressions for the demagnetization tensor for various geometries are based on Gauss's law for magnetism and Ampere's law.They are derived considering a single magnetized body Ω in the whole space and assuming no free currents, which means that the magnetic field vector H is irrotational everywhere.It can therefore be expressed as the gradient of a scalar field, termed the magnetic scalar potential, ϕ M , as The relation between the magnetic field H, the magnetic flux density B and the magnetization M is defined as where µ 0 is the vacuum permeability.By using this relation in the definition for the scalar potential, taking the divergence and remembering Gauss's law for magnetism, ∇ • B = 0, results in a Poisson-like equation We remark that magnetization is assumed to be nonzero only within the region Ω.Thus, the latter equation must satisfiy the following condition on the boundary of Ω, that is the surface S where we have denoted with n the outward normal on the surface S ′ and with [−∇ϕ M ] S ′ the jump of the vector field H = −∇ϕ M across S ′ .The above equations have the same form as those for the electrostatic field produced by assigned (volume and surface) charges densities and lead to identify the volume and surface equivalent magnetic charge The solution of Eqs.(3)-(4) for a homogeneously magnetized body, i.e. such that M is constant in Ω (and, consequently, ρ M = 0), is the scalar potential produced only by the surface charge density σ M = n • M and is given by [18]: The above equation has two sets of coordinates.The coordinates marked with a ′ are the coordinates of the face that creates the magnetic field, whereas the non-marked coordinates are to the point at which the field is evaluated.The magnetic field is obtained by combining Eq. 1 and Eq.5: Here it is worth noting that the gradient operator and the integration are with respect to different sets of coordinates: the gradient operator is with respect to the unprimed coordinates, while the integration is with respect to the primed coordinates, as can also be seen in the equation.Equation ( 6) can be formally written as the product of a tensor, called the demagnetization tensor, N, and the (constant) magnetization as The demagnetization tensor contains all the geometrical information on the generating tile, and can be calculated analytically for several geometrical objects or tiles, such as a prism [37], cylinder [28,36], tetrahedron [29], ellipsoid [31,40] and hollow sphere [34].For each of these the magnetic field at a point r can then be found if the geometrical specifications of the tile, i.e. its size, orientation and location in space, is known.In terms of this work, the interest is on the accuracy of the demagnetization tensor for the prism [37], tetrahedron [29] and cylindrical tiles [28,36], as shown in Fig. 1.
We can assess the error in the demagnetization tensor by considering the fact that far away from any magnetic object, the magnetic field H should approach that produced by a magnetic dipole.This can be seen by considering that the magnetic scalar potential can be expressed by means of the multipole expansion [20] in powers of 1/r as where the ϕ n M term is proportional to 1/r n+1 .In magnetostatics, one has ϕ 0 M = 0 (zero net magnetic charge) and, at distances far away from the source generating the field, it is usual practice to restrict oneself to the dipole term ϕ 1 M , as this will dominate the expansion above [2].
The dipole field is given as where the magnetic moment m is given by the product of the magnetization and the tile volume and r is the unit vector from the tile center to the point at which the field is evaluated.We can thus define the relative error in the magnetic field calculated using the analytical demagnetization tensor in a given point as i.e. as the norm difference in the field between the analytical demagnetization field and dipole field, normalized by the norm of the dipole field.

Geometry
In the following, we take the magnetization be M = [1, 1, 1] in arbitrary units.
The following three tile geometries were considered: prism, cylindrical tile and tetrahedron.These are shown in Fig. 1.Their specifics are discussed individually below.
For each tile considered we define a tile radius and tile center, which is used to quantify the distance from the tile to the point that the magnetic field is evaluated in.Then we generated 121 spheres with radii varying uniformly in logspace between 0 and 10,000 tile radius and with each sphere having 2452 points uniformly distributed on its surface.The field is then evaluated in each of these points and compared with the dipole field using Eq.(10).

Prism tile
The prism tile is defined by its center coordinates, [x 0 , y 0 , z 0 ] and its size, [a, b, c].Its center is simply the center of the prism tile and the tile radius is defined as r tile = (a/2) 2 + (b/2) 2 + (c/2) 2 .The analytical expression for the demagnetization tensor is taken from Ref. [37] in the implementation in the MagTense framework [7].

Cylindrical tile
The cylindrical tile is defined by its center position, (r 0 , θ 0 , z 0 ) and its extent in the radial, angular and height directions, (dr, dθ, dz).There is no analytical expression for the centroid of the cylindrical tile, so the centroid position and the tile radius is found using a Monte Carlo approach.
The analytical expression for the demagnetization tensor exists in different forms.The full solution in closed analytical form was recently derived in Ref. [36] and is implemented with all special cases in the Magpylib framework [30].However, a semi-analytical solution where only four components of the tensor, N xy , N xz , N yz and N zy , have been expressed fully analytically while the five remaining components contain integrals that have to be evaluated numerically, have also been published in Ref. [28].This solution in implemented in the MagTense framework [7].In the following the fully analytical solution from Magpylib [30] is used for the cylindrical tile unless otherwise stated.

J o u r n a l P r e -p r o o f
Journal Pre-proof 4. Results

Tile geometries
To show the accuracy of the demagnetization tensor, we initially calculated the relative error as defined above for the previously mentioned 2452 points on a spherical surface for 121 spheres with radius varying uniformly in logspace between 0 and 10,000 tile radius.For each of the 121 spheres, we computed the median relative error.This was done for the three different tile geometries for two different tile sizes, one with equal and one with unequal dimensions.
We consider a prism tile which has dimensions [a, b, c] = [1, 1, 1], and another prism tile with dimensions [a, b, c] = [1, 2, 3], i.e. tiles with different aspect ratios.For the cylindrical tile, we consider a tile with center position [r 0 , θ 0 , z 0 ] = [1, π/8, 0] and size [dr, dθ, dz] = [1, π/8, 1] as well as a tile with center position [r 0 , θ 0 , z 0 ] = [2, π/8, 0] and dimensions [dr, dθ, dz] = [2, π/8, 3].Finally, for the tetrahedron, we consider a tile with vertices The relative error as function of distance from the tile in terms of the previously defined tile radius is shown in Fig. 2a.As can be seen from the figure, the relative error in the field is initially high, as the field close to the tiles are not identical to the dipole field, simply due to the shape of the tiles themselves.As the radius at which the field is evaluated increases, the relative error decreases substantially until around ∼ 10 2 tile radii, at which point the relative error starts to increase for all tile types.Even at 10 4 tile radii the median relative error is below 10 −2 for all tiles.
The reason that the relative error increase after ∼ 10 2 tile radii is mostly due to numerical cancellations of two similar surfaces on the tiles.In fact, remembering that the demagnetization tensor is constructed by integration of Eq. 6 across the surface of a given tile, it is so that opposite surfaces will contribute with terms with opposite sign, due to the term n(r ′ ) • M(r ′ ), i.e. the dot product of the surface normal and the magnetization of the tile.For very large distances, two such opposing surface terms approach each other in magnitude, as the distance between the surfaces becomes much smaller than the distance to the field evaluation point, while having opposite sign.This can cause numerical cancellations to occur when the terms are summed, similarly to what is reported in Ref. [10].
Interestingly enough, the analytical demagnetization tensor expressions for the prism and tetrahedron tiles seems to be much more numerically robust compared to the cylindrical geometry.However, as it can also be seen from the figure, this is only the case for the tiles with an aspect ratio of one, which indicates that the numerical instabilities in the analytical expressions in the demagnetization tensor cancel out for these tiles.Note also that the cylinder tile reaches a constant relative error.This will be discussed subsequently.
Instead of considering the field at equal distances from the generating tile, we can also consider level surfaces (isosurfaces) of the field magnitude.The   norm of the dipole field is given by where θ is the spherical polar angle.Isolating for r and converting from spherical to cartesian coordinates, the coordinates for constant values of ||H dipole ||, normalized by the norm of the magnetization, can be determined.The median relative error evaluated at these points, for decreasing values of ||H dipole || corresponding to increasing distances from the tile, is shown in Fig. 2b.This figure displays the same trend as Fig. 2a, but now the value of ||H dipole ||/||M|| at which the dipole field and the demagnetization field diverge can be determined, which is around 10 −7 to 10 −10 .

Tile aspect ratios
It is of interest to investigate the influence of the tile dimensions on the relative error.Shown in Fig. 3 is the relative error at 10 3 tile radii for the different tiles as function of specific tile dimensions.
For the prism tile the a-and b-dimension are varied with respect to the c-dimension, such that the tile dimensions are [a/c, b/c, c].For the cylindrical tile we consider a tile with center position [r 0 , θ 0 , z 0 ] = [k r * 1, π/8, 0] and size [dr, dθ, dz] = [k r * 1, π/8, k z * 1].Finally, for the tetrahedron, the coordinates of the vertices v 3 and v 4 are varied independently, such that the tile is given by As can be seen from the figure, varying the size of the prism tile results in a symmetric distribution of the relative error with respect to the a and b dimensions of the tile.Generally, the more asymmetric the tile is, that larger the relative error becomes.For the cylindrical tile, the relative error becomes large for the tile with large radial extent and low thickness, and generally the tile should have about the same radial extent as its thickness in order for the relative error be to low.Finally, for the tetrahedron, again the more asymmetric the tile becomes, the larger the relative error is.

Singularities for the cylinder tile
As noted earlier, the cylindrical tile warrants special attention.Shown in Fig. 4 is the relative error in all field evaluation points for the cylindrical tile for both the Magpylib [30,36] and the MagTense implementations [28,7] for the cylindrical tile with center position [r 0 , θ 0 , z 0 ] = [1, π/8, 0] and size [dr, dθ, dz] = [1, π/8, 1].As can be seen the values are distributed fairly closely around the median value for r tile < 300 tile radius.However farther away from the tile than this very large relative error values appear, with the relative error being up to 10 5 .These values originate from the single evaluation point laying on the same z-value as the tile center, z 0 , i.e the evaluation point at the sphere northand south pole centered on the tile.For this point exactly the numerical error becomes very large.This happens both for the semi-analytical (MagTense), and the fully analytical implementations (Magpylib).

Numerical precision
It is also of interest to compare the numerical precision of the analytical demagnetization tensor when using single precision computations.Shown in Fig. 5 is the relative error for the tiles introduced earlier computed in both single and double numerical precision, i.e. a prism tile with dimensions [a, b, c] = [1, 1, 1], a cylindrical tile with center position [r 0 , θ 0 , z 0 ] = [1, π/8, 0] and size [dr, dθ, dz] = [1, π/8, 1] and a tetrahedron tile with vertices v As can be seen from the figure, reducing the numerical precision to single means that the lower relative error plateaus around 5 * 10 −8 for both the prism and tetrahedron tiles, while the cylindrical tiles are not affected by the lower numerical precision.

Data statement
The demagnetization field computed in all cases are directly available from Ref. [4].This includes the script used to compute the field from MagTense and Magpylib directly.

Conclusions
We have investigated the accuracy of the analytical expressions for the demagnetization tensor by calculating the stray magnetic field and comparing this to the field generated by a dipole, as these should be identical at distances far from the tile.The relative error, i.e. difference between the magnetic field of the tile calculated using the analytical expression and the dipole field, reaches below 10 −7 at 200 tile radii away from the generating tile for the prism and tetrahedron tiles, while the cylindrical tile has a median relative error of about 10 −5 at this distance.Even at 10 4 tile radii the median relative error is below 10 −2 for all tiles.
If the tiles have an aspect ratio very different from 1, the relative error can increase, although such geometries are typically not used for generating meshes.Finally, we show that using single precision for calculating the magnetic field is as accurate as using double precision except for prism and tetrahedron tiles at distances of around 10 2 tile radii, where a small gain in accuracy is possible using double precision.

Figure 1 :
Figure 1: The three tile geometries considered in this work.

Figure 2 :
Figure 2: a) The median relative error in the magnetic field as function of distance from the tile center for three different tiles with two different tile sizes as given in the legend and b) the median relative error in the magnetic field as function on surfaces of constant ||H dipole ||.

Figure 3 :
Figure 3: The log10 of the relative error as function of two of the tile dimensions for the prism, tetrahedron and cylindrical tiles.

Figure 4 :
Figure 4: The relative error as function of distance in units of the tile radius for both the Magpylib and MagTense tile implementations.The data points for each radius is uniformly distributed on a sphere.

Figure 5 :
Figure 5: The relative error as function of distance from the tile, for different tile types, in double and single numerical precision.