- Email: [email protected]
Smoothing errors in micromagnetic calculations
ELSEVIER
Physica B 233 (1997) 348 355
Smoothing errors in micromagnetic calculations David Lewis, Edward Della Torre* Institute for Magnetics Research, The George Washington University, Washington, DC 20052, USA
Abstract
Errors resulting from approximating the atomic crystalline medium with a continuous medium are explored. Vortices which occur in real materials are not permitted in smoothed micromagnetic calculations, but occur incorrectly in numerical calculations due to discretization. Some corrections are suggested. A method is proposed to permit large angles and multiple twists in intermediate steps of the calculations. Keywords: Smoothing error; Discretization; Micromagnetic calculations
1. Introduction
Micromagnetism [1] is the study of the magnetization distribution in materials by solving for the continuous curve that interpolates the magnetization between magnetic cells. Since all materials consist of discrete atoms, there is an error, which we will call the smoothin9 error, associated with approximating the magnetization distribution by a continuous function. It has been shown that discrete and continuous systems have, in some cases, significantly different behavior [2]. To solve these equations, the continuous system is usually discretized and a set of numerical micromagnetic equations is attained. Errors associated with the discretization are referred to as discretization errors I-3]. Thus, a discrete system is normally approximated by a continuous system which is then approximated by a discrete system.
* Corresponding author.
As the governing equations of the atomic model are a large set of difference equations, and the numerical formulation requires a smaller set of difference equations, this paper suggests an alternative approach to micromagnetics for reducing smoothing errors.
1.1. Errors associated with the steady state approximation
Since the micromagnetic equations exhibit hysteresis, the solution is path-dependant. It is possible that the computed solution, although an equilibrium solution, is not the correct equilibrium solution. The only way to guarantee that the solution is the correct one, is to take the same path that the physical system does as dictated by the dynamics. It is generally accepted that using a micromagnetic version of the Landau-Lifschitz equations, instead of simply minimizing the energy, would overcome this problem; however, this also has smoothing errors which might make the calculation take an incorrect path.
0921-4526/97/$17.00 © 1997 ElsevierScienceB.V. All rights reserved PII S092 1-45 26(97)00320-7
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
1.2. Errors associated with the continuum micromagnetic model
n/2
Two-dimensional vortices, that are often present in published numerical solutions to micromagnetic problems involving thin films, cannot exist in a continuum model. If the magnetization is constrained to lie in a plane, then a vortex involves a discontinuous change in the angle of magnetization at its center. The micromagnetic exchange energy involves the gradient of the magnetization which is infinite at a vortex. In contrast, adjacent atoms, with antiparallel magnetizations, have a finite exchange energy, thus making vortexlike formations energetically possible, although unlikely for real materials. To illustrate a similar effect, the energy of a discrete wall, whose shape is shown in Fig. 1, was compared with that of a continuous wall. Fig. 2 shows how the energy of the wall calculated by the continuum model increases significantly faster than that of the discrete model and diverges as the thickness goes to zero. In the continuous case, the average exchange energy density, Wco,t, is given by
i f rTA (~30"~2 dx, \axJ
Wc°nt = 2~
349
(1)
io
0
- n/2
-
i
o
25
50
100
75
Distance (lattice spacings)
Fig. 1. A graph of the angle verses thickness of the wall. The
same shape wall was used no matter how thin the wall was.
,20
.... *-'*
7;
o ,m
Discrete model Continuous model
_
60 where A = JS2/a, the exchange constant, was taken to be 1.3 × 10-11 J/m, J is the exchange integral, a is the lattice constant, S is the spin of the magnetic rotator, and T, half the distance of integration, must be large enough to enclose the entire wall whose classical width is 6. On the other hand, in the discrete case, the average energy density, Walsc, is given by
r~
40 o
20 ""A., i
0 2
3
4
5
--'-6
L 7
• 8
" 9
3_
10
Wall thickness (lattice cells)
j=N Wdisc =
J ~, {ISl2 - Sj. Si+a}
Fig. 2. A comparison of the wall energy density of the discrete
j=O
A
i=N
- (N + 1)a2j__~0 {1 - c o s [ A 0 j ] } ,
(2)
where a is the length of the magnetic cell, N is the number of cells within the chain, and A0j is the angle difference between adjacent cells at j. In the figures, a was taken to be 3 A, independent of the wall thickness 6. This form of the energy was chosen so that the exchange energy would be zero inside a domain.
model and the continuum model. The energy of the continuum model blows up significantly faster than that of the discrete model. Also, the energy density of the discrete model has a finite limit, unlike the continuum model.
In Fig. 3, the ratio of the continuous energy density to discrete energy density, R, given by R -
W
cont
Wdisc
(3)
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
350
~7 6 5 4 3
\
2 1
0
15 6 (number of lattice cells)
0 30
Fig. 3. The ratio of the continuum energy density over the discrete energy density,. Each lattice spacing or magnetic rotator was assumed to be 4 A long.
is graphed as a function of wall thickness. We note that the smoothing error decreases to less than 1% only when the wall thickness exceeds 104 lattice cells. To obtain an error of less than 0.1% one has to have a wall thickness that exceeds 1009 lattice cells. It is seen that the energy of the discrete model is always less than the energy of the continuum model. Since domain wall thickness is proportional to the square root of the exchange energy, we expect that because of smoothing errors the continuous micromagnetic approximation computes walls to be thicker than they should be. For example, in a micromagnetic calculation one may compute a vortex. Some of these vortices may not be real. To determine whether they are real, a finer discretization is necessary. For small magnetic films, closure domains will form in order to minimize the energy. Such a magnetic structure appears to have the same structure as a vortex, except instead of enclosing a point it might have a domain wall at its center. This wall might have a large number of twists. For example, instead of being a 180 ° wall it may be a 540 ° wall. If in the domain wall calculation one keeps only principal values of magnetization orientations, when the domain wall is of the same order as the discretization size, these walls are indistinguishable to the computer.
There are several errors which will not be considered here. The first is the fluctuation in the magnitude of magnetization. Another error which is not considered is that the micromagnetic equations are based on a Heisenberg magnet, which neglect next-nearest-neighbor interactions. Strictly speaking, a Heisenberg magnet is an insulator, without an electron gas. Therefore, without the electrons to break the spin symmetry, the anisotropy energy must be treated phenomenologically [4]. The continuum micromagnetic equations replace the quantum mechanical spins with classical spins. Finally, we shall not consider the error in micromagnetic calculations associated with the assumption that the magnetization of each computation cell is uniform and the assumption that the form of interactions between cells is the same as the form of atomic interactions. However, when there is a variation in the computation cell, the cell anisotropy energy has a different form from that of the atomic anisotropy energy, since the former energy is the three-dimensional integral of the latter. The standard derivation of the numerical micromagnetic equations hides some of the underlying physics by using an intermediate stage of differential equations. The approximation that the form of the atomic interaction is the same as the micromagnetic interaction is at times a gross approximation. The continuum micromagnetic model, in spite of all its shortcomings, is a very useful model, especially for analytical work. However, at high resolution the continuum model has errors that need to be reexamined.
1.3. An alternative approach to micromagnetics We now propose an alternate method of obtaining the magnetization distribution, which overcomes some of these problems. The magnetization of the atoms between computation points is obtained by interpolation. We obtain the magnetization distribution in the standard way by minimizing the total energy, which is the sum of the exchange energy, the anisotropy energy and the Zeeman energy. The region between the mesh points of the numerical model has many atoms that form the
351
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
magnetic compasses of the Heisenberg magnet. However, it is almost always inconvenient to work with individual atoms. In the present approach, the cells that are used should be small compared with the computation cells of the numerical problem, when modeling the contribution of the atoms. Preferably, these cells should be about the size of a magnetic unit cell. Since the magnetic cell is sometimes inconvenient for calculations, a fictitious cell, referred to as the magnetic rotator, which may contain many magnetic unit cells, will be the basic building block of the approach of this paper. In this paper the formulae for a cubic lattice are derived; however, the principles of this paper can be extended to magnetic rotators of other shapes. Alternatively, noncubic magnetic cells can be approximated by cubic magnetic rotators of about the same size. Let us discretize the magnet into computation cells. The size of the computation cell should be chosen according to what is convenient for computer calculations. Each computation cell is further subdivided into magnetic rotators. The energies of all the magnetic rotators of the computation cell are summed up to obtain the total energy of the computation cell. As in the standard approach, differentiating the total energy with respect to the angle of magnetization and setting this derivative equal to zero gives a set of numerical equations that describe an equilibrium configuration. This approach makes it easier to see where the errors come from and to interpret the results of a calculation. Although we will allow large angles between adjacent computation cells during the course of a calculation, we will not address the unlikely problem of large angles between magnetic rotators. If one degree of freedom of rotation is allowed, then the exchange energy density between two adjacent magnetic rotators whose magnetization differs in orientation by an angle A0 is given by %x -
A(A0)2 a2
(4)
For two adjacent magnetic rotators in a twodimensional lattice, having the same y coordinate,
AO becomes
I(~0 : a
020
~X -~- (X --
Xij)-~X 2
0207
+ (Y- Y'J)Sg yJ'
(5)
where i andj are the indices that define the computation cell. The differentiation will be performed by finite differences. A similar expression can be obtained for two adjacent magnetic rotators having the same x coordinate. The computation cell is divided into four quadrants. As shown in Fig. 4, the forward difference is used for the derivative of 0 with respect to x, for the positive x -x~j-half of the computation cell. The backward difference is used for the derivative of 0 with respect to x, for the negative x -xu-half of the computation cell. Similarly, the forward difference is used for the derivative of 0 with respect to y, for positive y -y~j, and the backward difference is used for the derivative of 0 with respect to y, for negative y -Yu. The central difference is used for the second derivatives of 0 throughout the cell. We obtain the exchange energy density of a computation cell, Wij, by integrating the
09_ Ox
0y forward difference
backward forward ifference differenc1
\
backward difference
J
computation cell
x Fig. 4. The flow chart compares the standard method of deriving numerical micromagnetic equations and the approach of this paper.
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
352
exchange energy density over the entire cell and dividing by L 3. That is
domain wall, for simplicity. In that case, the equation for the wall in a uniaxial material reduces to
W i j = ~ [ ( O i + lj -- Oij) 2 ~- (Oij -- Oi- l j) 2
2A-~-v = K,, sin 20 - / . t o M s a t
(~20 tTX-
+ (Oij+l -- Oij) 2 "~- (Oij -- Oij_ 1) 2 + ~6(Oi+1j + O i - l j + (Oij+ 1 + Oij-1 --
x ( - H x sin 0 + H r cos 0),
20ij) 2
where H~ and Hy are the x and y components of the local magnetic field, H. The local magnetic field includes both the applied field and the demagnetizing field. The boundary conditions for this equation
20ij) 2)
1 Jr 2-~4 0i+1j+1 -- O i - l j - 1 -- Oi+xj- 1
+ 0i- 1;- 1)2],
(7)
are
(6)
where L is the length of the computation cube. An angle difference greater than 2n is interpreted as the number of twists that the magnetization of the material makes between the mesh points. Thus the correction to first order, for the micromagnetic energy is made by storing the differences in angles between adjacent cells. It is noted that the ~ and ~s terms are second-order corrections which may sometimes be neglected. The other terms agree with the standard micromagnetic small angle calculations in algebraic form only; however, in these calculations we do not take the principal value of these angle differences, since the number of twists would then be lost. We will use the terms real angles or the real values of angles for these angle differences to emphasize the fact that we are not taking the principal values. The real angles are used to calculate the exchange energy. This way, when during the course of a numerical calculation, the direction of magnetization twists by n, the direction in which it should untwist is well-defined. Using the real value of the angle to calculate the exchange energy also eliminates the unphysical, unstable equilibrium associated with two adjacent computation cells having anti-parallel directions of magnetization.
1.4. Limitations on the validity of the linear interpolation
An upper limit for the size of the computation cell can be obtained by using a one-dimensional
(8)
01x=o = 01 and 0lx=L = 02.
Temporarily ignoring the anisotropy term and the Zeeman term, a zero order approximation of the solution to Eq. (10) is O(x)=01 +
02-01
L
(9)
x.
Now, the discretized form of Eq. (10) can be written as
K~L 2
AOi +1 - AOi = ~
//0 L2 _
sin 20i
-
~'-Msat
x ( - Hx sin Oi + Hv cos 0i).
(10)
Let us define 6 = nx/--A/Ku to be the domain wall thickness [5]. Thus, if the discretization length is small compared to the wall width, then the first term on the right-hand side of the equal sign in (7) is also small. Furthermore, if the ratio of the magnetic field to anisotropy field is small, then a linear approximation is valid. For discretizations which are too coarse for the linear approximation to be valid, the standard micromagnetic equations, using the principal value of the angle to calculate the exchange energy, are also not valid. To extend the interpolation to micromagnetic rotators up to the size of a classical domain wall, a different interpolating function should be used. This function would depend on the anisotropy constant and the magnetic field. In such a situation, however, using a higher-resolution discretization would be far better. Also, the intent of this interpolation scheme is
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
to allow both for angles greater than n at intermediate solutions and for slightly larger angles between cells than would ordinarily be allowed at the conclusion of the calculation. However, solutions containing twists greater than about one radian cannot be relied upon. Such solutions should be interpreted as an indication that a higher resolution is necessary.
1.5. Three-dimensional models This scheme can be extended to three dimensions. However, if the change in spherical coordinates angles, A0 and A~b, are stored in the same manner as they were in two dimensions, then the direction of magnetization always has to retrace the same path backwards to untwist. This would be physically incorrect [6]. Although being forced to retrace the wrong path may be better than not keeping track of the twist of the material between the points at all, an even better method of keeping track of the path exists. In general, another one or two independent variables are necessary to keep track of the twist. Let 7 represent the angle of twist between two adjacent micromagnetic computation cells. If the micromagnetic computation cells are kept small enough so that linear interpolations are still valid, the micromagnetic exchange energy density between two adjacent computation cells can be written as Wex =
A(Tk,t)2 L2 ,
(11)
where 7k,t is the real angle rather than the principal value of the angle, and the indices k and l refer to two computation cells, one labeled k and the other labeled I. The real value of 7k,~ obeys cos
7k.t = COS Ok COS 0t
+ sin Ok sin 0t cos(~bk -- q~t),
Lagrange multipliers and add the terms
Ck,, = 2kj(cos Ohcos 0l + sin Ok sin 01COS(~bk-- q~t) -- COS ~k,t)
(13)
to the energy functional, where the )~k,~ are the Lagrange multipliers, and the Ck,~ are the Lagrangian constraint terms. This approach has the fringe benefit of having a somewhat simplified principal value problem to solve. Although only the principal value of these 0's and ~b's are of interest, they appear only as arguments of cosine and sine functions, which automatically handle the principal value problem correctly. Higher-order expressions for the exchange energy density can be derived, which will be similar to the two-dimensional exchange energy density. This approach works only when the correction to the demagnetizing field, corresponding to the change in the magnetic moment due to the 7's, can be neglected. When the change in the demagnetizing field due to the change in the magnetic moment cannot be neglected, two new independent variables, the difference in the polar and azimuthal angles between adjacent rotators, A0 and A~b, are introduced. Then the exchange energy density between two classical spins of cubic lattice cells or two adjacent magnetic rotators is
Wex(X,y,z;x',y',z') =
A ~-5 cos f2x,y,z~x,,~,,,z,,
(14)
where f2x,y,z;x,y,~,is the angle between the spin vectors of the two cells {x,y,z} and {x',y',z'} which are represented by magnetic rotators. The spherical angles are interpolated by simple polynomials. Each octant of the micromagnetic rotator has a different set of polynomials. Using first-order polynomials, for example, in the octant where x, y, and z are all positive, for adjacent magnetic rotators at positions {x + a/2,y,z} and {x - a/2,y,z}, we will define
(15)
(12)
where the indices refer to two adjacent computation cells. When minimizing the energy functional, we will treat ~k,t as an independent variable and use (10) as a constraint. We then use the method of
353
and (16)
D. Lewis, E. Della Torre / Physica B 233 (1997) 348 355
354
In this case, the cosine of the angle between them, 12, becomes
coso( +
os(0+)cos(0 +sin(0+ ) sin(0-
,17,
After expanding the sines and cosines, using the small-angle approximation for Ao and A~, the micromagnetic exchange energy for two adjacent magnetic rotators becomes Wex x + ~, x -
=
[(Ao)2 + (A~)2sin2 0]
= Z-~[(A0)2 + (A~o)2sin2 0]. (18) This expression is identical to the standard smallangle expression for the energy density, except that A0 and A~b are now independent variables. The constraint term Ck,z now becomes
Taylor series about 0 o. When the energy density is integrated over the computation cell, one obtains an expression, which is accurate to third order in the change in angles, that has about 12 terms. A correction to the anisotropy energy can also be derived in a similar fashion. However, the magnitude of the anisotropy energy is typically much smaller than the exchange energy and demagnetizing energy. Therefore, it is questionable if it is worth spending extra computing time in order to calculate the anisotropy energy more accurately. Another choice of simple interpolating polynomials is the linear edge fit. This would have the advantage of eliminating discontinuities at the boundaries between computation cells in the functions representing the polar angles giving the direction of magnetization of a magnetic rotator. If sine and cosine functions are used as the interpolating functions, one can arrive at a Fourier representation for the micromagnetic equations. The primary differences between this Fourier representation and other representations will be the method by which they were derived and that there will be an explicit reference to the distance between magnetic rotators in the exchange torque.
Ck,~ = 2k.t[cos Okcos 0~
2. Summary and conclusions
+ sin 0k sin 01COS(qSk-- ~bl) - (cos(0~ + A0k,3cos 0z + sin(0z + A0k,3sin 0t cos(A~bk,3)].
(19)
The resulting equilibrium equations can be simplified by making use of the fact that the 7's or, alternatively, the A0's and Aq~'s are local variables. Consequently, the final equations can be written in terms of local coordinate systems in which the azimuthal axes are along the direction of magnetisation of the mesh points. Using Eq. (22), for the exchange energy density between two computation cells, corresponds to treating sin z 0 of the energy density as if it were a constant throughout the micromagnetic rotator; however, sin 2 0 varies in the computation cell. Expanding 0, even with a linear interpolation, and then integrating over all eight octants, results in a correction, with 64 terms. The number of terms in the result can be reduced at the expense of some accuracy, by expanding sin20 in a first-order
Errors associated with the standard continuum micromagnetic model were discussed. A method of deriving alternative numerical micromagnetic equations instead of using the continuum model was described. Some corrections to the numerical micromagnetic energy were derived. The properties of the technique presented here are the following: The differences in angles of more than 21t between adjacent micromagnetic rotators have meaning. They represent the total twists in the direction of the magnetization between adjacent micromagnetic rotators. A physically more accurate model keeps track of the differences in angles between cells and never take the principal value of those angles. This eliminates any ambiguity about the direction in which the magnetization vectors should untwist. In three dimensions the twist can be made an independent variable with a constraint. However, when the solution has a large angle of twist, a higher resolution should be used.
D. Lewis, E. Della Torre / Physica B 233 (1997) 348 355
355
Acknowledgements
References
We would like to thank Jason Eicke, Gary Kahler, Ferenc Vajda, Larry Bennett, Martha Pardavi-Horvath, of the Institute for Research in Magnetics, as well as William Parke, Roger Peverley, Robert McMichael, and Michael Donahue for their many helpful discussions. This work was supported by an A T P grant from the U.S. Department of Commerce through the National Storage Industries Consortium.
[1] W.F. Brown, Micromagnetics, (Interscience Publishers, New York 1963) p. 35. [2] R. Rogers, Presented at Micromagnetism and Hysteresis Workshop, May 1996, George Washington University, Ashburn, VA. [3] Y.D. Yan and E. Della Torre, IEEE Trans. Magn. Mag24(6) (1988) 2368. [4] N.W. Ashcroftand N.D. Mermin, Solid State Physics,W.B. Saunders Company, Philadelphia, 1976. 1-5] S. Chickazumi, The Physics of Magnetism (Wiley, New York, 1964). [6] R. McMichael and M. Donahue, private communication.
Physica B 233 (1997) 348 355
Smoothing errors in micromagnetic calculations David Lewis, Edward Della Torre* Institute for Magnetics Research, The George Washington University, Washington, DC 20052, USA
Abstract
Errors resulting from approximating the atomic crystalline medium with a continuous medium are explored. Vortices which occur in real materials are not permitted in smoothed micromagnetic calculations, but occur incorrectly in numerical calculations due to discretization. Some corrections are suggested. A method is proposed to permit large angles and multiple twists in intermediate steps of the calculations. Keywords: Smoothing error; Discretization; Micromagnetic calculations
1. Introduction
Micromagnetism [1] is the study of the magnetization distribution in materials by solving for the continuous curve that interpolates the magnetization between magnetic cells. Since all materials consist of discrete atoms, there is an error, which we will call the smoothin9 error, associated with approximating the magnetization distribution by a continuous function. It has been shown that discrete and continuous systems have, in some cases, significantly different behavior [2]. To solve these equations, the continuous system is usually discretized and a set of numerical micromagnetic equations is attained. Errors associated with the discretization are referred to as discretization errors I-3]. Thus, a discrete system is normally approximated by a continuous system which is then approximated by a discrete system.
* Corresponding author.
As the governing equations of the atomic model are a large set of difference equations, and the numerical formulation requires a smaller set of difference equations, this paper suggests an alternative approach to micromagnetics for reducing smoothing errors.
1.1. Errors associated with the steady state approximation
Since the micromagnetic equations exhibit hysteresis, the solution is path-dependant. It is possible that the computed solution, although an equilibrium solution, is not the correct equilibrium solution. The only way to guarantee that the solution is the correct one, is to take the same path that the physical system does as dictated by the dynamics. It is generally accepted that using a micromagnetic version of the Landau-Lifschitz equations, instead of simply minimizing the energy, would overcome this problem; however, this also has smoothing errors which might make the calculation take an incorrect path.
0921-4526/97/$17.00 © 1997 ElsevierScienceB.V. All rights reserved PII S092 1-45 26(97)00320-7
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
1.2. Errors associated with the continuum micromagnetic model
n/2
Two-dimensional vortices, that are often present in published numerical solutions to micromagnetic problems involving thin films, cannot exist in a continuum model. If the magnetization is constrained to lie in a plane, then a vortex involves a discontinuous change in the angle of magnetization at its center. The micromagnetic exchange energy involves the gradient of the magnetization which is infinite at a vortex. In contrast, adjacent atoms, with antiparallel magnetizations, have a finite exchange energy, thus making vortexlike formations energetically possible, although unlikely for real materials. To illustrate a similar effect, the energy of a discrete wall, whose shape is shown in Fig. 1, was compared with that of a continuous wall. Fig. 2 shows how the energy of the wall calculated by the continuum model increases significantly faster than that of the discrete model and diverges as the thickness goes to zero. In the continuous case, the average exchange energy density, Wco,t, is given by
i f rTA (~30"~2 dx, \axJ
Wc°nt = 2~
349
(1)
io
0
- n/2
-
i
o
25
50
100
75
Distance (lattice spacings)
Fig. 1. A graph of the angle verses thickness of the wall. The
same shape wall was used no matter how thin the wall was.
,20
.... *-'*
7;
o ,m
Discrete model Continuous model
_
60 where A = JS2/a, the exchange constant, was taken to be 1.3 × 10-11 J/m, J is the exchange integral, a is the lattice constant, S is the spin of the magnetic rotator, and T, half the distance of integration, must be large enough to enclose the entire wall whose classical width is 6. On the other hand, in the discrete case, the average energy density, Walsc, is given by
r~
40 o
20 ""A., i
0 2
3
4
5
--'-6
L 7
• 8
" 9
3_
10
Wall thickness (lattice cells)
j=N Wdisc =
J ~, {ISl2 - Sj. Si+a}
Fig. 2. A comparison of the wall energy density of the discrete
j=O
A
i=N
- (N + 1)a2j__~0 {1 - c o s [ A 0 j ] } ,
(2)
where a is the length of the magnetic cell, N is the number of cells within the chain, and A0j is the angle difference between adjacent cells at j. In the figures, a was taken to be 3 A, independent of the wall thickness 6. This form of the energy was chosen so that the exchange energy would be zero inside a domain.
model and the continuum model. The energy of the continuum model blows up significantly faster than that of the discrete model. Also, the energy density of the discrete model has a finite limit, unlike the continuum model.
In Fig. 3, the ratio of the continuous energy density to discrete energy density, R, given by R -
W
cont
Wdisc
(3)
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
350
~7 6 5 4 3
\
2 1
0
15 6 (number of lattice cells)
0 30
Fig. 3. The ratio of the continuum energy density over the discrete energy density,. Each lattice spacing or magnetic rotator was assumed to be 4 A long.
is graphed as a function of wall thickness. We note that the smoothing error decreases to less than 1% only when the wall thickness exceeds 104 lattice cells. To obtain an error of less than 0.1% one has to have a wall thickness that exceeds 1009 lattice cells. It is seen that the energy of the discrete model is always less than the energy of the continuum model. Since domain wall thickness is proportional to the square root of the exchange energy, we expect that because of smoothing errors the continuous micromagnetic approximation computes walls to be thicker than they should be. For example, in a micromagnetic calculation one may compute a vortex. Some of these vortices may not be real. To determine whether they are real, a finer discretization is necessary. For small magnetic films, closure domains will form in order to minimize the energy. Such a magnetic structure appears to have the same structure as a vortex, except instead of enclosing a point it might have a domain wall at its center. This wall might have a large number of twists. For example, instead of being a 180 ° wall it may be a 540 ° wall. If in the domain wall calculation one keeps only principal values of magnetization orientations, when the domain wall is of the same order as the discretization size, these walls are indistinguishable to the computer.
There are several errors which will not be considered here. The first is the fluctuation in the magnitude of magnetization. Another error which is not considered is that the micromagnetic equations are based on a Heisenberg magnet, which neglect next-nearest-neighbor interactions. Strictly speaking, a Heisenberg magnet is an insulator, without an electron gas. Therefore, without the electrons to break the spin symmetry, the anisotropy energy must be treated phenomenologically [4]. The continuum micromagnetic equations replace the quantum mechanical spins with classical spins. Finally, we shall not consider the error in micromagnetic calculations associated with the assumption that the magnetization of each computation cell is uniform and the assumption that the form of interactions between cells is the same as the form of atomic interactions. However, when there is a variation in the computation cell, the cell anisotropy energy has a different form from that of the atomic anisotropy energy, since the former energy is the three-dimensional integral of the latter. The standard derivation of the numerical micromagnetic equations hides some of the underlying physics by using an intermediate stage of differential equations. The approximation that the form of the atomic interaction is the same as the micromagnetic interaction is at times a gross approximation. The continuum micromagnetic model, in spite of all its shortcomings, is a very useful model, especially for analytical work. However, at high resolution the continuum model has errors that need to be reexamined.
1.3. An alternative approach to micromagnetics We now propose an alternate method of obtaining the magnetization distribution, which overcomes some of these problems. The magnetization of the atoms between computation points is obtained by interpolation. We obtain the magnetization distribution in the standard way by minimizing the total energy, which is the sum of the exchange energy, the anisotropy energy and the Zeeman energy. The region between the mesh points of the numerical model has many atoms that form the
351
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
magnetic compasses of the Heisenberg magnet. However, it is almost always inconvenient to work with individual atoms. In the present approach, the cells that are used should be small compared with the computation cells of the numerical problem, when modeling the contribution of the atoms. Preferably, these cells should be about the size of a magnetic unit cell. Since the magnetic cell is sometimes inconvenient for calculations, a fictitious cell, referred to as the magnetic rotator, which may contain many magnetic unit cells, will be the basic building block of the approach of this paper. In this paper the formulae for a cubic lattice are derived; however, the principles of this paper can be extended to magnetic rotators of other shapes. Alternatively, noncubic magnetic cells can be approximated by cubic magnetic rotators of about the same size. Let us discretize the magnet into computation cells. The size of the computation cell should be chosen according to what is convenient for computer calculations. Each computation cell is further subdivided into magnetic rotators. The energies of all the magnetic rotators of the computation cell are summed up to obtain the total energy of the computation cell. As in the standard approach, differentiating the total energy with respect to the angle of magnetization and setting this derivative equal to zero gives a set of numerical equations that describe an equilibrium configuration. This approach makes it easier to see where the errors come from and to interpret the results of a calculation. Although we will allow large angles between adjacent computation cells during the course of a calculation, we will not address the unlikely problem of large angles between magnetic rotators. If one degree of freedom of rotation is allowed, then the exchange energy density between two adjacent magnetic rotators whose magnetization differs in orientation by an angle A0 is given by %x -
A(A0)2 a2
(4)
For two adjacent magnetic rotators in a twodimensional lattice, having the same y coordinate,
AO becomes
I(~0 : a
020
~X -~- (X --
Xij)-~X 2
0207
+ (Y- Y'J)Sg yJ'
(5)
where i andj are the indices that define the computation cell. The differentiation will be performed by finite differences. A similar expression can be obtained for two adjacent magnetic rotators having the same x coordinate. The computation cell is divided into four quadrants. As shown in Fig. 4, the forward difference is used for the derivative of 0 with respect to x, for the positive x -x~j-half of the computation cell. The backward difference is used for the derivative of 0 with respect to x, for the negative x -xu-half of the computation cell. Similarly, the forward difference is used for the derivative of 0 with respect to y, for positive y -y~j, and the backward difference is used for the derivative of 0 with respect to y, for negative y -Yu. The central difference is used for the second derivatives of 0 throughout the cell. We obtain the exchange energy density of a computation cell, Wij, by integrating the
09_ Ox
0y forward difference
backward forward ifference differenc1
\
backward difference
J
computation cell
x Fig. 4. The flow chart compares the standard method of deriving numerical micromagnetic equations and the approach of this paper.
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
352
exchange energy density over the entire cell and dividing by L 3. That is
domain wall, for simplicity. In that case, the equation for the wall in a uniaxial material reduces to
W i j = ~ [ ( O i + lj -- Oij) 2 ~- (Oij -- Oi- l j) 2
2A-~-v = K,, sin 20 - / . t o M s a t
(~20 tTX-
+ (Oij+l -- Oij) 2 "~- (Oij -- Oij_ 1) 2 + ~6(Oi+1j + O i - l j + (Oij+ 1 + Oij-1 --
x ( - H x sin 0 + H r cos 0),
20ij) 2
where H~ and Hy are the x and y components of the local magnetic field, H. The local magnetic field includes both the applied field and the demagnetizing field. The boundary conditions for this equation
20ij) 2)
1 Jr 2-~4 0i+1j+1 -- O i - l j - 1 -- Oi+xj- 1
+ 0i- 1;- 1)2],
(7)
are
(6)
where L is the length of the computation cube. An angle difference greater than 2n is interpreted as the number of twists that the magnetization of the material makes between the mesh points. Thus the correction to first order, for the micromagnetic energy is made by storing the differences in angles between adjacent cells. It is noted that the ~ and ~s terms are second-order corrections which may sometimes be neglected. The other terms agree with the standard micromagnetic small angle calculations in algebraic form only; however, in these calculations we do not take the principal value of these angle differences, since the number of twists would then be lost. We will use the terms real angles or the real values of angles for these angle differences to emphasize the fact that we are not taking the principal values. The real angles are used to calculate the exchange energy. This way, when during the course of a numerical calculation, the direction of magnetization twists by n, the direction in which it should untwist is well-defined. Using the real value of the angle to calculate the exchange energy also eliminates the unphysical, unstable equilibrium associated with two adjacent computation cells having anti-parallel directions of magnetization.
1.4. Limitations on the validity of the linear interpolation
An upper limit for the size of the computation cell can be obtained by using a one-dimensional
(8)
01x=o = 01 and 0lx=L = 02.
Temporarily ignoring the anisotropy term and the Zeeman term, a zero order approximation of the solution to Eq. (10) is O(x)=01 +
02-01
L
(9)
x.
Now, the discretized form of Eq. (10) can be written as
K~L 2
AOi +1 - AOi = ~
//0 L2 _
sin 20i
-
~'-Msat
x ( - Hx sin Oi + Hv cos 0i).
(10)
Let us define 6 = nx/--A/Ku to be the domain wall thickness [5]. Thus, if the discretization length is small compared to the wall width, then the first term on the right-hand side of the equal sign in (7) is also small. Furthermore, if the ratio of the magnetic field to anisotropy field is small, then a linear approximation is valid. For discretizations which are too coarse for the linear approximation to be valid, the standard micromagnetic equations, using the principal value of the angle to calculate the exchange energy, are also not valid. To extend the interpolation to micromagnetic rotators up to the size of a classical domain wall, a different interpolating function should be used. This function would depend on the anisotropy constant and the magnetic field. In such a situation, however, using a higher-resolution discretization would be far better. Also, the intent of this interpolation scheme is
D. Lewis, E. Della Torre / Physica B 233 (1997) 348-355
to allow both for angles greater than n at intermediate solutions and for slightly larger angles between cells than would ordinarily be allowed at the conclusion of the calculation. However, solutions containing twists greater than about one radian cannot be relied upon. Such solutions should be interpreted as an indication that a higher resolution is necessary.
1.5. Three-dimensional models This scheme can be extended to three dimensions. However, if the change in spherical coordinates angles, A0 and A~b, are stored in the same manner as they were in two dimensions, then the direction of magnetization always has to retrace the same path backwards to untwist. This would be physically incorrect [6]. Although being forced to retrace the wrong path may be better than not keeping track of the twist of the material between the points at all, an even better method of keeping track of the path exists. In general, another one or two independent variables are necessary to keep track of the twist. Let 7 represent the angle of twist between two adjacent micromagnetic computation cells. If the micromagnetic computation cells are kept small enough so that linear interpolations are still valid, the micromagnetic exchange energy density between two adjacent computation cells can be written as Wex =
A(Tk,t)2 L2 ,
(11)
where 7k,t is the real angle rather than the principal value of the angle, and the indices k and l refer to two computation cells, one labeled k and the other labeled I. The real value of 7k,~ obeys cos
7k.t = COS Ok COS 0t
+ sin Ok sin 0t cos(~bk -- q~t),
Lagrange multipliers and add the terms
Ck,, = 2kj(cos Ohcos 0l + sin Ok sin 01COS(~bk-- q~t) -- COS ~k,t)
(13)
to the energy functional, where the )~k,~ are the Lagrange multipliers, and the Ck,~ are the Lagrangian constraint terms. This approach has the fringe benefit of having a somewhat simplified principal value problem to solve. Although only the principal value of these 0's and ~b's are of interest, they appear only as arguments of cosine and sine functions, which automatically handle the principal value problem correctly. Higher-order expressions for the exchange energy density can be derived, which will be similar to the two-dimensional exchange energy density. This approach works only when the correction to the demagnetizing field, corresponding to the change in the magnetic moment due to the 7's, can be neglected. When the change in the demagnetizing field due to the change in the magnetic moment cannot be neglected, two new independent variables, the difference in the polar and azimuthal angles between adjacent rotators, A0 and A~b, are introduced. Then the exchange energy density between two classical spins of cubic lattice cells or two adjacent magnetic rotators is
Wex(X,y,z;x',y',z') =
A ~-5 cos f2x,y,z~x,,~,,,z,,
(14)
where f2x,y,z;x,y,~,is the angle between the spin vectors of the two cells {x,y,z} and {x',y',z'} which are represented by magnetic rotators. The spherical angles are interpolated by simple polynomials. Each octant of the micromagnetic rotator has a different set of polynomials. Using first-order polynomials, for example, in the octant where x, y, and z are all positive, for adjacent magnetic rotators at positions {x + a/2,y,z} and {x - a/2,y,z}, we will define
(15)
(12)
where the indices refer to two adjacent computation cells. When minimizing the energy functional, we will treat ~k,t as an independent variable and use (10) as a constraint. We then use the method of
353
and (16)
D. Lewis, E. Della Torre / Physica B 233 (1997) 348 355
354
In this case, the cosine of the angle between them, 12, becomes
coso( +
os(0+)cos(0 +sin(0+ ) sin(0-
,17,
After expanding the sines and cosines, using the small-angle approximation for Ao and A~, the micromagnetic exchange energy for two adjacent magnetic rotators becomes Wex x + ~, x -
=
[(Ao)2 + (A~)2sin2 0]
= Z-~[(A0)2 + (A~o)2sin2 0]. (18) This expression is identical to the standard smallangle expression for the energy density, except that A0 and A~b are now independent variables. The constraint term Ck,z now becomes
Taylor series about 0 o. When the energy density is integrated over the computation cell, one obtains an expression, which is accurate to third order in the change in angles, that has about 12 terms. A correction to the anisotropy energy can also be derived in a similar fashion. However, the magnitude of the anisotropy energy is typically much smaller than the exchange energy and demagnetizing energy. Therefore, it is questionable if it is worth spending extra computing time in order to calculate the anisotropy energy more accurately. Another choice of simple interpolating polynomials is the linear edge fit. This would have the advantage of eliminating discontinuities at the boundaries between computation cells in the functions representing the polar angles giving the direction of magnetization of a magnetic rotator. If sine and cosine functions are used as the interpolating functions, one can arrive at a Fourier representation for the micromagnetic equations. The primary differences between this Fourier representation and other representations will be the method by which they were derived and that there will be an explicit reference to the distance between magnetic rotators in the exchange torque.
Ck,~ = 2k.t[cos Okcos 0~
2. Summary and conclusions
+ sin 0k sin 01COS(qSk-- ~bl) - (cos(0~ + A0k,3cos 0z + sin(0z + A0k,3sin 0t cos(A~bk,3)].
(19)
The resulting equilibrium equations can be simplified by making use of the fact that the 7's or, alternatively, the A0's and Aq~'s are local variables. Consequently, the final equations can be written in terms of local coordinate systems in which the azimuthal axes are along the direction of magnetisation of the mesh points. Using Eq. (22), for the exchange energy density between two computation cells, corresponds to treating sin z 0 of the energy density as if it were a constant throughout the micromagnetic rotator; however, sin 2 0 varies in the computation cell. Expanding 0, even with a linear interpolation, and then integrating over all eight octants, results in a correction, with 64 terms. The number of terms in the result can be reduced at the expense of some accuracy, by expanding sin20 in a first-order
Errors associated with the standard continuum micromagnetic model were discussed. A method of deriving alternative numerical micromagnetic equations instead of using the continuum model was described. Some corrections to the numerical micromagnetic energy were derived. The properties of the technique presented here are the following: The differences in angles of more than 21t between adjacent micromagnetic rotators have meaning. They represent the total twists in the direction of the magnetization between adjacent micromagnetic rotators. A physically more accurate model keeps track of the differences in angles between cells and never take the principal value of those angles. This eliminates any ambiguity about the direction in which the magnetization vectors should untwist. In three dimensions the twist can be made an independent variable with a constraint. However, when the solution has a large angle of twist, a higher resolution should be used.
D. Lewis, E. Della Torre / Physica B 233 (1997) 348 355
355
Acknowledgements
References
We would like to thank Jason Eicke, Gary Kahler, Ferenc Vajda, Larry Bennett, Martha Pardavi-Horvath, of the Institute for Research in Magnetics, as well as William Parke, Roger Peverley, Robert McMichael, and Michael Donahue for their many helpful discussions. This work was supported by an A T P grant from the U.S. Department of Commerce through the National Storage Industries Consortium.
[1] W.F. Brown, Micromagnetics, (Interscience Publishers, New York 1963) p. 35. [2] R. Rogers, Presented at Micromagnetism and Hysteresis Workshop, May 1996, George Washington University, Ashburn, VA. [3] Y.D. Yan and E. Della Torre, IEEE Trans. Magn. Mag24(6) (1988) 2368. [4] N.W. Ashcroftand N.D. Mermin, Solid State Physics,W.B. Saunders Company, Philadelphia, 1976. 1-5] S. Chickazumi, The Physics of Magnetism (Wiley, New York, 1964). [6] R. McMichael and M. Donahue, private communication.