Numerical simulations for magnetostatics

For the free-electron maser design we describe in Prototype free-electron maser design, it is essential to model the magnetic fields in the undulator as well as in the collimating magnets and magnetic lenses. This nearly conducting computational magnetostatics, since the magnetic fields do not have an exact analytical solution for all but the simplest cases. We will give a brief guide for how these simulations are performed.

Introduction to magnetostatics

We will start by giving a basic overview of magnetostatics. Magnetostatics describes the behavior of charges in time-independent magnetic fields. This is only an approximate description of electromagnetism, but works well when the fields vary very gradually (or not at all) in time. The two fundamental equations of magnetostatics are two of the Maxwell equations (Gauss’s law and Ampere’s law) in the time-independent regime, and take the form:

$$\nabla \cdot \mathbf{B}=0,\nabla \times \mathbf{H}=\mathbf{J}$$

Where $\mathbf{J}$ is the total current density, $\mathbf{B}$ is the magnetic field, $\mathbf{H}$ is the demagnetizing field, and the latter two are related by the constitutive relation of magnetostatics:

$$\mathbf{B}={\mu }_0(\mathbf{M}+\mathbf{H})$$

Where $\mathbf{M}$ is the magnetization of the material. Solving these equations for different boundary conditions allows us to find the fields, and therefore the trajectory of particles within these fields. The study of boundary-value problems (BVPs) in magnetostatics, and in particular the specific BVPs that we encounter in modelling free-electron lasers/masers, are discussed within this guide.

Magnetic fields for the undulator

Purely exterior case

To start, we are concerned with finding the magnetic field $\mathbf{B}$ outside the magnet. It is also possible to find the fields in the interior of the magnet, but this is less useful for our particular applications, so we first consider the external fields only.

Outside a permanent magnet, assuming there are no significant currents present and magnet is surrounded by a non-permeable medium ($\mu \approx {\mu }_0$), then $\nabla \times \mathbf{H}=\mathbf{J}$ simplifies to $\nabla \times \mathbf{H}=0$. Then it is possible to write $\mathbf{H}$ in terms of a potential ${\psi }_m$, where ${\psi }_m$ is the magnetic scalar potential, and satisfies $\mathbf{H}=-\nabla {\psi }_m$. By substituting into Gauss’s law ($\nabla \cdot \mathbf{B}=0$) where $\mathbf{B}={\mu }_0\mathbf{H}$ in the surrounding region, we arrive at:

$$\nabla \cdot \mathbf{B}=\nabla \cdot {\mu }_0(-\nabla {\psi }_m)=-{\mu }_0{\nabla }^2{\psi }_m=0$$

Therefore, after straightforward simplifications, we arrive at Laplace’s equation for the magnetic scalar potential:

$${\nabla }^2{\psi }_m=0$$

This can be solved in the exterior case to arrive at the following analytical solution¹:

$${\psi }_m(\mathbf{r})=\frac{1}{4\pi }{\oint }_S\frac{\mathbf{M}\cdot \hat{\mathbf{n}}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{A}^{\prime }$$

For a uniform permanent magnet, the magnitude of the magnetic field when no external field is applied (that is, $\mathbf{H}=0$) is defined to be equal to its magnetic remanence ${\mathbf{B}}_r$; thus from the constitutive relation $\mathbf{B}={\mu }_0(\mathbf{M}+\mathbf{H})$ for magnetic fields, we have²:

$${\mathbf{B}}_r={\mu }_0\mathbf{M}\leftrightarrow \mathbf{M}=\frac{{\mathbf{B}}_r}{{\mu }_0}\text{(inside magnet)}$$

Therefore, using this definition, the integral can be evaluated explicitly, though an analytical expression can only be found for simple geometries; the integral must be evaluated numerically (for instance, via Monte-Carlo methods or numerical quadrature) for all other cases. Here, we will take the example of a rectangular bar magnet, for which the surface integral can be done analytically. This corresponds with the magnetic geometry and type used in our undulator design, which uses permanent magnets only (as shown below).

Thus, in our case, $\mathbf{M}=(0,M_y,0)$ where $M_y=B_r/{\mu }_0$, and the integral becomes:

$${\psi }_m(\mathbf{r})=\frac{B_r}{4\pi {\mu }_0}{\oint }_S\frac{d{A}^{\prime }}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}=\frac{B_r}{4\pi {\mu }_0}{\oint }_S\frac{d{A}^{\prime }}{\sqrt{(x-x^{\prime }{)}^2+(y-y^{\prime }{)}^2+(z-z^{\prime }{)}^2}}$$

Let us assume that the magnets have length $L$, width $w$, and depth $d$, where the length is measured along the $y$ axis, the width is measured along the $z$ axis, and the depth is measured along the $x$ axis. Thus, the rectangular cuboid we are integrating over is defined by $x\in [-d/2,d/2],y\in [0,L],z\in [0,w]$. To make this integral more tractable, we can assume that the magnets are “thin” magnets and thus set $d\to 0$, giving us the slightly “simpler” integral:

$${\psi }_m(\mathbf{r})=\frac{B_r}{4\pi {\mu }_0}{\oint }_S\frac{d{A}^{\prime }}{\sqrt{x^2+(y-y^{\prime }{)}^2+(z-z^{\prime }{)}^2}}$$

Where here we need only integrate over one surface (by symmetry, whether we choose this to be the front or back of the magnet doesn’t matter) as opposed to the six integrals that would otherwise need to be done for $d\ne 0$. Thus the surface integral reduces to the following double integral:

$${\psi }_m(\mathbf{r})=\frac{B_r}{4\pi {\mu }_0}{\int }_0^L{\int }_0^w\frac{d{z}^{\prime }d{y}^{\prime }}{\sqrt{x^2+(y-y^{\prime }{)}^2+(z-z^{\prime }{)}^2}}$$

To solve this integral, we first integrate with respect to $z^{\prime }$ with the following substitution:

$$\begin{array}{c}{\int }_0^L{\int }_0^w\frac{d{z}^{\prime }d{y}^{\prime }}{\sqrt{x^2+(y-y^{\prime }{)}^2+(z-z^{\prime }{)}^2}}={\int }_0^L{\int }_z^{z-w}\frac{du}{\sqrt{R^2+u^2}}d{y}^{\prime }\\ u=z-z^{\prime },R^2=x^2+(y-y^{\prime }{)}^2\end{array}$$

Now we may use the following identity:

$$\int \frac{dx}{\sqrt{a^2+x^2}}=\ln \left|x+\sqrt{x^2+a^2}\right|+C=\mathrm{arcsinh}\left(\frac{x}{a}\right)+C,a>0$$

Which allows us to evaluate the integral in $u$ as follows:

$$\begin{array}{rl}{\int }_0^L{\int }_z^{z-w}\frac{du}{\sqrt{R^2+u^2}}d{y}^{\prime }& ={\int }_0^L\left[\mathrm{arcsinh}\left(\frac{z-w}{R}\right)-\mathrm{arcsinh}\left(\frac{z}{R}\right)\right]d{y}^{\prime }\end{array}$$

Now, substituting in the explicit form of $R$, we have:

$$\begin{array}{rl}{\int }_0^L\left[\mathrm{arcsinh}\left(\frac{z-w}{R}\right)-\mathrm{arcsinh}\left(\frac{z}{R}\right)\right]d{y}^{\prime }& ={\int }_0^L\mathrm{arcsinh}\left(\frac{z-w}{\sqrt{x^2+(y-y^{\prime }{)}^2}}\right)d{y}^{\prime }\\ & -{\int }_0^L\mathrm{arcsinh}\left(\frac{z}{\sqrt{x^2+(y-y^{\prime }{)}^2}}\right)d{y}^{\prime }\end{array}$$

With the substitution $v=y-y^{\prime }$ we obtain the two integrals in a slightly simpler form:

$${\int }_y^{y-L}\mathrm{arcsinh}\left(\frac{z}{\sqrt{x^2+v^2}}\right)dv-{\int }_y^{y-L}\mathrm{arcsinh}\left(\frac{z-w}{\sqrt{x^2+v^2}}\right)dv$$

Both of the resulting integrals in the sum are integrals of the form:

$$\int \mathrm{arcsinh}\left(\frac{a}{\sqrt{b^2+x^2}}\right)dx$$

This is an integral that can be evaluated using integration by parts or (more conveniently) by consulting a computer algebra system; the solution is given (up to a constant) by³:

$$K(a,b,x)=a\ln \left(\left|\sqrt{x^2+b^2+a^2}+x\right|\right)-b\arctan \left(\frac{ax}{b\sqrt{x^2+b^2+a^2}}\right)+x\mathrm{arsinh}\left(\frac{a}{\sqrt{x^2+b^2}}\right)$$

Where we will not write out $K$ explicitly in the rest of the solution for simplicity. Thus we have:

$$\begin{array}{rl}{\int }_y^{y-L}\mathrm{arcsinh}\left(\frac{z}{\sqrt{x^2+v^2}}\right)dv& -{\int }_y^{y-L}\mathrm{arcsinh}\left(\frac{z-w}{\sqrt{x^2+v^2}}\right)dv\\ & =K(z,x,v){|}_{v=y}^{v=(y-L)}-K(z-w,x,v){|}_{v=y}^{v=(y-L)}\\ & =\mathcal{F}(\mathbf{r},L)\end{array}$$

Where we define $\mathcal{F}(\mathbf{r},L)=\mathcal{F}(x,y,z,L)$ as:

$$\begin{array}{rl}\mathcal{F}(x,y,z,L)& =K(z,x,v){|}_{v=y}^{v=(y-L)}-K(z-w,x,v){|}_{v=y}^{v=(y-L)}\\ & =K(z,x,y-L)-K(z,x,y)\\ & -K(z-w,x,y-L)+K(z-w,x,y)\end{array}$$

And thus the complete solution can then be written succinctly as:

$${\psi }_m(\mathbf{r})=\frac{B_r\mathcal{F}(\mathbf{r},L)}{4\pi {\mu }_0}$$

While appearing deceptively-simple, this succinct answer expands out to 12 terms(!) in total, making it highly impractical to use for actually compute this solution. A much more numerically-usable solution can be found by not evaluating the integral in $y^{\prime }$, giving us:

$${\psi }_m(\mathbf{r})=\frac{B_r}{4\pi {\mu }_0}(G(\mathbf{r},L,w)-G(\mathbf{r},L,0))$$

Where the function $G=G(\mathbf{r},h,\omega$ is defined by the following integral:

$$G(\mathbf{r},h,\omega )={\int }_0^h\mathrm{arcsinh}\left(\frac{z-\omega }{\sqrt{x^2+(y-s{)}^2}}\right)ds$$

And thus the magnetic field can be found (via $\mathbf{B}={\mu }_0\mathbf{H}=-{\mu }_0\nabla {\psi }_m$) to be:

$$\mathbf{B}=\frac{B_r}{4\pi }(\nabla G(\mathbf{r},L,0)-\nabla G(\mathbf{r},L,w))$$

The numerical derivatives here can be found by finite differences, but the analytical derivative can also be found by differentiating under the integral sign; that is:

$$\begin{array}{rl}\nabla G(\mathbf{r},h,\omega )& ={\int }_0^h\nabla \left[\left(\mathrm{arcsinh}\left(\frac{z-\omega }{\sqrt{x^2+(y-s{)}^2}}\right)\right)\right]ds\\ & ={\int }_0^h(z-\omega ){\left(1+\frac{(z-\omega {)}^2}{x^2+(y-s{)}^2}\right)}^{-1/2}\nabla \left(\frac{1}{\sqrt{x^2+(y-s{)}^2}}\right)\end{array}$$

Note that $\mathrm{arcsinh}(x)\approx x$, and thus for small $x$ and $y$ (that is, reasonably close to the magnet’s surface) we have:

$$\begin{array}{rl}G(\mathbf{r},h,\omega )\approx {\int }_0^h\frac{z-\omega }{\sqrt{x^2+(y-s{)}^2}}ds& =(z-\omega )\mathrm{arcsinh}\left(\frac{s-y}{x}\right){|}_0^h\\ & =(z-\omega )\left[\mathrm{arcsinh}\left(\frac{y}{x}\right)+\mathrm{arcsinh}\left(\frac{h-y}{x}\right)\right]\end{array}$$

If we then assume that $y$ is small, we can use the $\mathrm{arcsinh}(x)\approx x$ approximation once again, giving us:

$$G(\mathbf{r},h,\omega )\approx (z-\omega )\left[\mathrm{arcsinh}\left(\frac{y}{x}\right)+\mathrm{arcsinh}\left(\frac{h-y}{x}\right)\right]\approx (z-\omega )\left(\frac{y}{x}+\frac{h-y}{x}\right)=\frac{(z-\omega )h}{x}$$

Substituting this result into our analytical expression for ${\psi }_m$ gives us:

$$\begin{array}{rl}{\psi }_m(\mathbf{r})& =\frac{B_r}{4\pi {\mu }_0}(G(\mathbf{r},L,w)-G(\mathbf{r},L,0))\\ & =\frac{B_r}{4\pi {\mu }_0}\left[\frac{(z-w)L}{x}-\frac{zL}{x}\right]\\ & =-\frac{B_r}{4\pi {\mu }_0}\frac{wL}{x}\end{array}$$

From which we may find $\mathbf{B}={\mu }_0\mathbf{H}=-{\mu }_0\nabla {\psi }_m$ to be:

$$\mathbf{B}=\frac{B_{r}wL}{4\pi x^2}$$

A validation test for the solution is that when taking the other limit (for very far distances), the solution should reduce to a simple dipole:

$${\psi }_m(\mathbf{r})\approx \frac{\mathbf{m}\cdot \mathbf{r}}{4\pi r^3}$$

And therefore the magnetic field would fall off by the inverse cube of the distance. This is what allows permanent magnets to be well-approximated by perfect magnetic dipoles.

TODO: Make a 3D plot of the exact solution and the field lines for the $\mathbf{H}$ and $\mathbf{B}$ fields.

Solving with a CAS

In theory, the surface integral can also be evaluated with the following Mathematica code:

yprime = Subscript[p, y]
zprime = Subscript[p, z]
Integrate[
	 1/Sqrt[x^2 + (y - yprime)^2 + (z - zprime)^2],
	 {yprime, 0, L}, {zprime, 0, w}
 ]

Or with the following SymPy code:

from sympy import *
x, y, z = symbols("x y z", real=True)
Br, w, L = symbols("B_r w L", constant=True, 
		real=True, positive=True)
mu0 = Symbol(r"\mu_0")
yprime, zprime = symbols("y' z'")
prefactor = Br/(2*pi*mu0)
integrand = 1/sqrt(x**2 + (y-yprime)**2 + (z - zprime)**2)
integral = Integral(integrand, (yprime, 0, L), (zprime, 0, w))
psi = prefactor * integral.doit()
Eq(Symbol(r"\psi_m"), psi)

The issue is that both SymPy and Mathematica time out when trying to compute the integral; it is likely that neither can compute it, even if they are given longer times.

Interior and exterior case

We will now consider the more complicated case where we solve for both the interior and exterior magnetic fields. This involves a more complicated boundary-value problem that generally cannot be solved analytically, even in integral form. Instead, numerical methods are necessary for most cases.

To formulate the problem, we must consider the pertinent boundary conditions. First, we impose the boundary condition that the magnetic scalar potential vanish at infinity, that is:

$$\underset{\mathbf{r}\to \infty }{\lim }{\psi }_m(\mathbf{r})=0$$

In practice, this boundary condition is only useful analytically; simulating an infinitely-large domain numerically is impractical (though there are some numerical methods that can do so, such as the boundary element method), so we can simply set a cutoff radius $R$ (which is large compared to the size of the magnet), giving us:

$${\psi }_m(R)\approx 0$$

Next, we consider the boundary conditions at the interface between the interior of the magnet (with permeability $\mu$) and the exterior of the magnet (with permeability $\approx {\mu }_0$). When the relative permeability ${\mu }_r=\mu /{\mu }_0\to \infty$ (which is approximately true for many ferromagnetic materials), then just outside the surface of the magnet⁴, $H_{\parallel }$ (the tangential component) becomes zero and $H_{\perp }$ (the normal component) becomes the only component of $\mathbf{H}$; this is exactly analogous to the perfect conductor in electrostatics. Therefore, at the surface of the magnet, we obtain the following Neumann boundary condition:

$$H_{\parallel }=-\frac{\partial {\psi }_m}{\partial \tau }{|}_{\partial \Sigma }=0$$

Where $\partial \Sigma$ denotes the surface of the magnet, $\tau$ denotes the coordinate tangent to the surface, and where the solution is valid for the entirety of the exterior of the magnet. Meanwhile, we may determine $H_{\perp }$ using the standard boundary conditions for magnetostatics, given by:

$$\mathbf{K}=\mathbf{M}\times \hat{\mathbf{n}},\hat{\mathbf{n}}\times ({\mathbf{H}}^{(2)}-{\mathbf{H}}^{(1)})=\mathbf{K}$$

Where $\mathbf{K}$ is the effective surface current density, ${\mathbf{H}}^{(1)}$ is the field just underneath the surface of the magnet, and ${\mathbf{H}}^{(2)}$ is the field just outside the magnet (assumed to be vacuum or another non-permeable medium). Then by substitution, one has:

$$\begin{array}{c}\hat{\mathbf{n}}\times ({\mathbf{H}}^{(2)}-{\mathbf{H}}^{(1)})=-(\hat{\mathbf{n}}\times \mathbf{M})\\ {\mathbf{H}}^{(1)}-{\mathbf{H}}^{(2)}=\mathbf{M}\end{array}$$

More precisely, the above is only true if ${\mathbf{H}}^{(1)},{\mathbf{H}}^{(2)}$ are nearly normal to the surface (which occurs in the limit of very high relative permeabilities). In the general case, we have:

$${\mathbf{H}}_{\perp }^{(1)}-{\mathbf{H}}_{\perp }^{(2)}=\mathbf{M}$$

Where ${\mathbf{H}}_{\perp }^{(1)},{\mathbf{H}}_{\perp }^{(2)}$ are the normal components of the ${\mathbf{H}}^{(1)},{\mathbf{H}}^{(2)}$ fields, respectively. Now, to determine ${\mathbf{H}}^{(1)}$, we use the constitutive relation $\mathbf{B}={\mu }_0(\mathbf{M}+\mathbf{H})$. Here it is implied that $\mathbf{M},\mathbf{H},\mathbf{B}$ represent the interior fields and magnetizations; that is, $\mathbf{M}={\mathbf{M}}_{interior}$, $\mathbf{H}={\mathbf{H}}^{(1)}$, $\mathbf{B}={\mathbf{B}}_1$. We must now impose a constitutive relation between $\mathbf{B}$ and $\mathbf{H}$ to solve. In general, this relationship is nonlinear for permanent magnets, and thus $\mathbf{B}=f(\mathbf{H})$. However, for many permanent magnets⁵ we can make the linear approximation given by⁶:

$$\mathbf{B}={\mathbf{B}}_r+{\mu }_0{\mu }_{rec}\mathbf{H}$$

Here, ${\mu }_{rec}$ is the recoil permeability of the material and is a material constant, where we can often make the assumption ${\mu }_{rec}\approx 1$ (which is true for a perfect magnet). From these boundary conditions, we can then solve for the fields both inside and outside the magnet.

Magnetic fields for the magnetic lens

The field of a sextupole magnet can be written in spherical coordinates (assuming azimuthal symmetry, or in this case, $\varphi =0$) via the multipole expansion, which, up to some constant factors, takes the form:

$$B(r,\theta )=B_0{P}_0(\cos \theta )+B_1{P}_1(\cos \theta )r+B_2{P}_2(\cos \theta )r^2+\dots$$

Where $B_0,B_1,B_2,\dots$ are some constant factors, and $P_n(\cos \theta )$ is a Legendre polynomial. The first three Legendre polynomials are given by:

$$\begin{array}{rl}{P}_0(x)& =1\\ P_1(x)& =x\\ P_2(x)& =\frac{1}{2}(3x^2-1)\end{array}$$

Thus, we can write $B(r,\theta )$ in explicit form as:

$$B(r,\theta )=B_0+B_{1}r\cos \theta +\frac{1}{2}{B}_2(3{\cos }^2\theta -1)r^2+\dots$$

Note: For the purposes of our discussion, the precise forms of $B_0,B_1,B_2$ don’t matter, since we’re only interested in the qualitative characteristics of the solution, which are not affected by the values of the constants.

To describe a sextupole magnet, we truncate the series after three terms, hence the name “sextupole” (which comes after the dipole (first) term and the quadrupole (second) term). An interactive Desmos visualization shows the field magnitude as a function of angle. The strong fields near the edge of the sextupole magnet are what cause charges to be repelled towards the center of the magnet, which allows the sextupole magnet to serve the function of a magnetic lens. Borrowing an idea from transmission electron microscopes⁶, if we use two sextupole magnets of opposite polarities (that is, one flipped so all N-S pairs become S-N pairs), the two lenses can act as a combined converging-diverging lens. This means that rather than diverging the beam (as with a single magnet lens), the two magnets both collimate the beam and restore the beam back to a circular shape.

Numerical methods

The typical numerical method used for magnetostatics is finite-element analysis,. To begin, we start with the basic differential equation of permanent magnets⁷:

$$\nabla \cdot (\mu (-\nabla {\psi }_m)+\mathbf{M})=0$$

Assuming a constant permeability $\mu$, one therefore obtains:

$$-\mu {\nabla }^2{\psi }_m+\nabla \cdot \mathbf{M}=0$$

To derive the weak form, we multiply by test function $v$, and integrate both sides (see a more in-depth description of the finite element method in RF cavity simulation with aperture). This gives us:

$$\begin{array}{r}{\int }_{\Omega }v(-\mu {\nabla }^2{\psi }_m+\nabla \cdot \mathbf{M})dV=0\\ -\mu {\int }_{\Omega }v{\nabla }^2{\psi }_{m}dV+{\int }_{\Omega }v\nabla \cdot \mathbf{M}=0\end{array}$$

We can simplify this by using integration by parts or by simply consulting a list of vector calculus integral identities. This gives us:

$$\mu {\int }_{\Omega }\nabla v\cdot \nabla {\psi }_{m}dV-\mu {\int }_{\partial \Omega }v\nabla {\psi }_m\cdot d\vec{S}-{\int }_{\Omega }\mathbf{M}\cdot \nabla vdV+{\int }_{\partial \Omega }v\mathbf{M}\cdot d\vec{S}=0$$

Combining the volume and surface terms gives us:

$$\mu {\int }_{\Omega }(\nabla v\cdot \nabla {\psi }_m-\mathbf{M}\cdot \nabla v)dV-\mu {\int }_{\partial \Omega }v(\nabla {\psi }_m-\mathbf{M})\cdot d\vec{S}=0$$

Note that $\mathbf{M}$ must be pre-provided. In our case for permanent magnets (as detailed previously), this is a piecewise function:

$$\mathbf{M}=\{\begin{array}{ll}\frac{B_r}{{\mu }_0}\hat{\mathbf{z}},& \text{inside magnet}\\ 0,& \text{outside magnet}\end{array}$$

Where $B_r$ is the remanence, as before. This can then be solved for suitable boundary conditions (as described above). There are dedicated open-source magnetostatics codes such as Radia as well as general-purpose FEM packages (such as FreeFEM and MFEM) that can be used to solve the resulting BVP numerically.

Footnotes

1. See Jackson, Classical Electrodynamics (3rd. ed.), pg. 207, eq. 5.100. We assume a permanent magnet with homogeneous magnetization, and thus the volume integral term drops out (since $\nabla \cdot \mathbf{M}=0$). ↩

2. Please see Jiles, Introduction to Magnetism and Magnetic Materials (2nd ed.), 1998. Originally from this article. ↩

3. This was obtained using MathDF but can also be found using a dedicated computer algebra system. ↩

4. This applies locally, no matter the geometry of the magnet itself; for more information, see Jackson, Classical Electrodynamics (3rd. ed.), pg. 212 ↩

5. Specifically for “hard” magnets like neodymium or ferrite magnets, which keep their magnetization near-permanently, independent of any applied fields. ↩

6. From Campbell, Permanent Magnet Materials and Their Application, 1994 (p. 94-96). Originally from this Physics SE answer. ↩ ↩²

7. Referenced from the Wolfram guide on magnetostatic modelling for permanent magnets. ↩