Strayfield calculation for micromagnetic simulations using true periodic boundary conditions

We present methods for calculating the strayfield in finite element and finite difference micromagnetic simulations using true periodic boundary conditions. In contrast to pseudo periodic boundary conditions, which are widely used in micromagnetic codes, the presented methods eliminate the shape anisotropy originating from the outer boundary. This is a crucial feature when studying the influence of the microstructure on the performance of composite materials, which is demonstrated by hysteresis calculations of soft magnetic structures that are operated in a closed magnetic loop configuration. The applied differential formulation is perfectly suited for the application of true periodic boundary conditions. The finite difference equations can be solved by a highly efficient Fast Fourier Transform method.

Micromagnetic simulation are often used for the characterization of magnetic materials with a certain microstructure. Since the magnetic samples are very large, only a small part of the material can be simulated. A naive truncation of the magnetic domain would lead to strong shape anisotropy originating from surface effects. Periodic boundary conditions (PBCs) allow to eliminate this influence of the surface by modeling periodic images of the primary supercell.
Most well-known micromagnetic finite difference simulation packages like OOMMF 1 , MuMax 32 , magnum. fd 3 , magnum.af 4 and Fidimag 5 calculate the demagnetization field without PBCs using an analytic expression of the demagnetization tensor of homogeneously magnetized cubes 6 combined with an efficient FFT method making use of the convolution theorem 7 . Since this method is based on an integral formulation of the magnetic strayfield equations, incorporating PBCs requires the summation over an infinite sum of periodic images. Solutions have been proposed for 1D and 2D problems 8,9 , however an extension to 3D is not possible since the occurring sums are not absolutely convergent. Using point-dipoles instead of finite magnetized cubes seems to overcome this limitation and allows true 3D periodic boundary conditions 10 . In contrast to true PBCs which require an infinite summation, some codes use pseudo PBCs where the summation is truncated after a finite number of periodic images 2,3 . Apart from the easier implementation those methods are well suited for systems of intermediate size, where a finite number of periodic images is sufficient to model the complete sample. This principle has even been applied to FEM simulation where it is called macro-geometry 11 . Since finite element calculations are based on the differential form of the strayfield equation (real) PBCs can be directly applied by providing a proper cell-connectivity which represents the periodic structure.
We propose the application of PBCs for 3D problems based on a differential form of the strayfield equations both for finite elements (FE) and finite differences (FD). We focus on the efficient strayfield calculation because it is the most time-consuming part of micromagnetic simulations. Due to the long-range interaction the strayfield is usually solved using integral formulations. Since we use a differential formulation, the discretization using FE or FD with PBCs is straightforward, however it requires the solution of a sparse system of equations. In case of FD the use of a Fourier space method allows direct inversion of the system and offers a significant speed up.

Strayfield calculation using PBCs
The magnetic strayfield h of a given magnetization m can be calculated by means of magnetostatic Maxwell's equations. Since the magnetic strayfield is curl-free a scalar potential formulation can be used: www.nature.com/scientificreports/ where u is the magnetic scalar potential and the magnetic field can be calculated as h = −∇u . Proper boundary conditions need to be defined in order to obtain a unique solution.
If the magnetization is localized in a magnetic region ⊂ R 3 the problem is called an open-boundary problem with u = O( 1 r ) as r → ∞ . Since the boundary condition are not known at the surface of the magnet, direct use of the differential formulation (1) would require the discretization of an (infinite) air domain outside of the magnet. Accurate and efficient methods for solving open-boundary problems are often based on a corresponding integral formulation. In finite-element micromagnetics, the strayfield is usually solved in the differential form (1) and a hybrid method by Fredkin and Koehler 12 is used in order to resolve the boundary conditions on the magnetic surface by means of the boundary element method. In finite-difference micromagnetics, the direct integration of the strayfield tensor 6 combined with an efficient FFT based convolution is commonly used.
When dealing with large periodic structures, only a fraction of the magnetic sample can be discretized. Assuming open-boundary conditions is no longer valid in this case. Instead, a periodic magnetization can be assumed (see Fig. 1) and thus PBCs can be used as proper boundary conditions for u.
Finite element discretization. The finite element formulation is based on the weak form of the differential equation (1) with proper test functions φ i . ∂� would be the domain boundary with a corresponding unit normal n . In case of PBCs those surface intergrals vanish, because there is no physical domain boundary. For the application of PBCs a mapping of boundary nodes to their periodic images needs to be provided in order to eliminate the corresponding degree of freedoms. We used finite element packages FEniCS 13 or firedrake 14 (via firedrake-periodicity), which offer capabilities to define PBCs.
One difficulty when dealing with PBCs in FEM is the creation of a periodic mesh. Sophisticated periodic grain structures can be created with neper (see for example Fig. 2). Compared with a finite difference discretization the finite element model provides a better geometry representation, but requires the solution of an (unstructured) sparse system of equations, which significantly increases execution time by some orders of magnitude.
Finite difference discretization. The finite difference method is based on a rectangular N x × N y × N z mesh. In order to use an efficient FFT method for the solution of the occurring system of equations, we assume that the mesh is equidistant. Thus the index set (i, j, k) is sufficient to identify each vertex of the mesh.
In the following, we use the convention that u and div (m) are defined on the grid vertices, whereas m and h are defined at the cell centers (see Fig. 3). In the case of PBCs, the number of grid points and cell centers are equal and the cell centers form a shifted grid. www.nature.com/scientificreports/ The chosen convention allows to use central differences and leads to the following discrete approximations of the continuous differential operators Using these sparse discrete operators allows to solve the discrete version of Eq. (1) Efficient implementation using a fourier space method. The discrete system (5) can be solved directly, as it has to be done in the finite element method. Due to the regular (and equidistant) grid, the FD system can also be solved in Fourier space, where all differential operators become algebraic and the system can be directly inverted. The potential u i,j,k can be represented by means of the corresponding Fourier-space potential ũ i,j,k using the Discrete Fourier Transform (DFT)  www.nature.com/scientificreports/ where is the imaginary unit and the wave-vectors k l , k m , and k n are defined by A similar Ansatz is used for the other fields m , and h . When substituting the Fourier-space representation (6) into the definition of the discrete operators (4) the spatial indices i, j, k only occur within the exponent, which results in a simple multiplicative phase factor for all neighbouring cells. Simplifying all occurring prefactors finally yields One can see that the non-local terms within the operator lead to local pre-factors within Fourier-space. Using the fact that the Fourier basis functions are linearly independent of each other, allows to explicitly express the Fourier coefficients of the strayfield h l,m,n as a function of the Fourier coefficients of the magnetization m l,m,n : Note that evaluation of ũ 0,0,0 , which represents the constant part of the potential, would lead to a division by 0. However it can be set to zero since it has no influence on the magnetic field.
Due to the use of the FFT, the resulting algorithm is very efficient. Compared with the non-periodic FFT strayfield calculation, which is based on the integral formulation, the assembly and storage of the demagnetization tensor can be avoided. Since no zero-padding is necessary the system size gets even smaller in case of PBCs. As a further optimization one can use a real-FFT, since the input magnetization as well as the resulting strayfield is purely real-valued. This leads to a further speedup by a factor of 2, as well as to a reduced storage size. The presented method is well suited for parallel execution on modern GPUs.
The FD method has been implemented within the micromagnetic code magnum.af and can be used on CPU or GPU, respectively. The FE method utilizing true PBCs has been implemented using magnum.pi. Table 1 summarizes timings of the strayfield calculation using the presented methods.

Numerical experiments
The presented method is validated by comparison with analytical calculations. The trivial case of a homogeneously magnetized bulk material leads to zero demagnetization field according to Eq. ((1)) since div m = 0 everywhere.
A non-trivial solution can be found for an infinite number of infinitely extended thin-films with thickness d 1 and a spacing between each thin-film of d 0 . Each thin-film is magnetized perpendicular to the film plane (in −x direction) with a saturation magnetization of M s . Due to symmetry considerations the resulting field only points (6)  Note that in the limit d 0 → ∞ one ends up with the well known result for a single infinite thin-film with a field −M s inside of the film and 0 outside. It is remarkable that while a single thin-film shows no external field, an infinite number of films do. The calculated fields and the corresponding potential is visualized in Fig. 4 and compared with simulation results using the presented FD method. The performance of the presented method should be further demonstrated by the calculation of the hysteresis loop of a soft-magnetic-composite (SMC) material. The material consists of isolated particles, with each particle itself consisting of several magnetic grains. The simulation is restricted to a primary cell containing only one magnetic particle, consisting of 3 × 3 × 3 grains with a size of 300 nm × 300 nm × 300 nm , as well as a nonmagnetic interparticle layer (see Fig. 5). PBCs are used to mimic interparticle interactions and to avoid surface effects. The width of the interparticle layer w gap is varied and its influence on the magnetic hysteresis is studied.
For the FE model a (non-equidistant) regular mesh is used. Each dimension is divided into N = 3N i + 1 parts, where N i is the number of divisions of each grain. The grid spacing within the grains is constant, whereas the thickness of the interparticle layer can be adjusted as desired. In contrast, the FD discretization requires using an equidistant grid, which limits the possible interparticle thicknesses to integer divisors of the grain-size. Furthermore many FFT libraries require that largest prime factor of the system size is smaller than a certain value. This is based on the fact, that the FFT performs best for system sizes N = 2 M for integer M. Performance www.nature.com/scientificreports/ decreases dramatically (more than one order of magnitude) if system sizes with much higher prime-factors are used 2 . This general limitation of the FD method also results in very large system sizes for small width of the inter-particle layer and makes the more flexible FE method still competitive. The used material parameters can be found in Table 2. The magnetic anisotropy axes within the 27 grains are randomly distributed. A homogeneous external magnetic field is applied and linearly varied from −100 mT to 100 mT with a frequency of 100 MHz . The resulting hysteresis loops for FE and FD method can be found in Fig. 6. It can be seen that a larger interparticle layer leads to a smaller hysteresis and thus reduces the hysteresis losses of the material. The dramatic influence of self-demagnetization without using true PBCs is demonstrated in Fig. 7. Due to the strong demagnetization effects, the external field range needs to be extended to −2.5 T to 2.5 T . Without PBCs the subtle effect of varying interparticle layer widths is superimposed by a much stronger finite size-effect, which additionally depends on the shape of the boundary. Since the influence of the boundary can be subtracted out only in average, extracting the desired macroscopic material properties will be much harder and less accurate.
The FE and FD simulations with true PBCs have been performed with magnum.pi, or magnum.af, respectively. The FD simulations without PBCs have been performed with magnum.af and validated with MuMax 3 . The FD simulations using pseudo PBCs have been performed with MuMax 3 .

Conclusion
The importance of using true PBCs for the calculation of material properties without the influence of surface effects has been pointed out. An efficient FFT-based FD strayfield calculation providing true 3D PBCs has been presented. This method perfectly complements methods for 1D and 2D periodic boundary conditions. Those Figure 5. Simplified geometry of the SMC material consisting of 3 × 3 × 3 magnetic grains separated with one non-magnetic interparticle layer (dark blue). The size of each grain is 300 nm × 300 nm × 300 nm , whereas the interparticle thickness is varied. Table 2. Micromagnetic material parameters used within the magnetic grains as well as for the non-magnetic interparticle layer. The easy axes of the uniaxial anisotropy are randomly distributed. www.nature.com/scientificreports/  Finite difference hysteresis curves of the SMC material for interparticle widths w gap = 23.08 nm with different boundary conditions. Without true PBCs the external field range has to be extended to ±2.5 T , in order to fully saturate the material.