2.2.3.2 Magnetostatic potential. The DADI method

The micromagnetism poses particular emphasis in the exact calculation of the demagnetizing field. In the continuous micromagnetic approach it is described by the Maxwell equation in the absence of currents $\nabla \times \mathbf{H}=0$ and, therefore, there exists a potential called scalar magnetostatic potential that satisfies:

$$\mathbf{H_{dem}}=-\nabla U_{mag}.$$ (2.26)

This potential satisfies the Laplace equation outside the material and the Poisson one inside the material:

$$\begin{array}{cc} \nabla^2 U_{mag}=0 &\mbox{Outside ferromag... ...thbf{M(r)} & \mbox{Inside ferromagnetic material} \end{array}$$ (2.27)

with the boundary conditions:

$$\frac{\partial U_{mag,in}}{\partial \mathbf{n}}-\frac{\partial U_{mag,out}}{\partial \mathbf{n}}=4 \pi \mathbf{M}\cdot \mathbf{n}$$ (2.28)

where $\mathbf{n}$ is the unit normal to surface of the ferromagnetic sample.

The sources of the magnetostatic field are called magnetostatic charges. The surface charge is defined as $\sigma =\mathbf{M}\cdot \mathbf{n}$ and the volume charge as $\rho=-\nabla \cdot \mathbf{M}$. The magnetostatic energy integral extends to all the space, not only the ferromagnetic body, and the minimization of the magnetostatic energy will seek the reduction of the field outside the ferromagnetic body. This can be achieved reducing of the magnetic charges and that fact is stated as the pole avoidance principle. This field gives origin to domains and magnetization inhomogeneities. The generic solution of a field from its sources is found in traditional electromagnetism

$$U_{mag}=\int_V \frac{\rho (\mathbf{r'})} {\vert\mathbf{r}-\mathb... ...'+\oint_S \frac{\sigma (\mathbf{r'})} {\vert\mathbf{r}-\mathbf{r'}\vert} d^3r'$$ (2.29)

and from Eq. (2.26)

$$\mathbf{H_{dem}}=-\int_V \frac{(\mathbf{r}-\mathbf{r'})\rho (\ma... ...-\mathbf{r'})\sigma (\mathbf{r'})} {\vert\mathbf{r}-\mathbf{r'}\vert^3} d^3r'.$$ (2.30)

In micromagnetic calculations most of the computation time is spent in the calculation of the magnetostatic field. The dipolar interaction is non-local and long ranged. The field in dipole-dipole representation is a sum of $3N$ elements, where N is the number of atomic moments and the numerical integration of Eq. (2.30) scales like $N^2$ and this represents too much computational time even for today advanced computers. Due to this reason, several advanced methods have been used to solve the magnetostatic problem. The methods used are Fast Fourier Transform(FFT) [Yuan 92], Finite Element Method (FEM) [Chen 97], hybrid FEM - Boundary Element Method (FEM-BEM) [Fredkin 90], Fast Multipole Method (FMM) [Brown 04] and other advanced methods. The efficiency of each method depends on the implementation and on the simulated physical problem. We use a different method called the Dynamic Alternating Direction Implicit (DADI) [Gibbons 98]. The underlying principle of this method is to add a fictitious time derivative to Eq. (2.27):

$$-\frac{\partial \psi}{\partial t}=L\psi+\rho$$ (2.31)

with $L=\nabla^2$. Next, we integrate the resulting equation forward in time to the time-asymptotic state. At this point the fictitious time derivative goes to zero and we have the desired solution $\psi=U_{mag}$. The ideal case would be the performance of only one step of integration of Eq. (2.31). The time step $1/\omega$ should be chosen to achieve the convergence in as few steps as possible. In a hysteresis loop the solution is obtained from the solution of the previous field step, thus the test potential is very close to the solution and the method converges in few steps. Near coercivity large changes of magnetization and scalar potential are expected. In this situation, the complete convergence is needed.

The operator $L$ is replaced by its finite difference equivalent. For example the $L_x$ contribution is

$$L_x=\frac{\psi_{i+1,j,k}+\psi_{i-1,j,k}-2\psi_{i,j,k}}{\Delta^2}$$ (2.32)

with $\Delta$ - the size of the spatial discretization and $\psi_{i,j,k}$ the value of the scalar potential in the element on the node $(i,j,k)$. The method consists of three passes

$$\begin{array}{cc} (-\omega+L_x/2)\psi^{n+1/3}=(-\omega+L_x/2... ...+1}=-\omega\psi^{n+2/3}+(L_z/2)\psi^n & zpass. \ \end{array}$$ (2.33)

The magnetostatic volume charges are included replacing the operator divergence by its finite difference equivalent. In the surface cells we have to change this term by the surface charges. For example the contribution of a cell which occupy the position (i,j,k) will be:

$$\rho=\frac{1}{2}(M_{i+1,j,k}-M_{i-1,j,k}+M_{i,j+1,k}-M_{i,j-1,k}+M_{i,j,k+1}-M_{i,j,k-1}).$$ (2.34)

The solution of these equations involve tridiagonal matrices, which are linear systems easy to solve. This can be easily seen from Eq. (2.32) because the left hand side only includes the values of $\psi_{i+1,j,k}$, $\psi_{i-1,j,k}$ and $\psi_{i,j,k}$. The evaluation of the potential is the successive solution of Eqs. (2.33). From the current value $\psi^n$, we obtain $\psi^{n+1/3}$, and this yields the value $\psi^{n+2/3}$ and, finally, the value of the potential in the next step $\psi^{n+1}$.

This method needs the addition of free (zero magnetization) space in the simulated region to allow the potential to slowly tend to its infinity value. This value is arbitrary because the equation is differential and we choose it, for convenience, to be zero. This is also the case in FEM and FFT methods(zero-padding). The hybrid FEM-BEM method presents the advantage of not requiring discretization elements outside the simulated sample. In the DADI method periodic boundary is easily implemented if we consider $\psi(1)=\psi(N+1)$. The resulting tridiagonal system is replaced then by a periodic tridiagonal system, which is also very simple to solve.