Numerical modelling of electron beams

Free-electron lasers and masers use a relativistic electron beam as their light source, meaning that being able to study and understand the behavior of free-electron lasers requires being able to model electron beams. Here, we give a description of the mathematics and particularly the numerical methods used for modelling electron beams, electron guns, undulators, and other related systems.

Note: Also see Guide to modelling lasers and masers for more information about free-electron lasers/masers (as well as other types of lasers/masers).

Modelling of undulator electron beams

Here, we will use slightly different coordinate conventions compared to Free-electron maser physics. In our new conventions, $z$ is the optical axis, $y$ is the direction parallel to the north-south axis of the magnets, and $x$ the direction perpendicular to the magnets. The electrons travel in the $xz$ plane in the undulator.

The electron beam trajectory $\mathbf{r}(t)$ follows the relativistic Lorentz force equation for magnetic fields:

$$\frac{d\mathbf{p}}{dt}=\frac{d}{dt}(\gamma m\dot{\mathbf{r}})=q\mathbf{v}\times \mathbf{B}$$

This is a highly-nonlinear differential equation, one we will need to solve (either through analytical or numerical means) to determine the trajectory of the electron beam through undulator.

A heuristic derivation of the beam energy threshold

A good way to try to estimate the beam energy necessary is to consider the following case. Let $\mathbf{B}=B_0\hat{z}$ be a uniform field oriented transverse to the motion of a charged particle. By Newton’s second law, we know that the particle will follow a circular path if its velocity $\mathbf{v}$ satisfies:

$$\frac{m{v}^2}{R}=qv{B}_0\sin \theta ,v=|\mathbf{v}|$$

Thus upon rearrangement we have:

$$v=\frac{q{B}_0}{m}R\sin \theta$$

Therefore, only a particle with an “escape velocity” $v_{esc}>v$ would not end up in a circular orbit. Note that we can rearrange in terms of $\beta =v/c$ as follows:

$$\beta =\frac{q{B}_0}{mc}R\sin \theta$$

For an electron, we can rearrange in terms of natural constants $e,m_e$ to obtain:

$${\beta }_{min}=-\frac{e{B}_0}{m_{e}c}R\sin \theta =-\frac{e{B}_{0}c}{E_{electron}}R\sin \theta$$

Where $E_{mass}=0.511\mathrm{MeV}$ is the electron mass-energy. Since the electron beam should be fairly paraxial, we know that $R/w\ll 1$ (where $w$ is the width of the undulator), which allows us to obtain a rough order-of-magnitude estimate for $R$ (in addition, we can take known values of $R$ for other free-electron maser and laser installations). At a $0.1\mathrm{T}$ field and a $1\mathrm{mm}$ radius beam, this yields a minimum velocity of $\beta >0.058$ (or around 6% of the speed of light, corresponding to a beam energy of around 1 keV) for $\theta =0$. Note, however, that this is less for off-axis electrons, i.e. for $\theta >0$. While this calculation is ultimately a very rough approximation, it does set a theoretical minimum for the energy necessary for an electron to escape the rotational pull of the magnetic field.

We may obtain a more accurate approximation from $a=\frac{{\omega }_0^{\prime }}{k_u^{2}c}$, the analytical expression for the displacement of an electron beam from the centerline in an ideal undulator (from Jackson, Classical Electrodynamics), where ${\omega }_0^{\prime }\equiv {\omega }_0/\gamma$, $k_u=2\pi /L_u$, and ${\omega }_0=e{B}_0/m_e$. Thus after the substitution $R\to a$, one has:

$$\begin{array}{rl}{\beta }_{min}& =-\frac{{\omega }_0^{\prime }}{k_u^{2}c}\frac{e{B}_{0}c}{E_{electron}}\sin \theta \\ & =-\frac{(e{B}_0{L}_u{)}^2}{m_e^2{c}^2}\frac{\sin \theta }{4{\pi }^2\gamma }\\ & \approx (8718.5{\mathrm{C}}^2{\mathrm{s}}^2\mathrm{k}{\mathrm{g}}^{-2}{\mathrm{m}}^{-2})\times \frac{(B_0{L}_u{)}^2\sin \theta }{\gamma }\end{array}$$

This result gives around 1% the speed of light for $L_u=1\mathrm{cm}$, and around 8% the speed of light for $L_u=3\mathrm{cm}$, rising to around 14% the speed of light for $L_u=4.5\mathrm{cm}$ due to the quadratic dependence on $L_u$.

Approximate solution in the highly-relativistic regime

The analytical solution to the problem is based on the famous EM textbook Classical Electrodynamics by Jackson. On pg. 684 (in chapter 14.7, 3rd. ed.), the solution for the beam trajectories is given as $x\approx a\sin (2\pi z/{\lambda }_0)$ (here Jackson uses $k_0\equiv 2\pi /{\lambda }_0$ but we would write this as $x=a\sin (k_{u}z)$ in our notation). He then claims that one can divide between wigglers (which emit incoherent radiation) and undulators (which emit coherent radiation) via the value of $K$. For a wiggler, $K\gg 1$, but for an undulator, $K\ll 1$. This is a result that we have regarded as implicitly-known.

On pg. 687, Jackson then defines $K=\gamma k_{0}a$ (in our notation, $K=\gamma k_{u}a$) and $K=\frac{e{B}_0{L}_u}{2\pi m{c}^2}$. This definition of $K$ is clearly in Gaussian units due to his extra factor of $c$, since in SI we have $K=\frac{e{B}_0{L}_u}{2\pi mc}$, which we can also write in our notation as $K=\frac{{\omega }_0}{k_{u}c}$. Rearranging gives us $a=K/\gamma k_u=\frac{{\omega }_0}{\gamma c{k}_u^2}$. If we then define ${\omega }_0^{\prime }\equiv {\omega }_0/\gamma$ we can write this more elegantly as $a=\frac{{\omega }_0^{\prime }}{k_u^{2}c}$. Thus (in our notation):

$$x\approx a\sin (k_{u}z)=\frac{{\omega }_0^{\prime }}{k_u^{2}c}\sin (k_{u}z)$$

Thus we see that the electrons exhibit oscillatory “wiggling” motion with period $L_u$, and a maximum displacement proportional to ${\omega }_0^{\prime }$. Since ${\omega }_0^{\prime }\sim {\gamma }^{-1}$, the maximum displacement of the electron oscillations decreases for faster electrons, but the period remains constant, which is the result of relativistic length contraction.

Note: For the electrons to not smash into the undulator walls, we obviously need $a<w/2$, where $w$ is the width of the undulator. Note that $a\propto B_0{L}_u^2/\gamma$ so the maximum displacement of the electrons from the $z$ axis is strongly dependent on the separation between the magnets. In practice, even with very strong fields of $B_0=\text{1 T}$ and a long magnet separation of $L_u=\text{5 cm}$, $a$ is only around 3.7 cm.

However, exactly how Jackson’s initial claim that $x\approx a\sin (2\pi z/{\lambda }_0)$ is derived is not directly shown; we can assume that it is just an ansatz. A hint comes from when Jackson states that:

“If the periodicity of the magnetic field structure is ${\lambda }_0$, the particle’s path will be approximately sinusoidal in the transverse direction with the same period.” (pg. 684)

He does explain that (on pg. 686) that working from the Lorentz force equation the $x$-component can be written as $\ddot{x}=-e{\beta }_z{B}_y/\gamma m$, which is the only meaningful component since ${\beta }_y{B}_z$ is assumed to be “negligible or zero”. Note that since he is working in Gaussian units, the relativistic Lorentz force equation reads as $\frac{d}{dt}(\gamma v)=\frac{q}{mc}\vec{v}\times \vec{B}=\frac{q}{m}\vec{\beta }\times \vec{B}$, which has an extra factor of $1/c$ compared to SI units. He assumes that $|\beta |\approx \text{const.}$ (its components may change but its magnitude stays relatively constant) so $\gamma \approx \text{const.}$ and he can factor it out, giving $\dot{v}=\frac{q}{\gamma m}\vec{\beta }\times \vec{B}$.

With the approximations $z\approx ct$ and ${\beta }_z\approx 1$, he arrives at an expression for the magnetic field, where $B_y\sim \sin (k_{0}z)$. Afterwards, he acknowledges that “an actual magnet structure will be periodic, but not sinusoidal” and that $B_y$ is better represented as a Fourier series in $k_0$, though he keeps only the fundamental mode of the series for simplicity.

On the following page, he gives ${\beta }_x\approx k_{0}a\cos (k_{0}ct)$ (which can be written as ${\beta }_x\approx k_{0}a\cos (k_{0}z)$ using his previous approximation), so at $t=0$ we have ${\beta }_x(0)=k_{0}a$. And since it is assumed that ${\beta }_y=0$ then ${\beta }_z^2={\beta }^2-{\beta }_x^2$, where he claims that $|{\beta }_x|\ll {\beta }_z$ and thus he uses a binomial approximation to derive ${\beta }_z$ (specifically, ${\beta }_z\approx |\beta |-\frac{1}{2}{\beta }_x^2$). We may infer that it is reasonable to assume that the electrons are initially moving nearly purely along the $z$ axis as they enter the undulator (by definition the direction of motion points in the same direction as the velocity vector, and thus points along $\vec{\beta }$).

However, note that the solution from Jackson is only valid under a wide range of simplifying assumptions, including (but not limited to) a magnetic field that is close to a pure harmonic, highly-relativistic electrons, $\beta \approx \text{const.}$, and $K\ll 1$. These are assumptions that generally do not hold in the mildly-relativistic regime that free-electron masers typically operate at, so we need better approximations.

Derivation of Jackson’s analytical solution

Since Jackson does not derive his analytical solution for $x(z)$ directly, we will sketch out a derivation here. To start, we know (or at least Jackson tells us) that the magnetic field inside the undulator is given by:

$$B_x=B_z=0,B_y=B_0\sin k_{u}x$$

(Note that here $k_u$ is the same thing as $k_0$; it is just a different notation that we’ll use interchangeably). Here, $B_0$ is the magnetic field amplitude, $k_u=2\pi /L_u$ is the wavenumber, and $L_u$ is the spatial period (distance between every pair of magnets with same polarity).

To solve the problem approximately, we assume that the initial velocity $v_0$ of the electrons (when they are fed into the undulator) is almost entirely along the $z$ axis. Since there is no $B_x$ or $B_z$ component, by substituting and expanding the cross-products in the Lorentz force formula, we have:

$$\begin{array}{r}\frac{d}{dt}(\gamma m\dot{x})=q({\cancel{v_y{B}_z}}^0-v_z{B}_y)\\ \frac{d}{dt}(\gamma m\dot{y})=q({\cancel{v_z{B}_x}}^0-\cancel{v_x{B}_z})\\ \frac{d}{dt}(\gamma m\dot{z})=q(v_x{B}_y-{\cancel{v_y{B}_x}}^0)\end{array}$$

Thus, we find that the electrons move purely in the $xz$-plane, and their equations of motion are given by:

$$\begin{array}{rl}\frac{d}{dt}(\gamma m\dot{x})& =-q{v}_z{B}_y\\ \frac{d}{dt}(\gamma m\dot{z})& =q{v}_x{B}_y\end{array}$$

Now, we will make the assumption that the $z$-component of the electron’s velocity is assumed to be essentially constant ($v_z\approx v_0$). This means that the second differential equation becomes $\frac{d}{dt}(\gamma m\dot{z})=\frac{d}{dt}(\gamma m{v}_z)=\frac{d}{dt}(\gamma m{v}_0)=0$, which has the trivial solution of $z=v_{0}t$. It also means that we have just one differential equation left, that being the one for $x$. But given that we know $v_z$, the first differential equation simplifies to:

$$\frac{d}{dt}(\gamma m\dot{x})=-q{v}_0{B}_y$$

This is still a nonlinear differential equation, and indeed if we expand it it takes the rather-ugly-looking form¹:

$$\dot{\gamma }\frac{dx}{dt}-\gamma \frac{d^{2}x}{d{t}^2}={\omega }_0{v}_0\sin (k_u{v}_{0}t),{\omega }_0\equiv q{B}_0/m$$

Where we once again define ${\omega }_0\equiv q{B}_0/m$. This is made even worse by the fact that the Lorentz factor $\gamma$ itself depends on $\dot{x}$, making this even harder to solve. It certainly cannot be solved analytically in this form. However, if we assume that $\gamma \approx \text{const.}$, then it can be factored away, giving us:

$$\frac{d^{2}x}{d{t}^2}=-({\omega }_0{v}_0/\gamma )\sin (k_u{v}_{0}t)$$

Solving is then simply a matter of twice-integrating the above, giving us:

$$x(t)=-\frac{{\omega }_0}{\gamma v_0{k}_u^2}\sin k_u{v}_{0}t=-\frac{{\omega }_0}{\gamma v_0{k}_u^2}\sin (k_{u}z)$$

Defining ${\omega }_0^{\prime }={\omega }_0/\gamma$ as the “relativistic” form of ${\omega }_0$ gives us:

$$x=-\frac{{\omega }_0^{\prime }}{v_0{k}_u^2}\sin k_{u}z=a\sin k_{u}z,a=-\frac{{\omega }_0^{\prime }}{v_0{k}_u^2}$$

Note that with the approximation $v_0\approx c$ (the strongly-relativistic limit) this is the same result as Jackson (where $a=\frac{{\omega }_0^{\prime }}{k_u^{2}c}$), up to a sign difference; the sign difference comes from the negative charge of the electrons, so ${\omega }_0\to -{\omega }_0$. This result clearly shows that the electrons move in sinusoidal trajectories, “wiggling” from left to right along the $xz$ plane - which is indeed what we would expect. However, remember that this is an approximate solution!

Approximate solution in the non-relativistic and mildly-relativistic regime

A more rigorous approximate solution can be found by linearizing the system about the critical points. This solution works in the non-relativistic and mildly-relativistic limit, where we can assume that $\gamma \approx 1$. However, it is not too difficult to extent it to the strongly-relativistic limit as long as we assume that $\gamma \approx \text{const}$.

To start, let us begin with the exact equations of motion with no simplifications:

$$\begin{array}{rl}\frac{d}{dt}(\gamma m\dot{x})& =-q\dot{z}{B}_y\\ \frac{d}{dt}(\gamma m\dot{z})& =q\dot{x}{B}_y\end{array}$$

As a reminder, these are highly-nonlinear differential equations that are not easy to solve using SciPy’s native differential equation solver. However, we can quasi-linearize them by taking $\gamma \approx 1$ which holds in the non-relativistic and mildly-relativistic limit (meaning it is valid to around $v=0.5c$, more than enough for our purposes). This gives us:

$$\begin{array}{rl}\ddot{x}& =-\dot{z}q{B}_y/m\\ \ddot{z}& =\dot{x}q{B}_y/m\end{array}$$

Or in vector form:

$$\frac{d}{dt}\left(\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right)=\frac{q{B}_y}{m}\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right)$$

Note that if $B_y=\text{const.}$ (we may call this constant field $B_0$) then it would be possible to solve this system of differential equations analytically. We would then have $\frac{q{B}_y}{m}=\frac{q{B}_0}{m}={\omega }_0$, and the system simplifies to:

$$\frac{d}{dt}\left(\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right)={\omega }_0\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right)$$

Using a computer algebra system or by just using the general methods of solving ODEs, we can solve the (simplified) system exactly; the solution is simply in terms of sines and cosines. An example implementation for solving the system symbolically using SymPy is given below:

# first show symbolic solution
import sympy as sp
from sympy import symbols, Eq, Function, re, Matrix
from sympy.solvers.ode.systems import dsolve_system
t, q, By, m = symbols("t q B_0 m", real=True)
x = Function("x")(t)
z = Function("z")(t)
 
# define derivatives
xdot = x.diff(t)
xddot = xdot.diff(t)
zdot = z.diff(t)
zddot = zdot.diff(t)
w0 = q*By/m
 
# define the simplified system of ODEs
ode1 = Eq(xddot, -w0*zdot)
ode2 = Eq(zddot, w0*xdot)
 
# solve the simplified system
symbolic_sols = dsolve_system([ode1, ode2])
 
# print out the solutions from SymPy
symbolic_sols[0][0]
symbolic_sols[0][1]

The solutions (from SymPy or worked-out manually) are given by:

$$\begin{array}{rl}x(t)& =C_1-\frac{C_2}{{\omega }_0}\sin ({\omega }_{0}t)+\frac{C_3}{{\omega }_0}\cos ({\omega }_{0}t)\\ z(t)& =C_4+\frac{C_2}{{\omega }_0}\cos ({\omega }_{0}t)+\frac{C_3}{{\omega }_0}\cos ({\omega }_{0}t)\end{array}$$

Thus we find that if $B_y=B_0=\text{const.}$ (that is, if it were a uniform field) we have $x(t)\sim \sin ({\omega }_{0}t),z(t)\sim \cos ({\omega }_{0}t)$ where here ${\omega }_0=B_{0}q/m$ is an angular frequency. Therefore, in the case of a uniform field (or alternatively, in the limit of $L_u\to 0$), we have formally shown the electrons just orbit in circles with characteristic orbital frequency ${\omega }_0$.

We can also use a classical method for solving nonlinear ODEs to obtain an approximate solution without assuming $B_y=\text{const}$. This involves linearizing the system about its steady-states. To do so, let us start with the system with only the simplification $\gamma \approx 1$. This gives us:

$$\left(\begin{array}{c}\ddot{x}\\ \ddot{z}\end{array}\right)=\sin (k_{u}z)\left(\begin{array}{cc}0& -{\omega }_0\\ {\omega }_0& 0\end{array}\right)\left(\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right)$$

The system can then be rewritten in the following (phase-space) form, where $\mathbf{y}=(x,z,v_x,v_z)$:

$$\dot{\mathbf{y}}=\left(\begin{array}{c}\dot{x}\\ \dot{z}\\ {\dot{v}}_x\\ {\dot{v}}_z\end{array}\right)=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 0& -{\omega }_0\sin (k_{u}z)\\ 0& 0& {\omega }_0\sin (k_{u}z)& 0\end{array}\right)\left(\begin{array}{c}x\\ z\\ v_x\\ v_z\end{array}\right)$$

The steady-states of the system can be found by setting the RHS of the above equation to zero; they are at $z=0$ and $k_{u}z=\frac{n\pi }{2}$ (for odd $n$), with $v_x=v_z=0$ in both cases. Since the $z=0$ steady-state is trivial, the first non-trivial steady state occurs at $z^{\ast }=\frac{\pi }{2k_u}$. Now expanding about this steady state is equivalent to evaluating the Jacobian at $z^{\ast }=\frac{\pi }{2k_u}$. The Jacobian in this case is given (in general) by:

$$J(x,z,v_x,v_z)=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& -k_u{\omega }_0\cos (k_{u}z)v_z& 0& -{\omega }_0\sin (k_{u}z)\\ 0& k_u{\omega }_0\cos (k_{u}z)v_x& {\omega }_0\sin (k_{u}z)& 0\end{array}\right)$$

Evaluating at $z^{\ast }=\pi /2k_u$ gives us:

$$J(z^{\ast })=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 0& -{\omega }_0\\ 0& 0& {\omega }_0& 0\end{array}\right)$$

The linearized version of the differential equation is then given by $\dot{\mathbf{y}}=J(z^{\ast })\mathbf{y}$. The solution can be found by first finding the eigenvalues and eigenvectors of the Jacobian $J(z^{\ast })$ (which we’ll notate as $J^{\ast }$ from now on for short). They can be manually, or (more conveniently) with a computer algebra system; an example code using SymPy is shown below:

from sympy import symbols, Matrix
# calculate eigenvalues & eigenvecs of Jacobian
w0 = symbols("w_0")
J_star = Matrix([
    [0, 0, 1, 0],
    [0, 0, 0, 1],
    [0, 0, 0, -w0],
    [0, 0, w0, 0]
])
J_eigenvecs = J_star.eigenvects()
print(J_eigenvecs)

Whether calculated manually or on computer, the result should be the same: there are three sets of eigenvalues and eigenvectors. The first is the trivial eigenvalue $\lambda =0$ with multiplicity 2, with associated eigenvectors $(1,0,0,0{)}^T$ and $(0,1,0,0{)}^T$. The second is the eigenvalue $-i{\omega }_0$ with multiplicity 1, with associated eigenvector $(1,i,-i{\omega }_0,{\omega }_0{)}^T$. The third is simply the conjugate of the second, with the eigenvalue of $i{\omega }_0$ and associate eigenvector $(1,-i,i{\omega }_0,w_0)$. Thus the general solutions are:

$$\left[\begin{array}{c}x(t)\\ z(t)\\ v_x(t)\\ v_z(t)\end{array}\right]=\left[\begin{array}{c}{C}_1\\ C_2\\ 0\\ 0\end{array}\right]+C_3\left[\begin{array}{c}1\\ i\\ -i{\omega }_0\\ {\omega }_0\end{array}\right]e^{-i{\omega }_{0}t}+C_4\left[\begin{array}{c}1\\ -i\\ i{\omega }_0\\ {\omega }_0\end{array}\right]e^{i{\omega }_{0}t}$$

Simplifying with some new constants ${\tilde{C}}_3,{\tilde{C}}_4$ gives us:

$$\begin{array}{r}x(t)=C_1+{\tilde{C}}_3\cos {\omega }_{0}t\\ z(t)=C_2+{\tilde{C}}_4\sin {\omega }_{0}t\end{array}$$

Which…is a fairly underwhelming result, since we have essentially arrived at our circling-particle solution once again, although this time in a much more mathematically-formal manner. However, this result is not useless; an important aspect of any numerical solver is that it must reproduce this solution in the appropriate limits.

Numerical solution of electron beam trajectories

To solve the highly-nonlinear ODE system, we need to use numerical methods for solving differential equations. While there are commercial software that can be used like EGUN, we will describe a general treatment that is (mostly) software-agnostic. For example code, we will use Python, but any other programming language works fine.

To begin, we will start with the exact, relativistic Lorentz force equations:

$$\begin{array}{rl}\frac{d}{dt}(\gamma m\dot{x})& =-q\dot{z}{B}_y\\ \frac{d}{dt}(\gamma m\dot{z})& =q\dot{x}{B}_y\end{array}$$

The presence of the Lorentz factor $\gamma$ makes them incredibly difficult to solve. However, if we take $\gamma \approx 1$ which holds in the non-relativistic and mildly-relativistic limit (meaning it is valid to around $v=0.5c$, more than enough for our purposes), we have:

$$\begin{array}{rl}\ddot{x}& =-\dot{z}q{B}_y/m\\ \ddot{z}& =\dot{x}q{B}_y/m\end{array}$$

Or in vector form:

$$\frac{d}{dt}\left(\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right)=\frac{q{B}_y}{m}\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right)$$

Since the electron charge and mass are so small while the speed of light is so big (in SI units), substituting in the raw values will lead to numerical overflows and break the solvers. To avoid this we will perform nondimensionalization. To begin, the (non-relativistic) Lorentz force law with zero electric field can be written in the form:

$${\ddot{x}}^i=\frac{q{B}_{\perp }}{m}{\dot{x}}^j={\omega }_{0}f(x^k){\dot{x}}^j,B_{\perp }=B_{0}f(x^k)$$

Where ${\dot{x}}^j$ is the $j$-th component of the velocity vector and ${\ddot{x}}^i$ is the $i$-th component of the acceleration vector (in tensor notation), and the two are always perpendicular to each other. Now, we will make a change of variables by defining two new variables $\chi ,\tau$, where $x=\chi x_c$ and $t=\tau t_c$. Then:

$$\frac{d}{dt}=\frac{d\tau }{dt}\frac{d}{d\tau }=\frac{1}{t_c}\frac{d}{d\tau },\frac{d{x}^j}{dt}=\frac{x_c}{t_c}\frac{d{\chi }^j}{d\tau }$$

Likewise:

$$\frac{d^2{x}^i}{d{t}^2}=\frac{d}{dt}\left(\frac{d{x}^i}{dt}\right)=\frac{x_c}{t_c^2}\frac{d^2{\chi }^i}{d{\tau }^2}$$

Substituting back into the Lorentz force law gives us:

$$\frac{x_c}{t_c^2}\frac{d^2{\chi }^i}{d{\tau }^2}={\omega }_{0}f(x_c{\chi }^k)\frac{x_c}{t_c}\frac{d{\chi }^j}{d\tau }$$

Simplifying and rearranging, we have:

$$x_c\frac{d^2{\chi }^i}{d{\tau }^2}=t_c{\omega }_0{x}_{c}f(x_c{\chi }^k)\frac{d{\chi }^j}{d\tau }$$

If we make the identification $t_c=1/{\omega }_0$ and define $\tilde{f}(\xi )=f(x_c\xi )$, we have thus nondimensionalized the differential equation:

$$\frac{d^2{\chi }^i}{d{\tau }^2}=\tilde{f}({\chi }^k)\frac{d{\chi }^j}{d\tau }$$

Note: The choice for $x_c$ is somewhat arbitrary since it is divided from both sides. This freedom in specifying $x_c$ will be important soon.

In its raw form, this differential equation can be numerically-unstable for very fast-moving electrons, since the electrons move so fast that $\frac{d{\chi }^j}{d\tau }$ blows up, leading to numerical overflows and errors. Thus, we need to perform some manipulations to make the differential equation stable enough for numerical integration. First, note that by multiplying both sides by $1/c$ and defining ${\tilde{\beta }}^i=\frac{1}{c}\frac{d{\chi }^i}{d\tau }$, one may obtain an equivalent differential equation, though now expressed in terms of ${\chi }^k$ and ${\tilde{\beta }}^k$:

$$\begin{array}{rl}\frac{1}{c}\frac{d^2{\chi }^i}{d{\tau }^2}& =\frac{1}{c}\tilde{f}({\chi }^k)\frac{d{\chi }^j}{d\tau }\\ \frac{d{\tilde{\beta }}^i}{d\tau }& =\tilde{f}({\chi }^k){\tilde{\beta }}^j\end{array}$$

This has the advantage of keeping our velocities expressed as reasonably small values, instead of being massive numbers on the order of $v\sim 10^8$! For electron beams with electrons travelling at a substantial percentage of the speed of light, which is standard for free-electron lasers and masers, this is indeed essential.

Note that this expression is very general since any magnetic field can be written in the form $\mathbf{B}={\mathbf{B}}_{0}f(\mathbf{x})$, and we need only perform a simple transformation $\tilde{f}(\xi )=f(x_c\xi )$ to obtain $\tilde{f}$ from $f$.

We will now substitute our particular case, where the only nonzero component of the field is $B_{\perp }=B_y$ and we have $B_y=B_0\sin (k_{u}z)$, and thus $f({\chi }^k)=\sin (k_u{\chi }^k)$ and $\tilde{f}({\chi }^k)=\sin (x_c{k}_u{\chi }^k)$. Thus in our case, the previous equation takes the form:

$$\frac{d{\tilde{\beta }}^i}{d\tau }=\sin (x_c{k}_u{\chi }^k){\tilde{\beta }}^j$$

Meanwhile, since $x_c$ is arbitrary, we may define $x_c$ freely. We thus choose $x_c=1/(k_{u}c)$,. This choice is convenient because if we then also define ${\tilde{\chi }}^k={\chi }^k/c$, we thus have ${\tilde{\beta }}^j=\frac{d{\tilde{\chi }}^k}{d\tau }$, allowing us to get rid of all the remaining factors of $c$ in the differential equation, and thus avoid all the numerical issues that having a factor of $c\sim 10^8$ causes. Thus we have the fully nondimensionalized variant of our differential equation:

$$\frac{d{\tilde{\beta }}^i}{d\tau }=\sin ({\tilde{\chi }}^k){\tilde{\beta }}^j$$

In our specific case, we have ${\tilde{\chi }}^k=\tilde{Z}$ where $\tilde{X}=\tilde{x}/c$ and $\tilde{Z}=\tilde{z}/c$, ${\tilde{\beta }}^i=({\tilde{\beta }}^x,{\tilde{\beta }}^z{)}^T$, and ${\beta }^j=(-{\tilde{\beta }}^z,{\tilde{\beta }}^x{)}^T$. Thus, we may write the above in explicit vector form as:

$$\frac{d}{d\tau }\left(\begin{array}{c}{\tilde{\beta }}^x\\ {\beta }^z\end{array}\right)=\sin (\tilde{Z})\left(\begin{array}{c}-{\tilde{\beta }}^z\\ {\beta }^x\end{array}\right)$$

Or equivalently:

$$\frac{d\tilde{\mathbf{\beta }}}{d\tau }=\mathbf{A}\tilde{\mathbf{\beta }},\tilde{\mathbf{\beta }}=\left(\begin{array}{c}{\tilde{\beta }}^x\\ {\tilde{\beta }}^z\end{array}\right),\mathbf{A}=\left(\begin{array}{cc}0& -\sin \tilde{Z}\\ \sin \tilde{Z}& 0\end{array}\right)$$

Note: assuming that the electrons move reasonably close to the beam axis and that the magnitude of their velocity stays reasonably constant (as Jackson assumes), one needs only to replace $\mathbf{A}$ with $\frac{1}{\gamma }\mathbf{A}$ to obtain a fully-relativistic result that applies in the strongly-relativistic regime (where $\gamma$ is, as usual, the Lorentz factor).

Note that to convert the above back to standard SI units (that is, to recover $x^j$ and ${\beta }^j$ from the numerically-computed ${\tilde{\chi }}^j$ and ${\tilde{\beta }}^j$), we make use of the conversions:

$$\begin{array}{rl}{x}^j& =x_c{\chi }^j=\frac{c}{k_{u}c}{\tilde{\chi }}^j=\frac{1}{k_u}{\tilde{\chi }}^j\\ {\beta }^j& =\frac{1}{c}\frac{d{x}^j}{dt}=\frac{1}{c}\frac{x_c}{t_c}\frac{d{\chi }^j}{d\tau }=\frac{x_c}{t_c}{\tilde{\beta }}^j=\frac{{\omega }_0}{k_{u}c}{\tilde{\beta }}^j\end{array}$$

Where the former expression is derived from ${\tilde{\chi }}^j={\chi }^j/c$, and the latter expression is derived from combining ${\dot{x}}^j=\frac{x_c}{t_c}\frac{d{\chi }^j}{d\tau }$, ${\beta }^j={\dot{x}}^j/c$, and ${\tilde{\beta }}^j=\frac{1}{c}\frac{d{\chi }^j}{d\tau }$. Equivalently, in vector notation, we have:

$$\begin{array}{rl}\mathbf{x}& =\frac{1}{k_u}\tilde{\mathbf{\chi }}\\ \mathbf{\beta }& =\frac{{\omega }_0}{k_{u}c}\tilde{\mathbf{\beta }}\end{array}$$

One must be careful to recognize that ${\omega }_0$ is actually negative to due to the fact that the electron has negative charge, so we must specify the magnitude of ${\omega }_0$ when discussing its value. For an undulator spatial period of $L_u=4.5\mathrm{cm}$ and a magnetic field strength of $B_0\approx 20\mathrm{mT}$ we have $|{\omega }_0|\approx 3.5\mathrm{GHz}$ and $k_u\approx 1.4\mathrm{c}{\mathrm{m}}^{-1}$. Thus one finds that:

$$\begin{array}{rl}{x}^j& \approx (0.716\mathrm{cm})\cdot {\chi }^j\\ {\beta }^j& \approx 0.084{\tilde{\beta }}^j\end{array}$$

Note that since ${\beta }^j$ is dimensionless, ${\tilde{\beta }}^j$ must also be for dimensional consistency; in addition, it is not difficult to see via dimensional analysis that ${\chi }^j$ has units of $\mathrm{c}{\mathrm{m}}^{-1}$ (inverse length). We will now describe several methods for solving the system using different software.

Approximate solution with leapfrog integration

We’ll start by writing a custom solver to numerically integrate the differential equation. It turns out that this is actually easier than using a ODE solving library, for reasons that will become apparent later. To do so, we use a method that is commonly-used in plasma physics, known as leapfrog integration. Leapfrog integration starts by approximating the continuous solution ${\tilde{\chi }}^j(\tau )$ as a series of points ${\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_n$ separated by finite intervals. As the number of intervals grows larger, the approximation grows more and more accurate.

The leapfrog integrator is a type of symplectic integrator, which roughly means that it ensures that physical systems that must follow conservation laws have generally-stable solutions that (approximately) preserve conserved quantities. It is an iterative algorithm that computes the position ${\mathbf{X}}_i$ and velocity ${\mathbf{V}}_i$ of a particle (or in our case, ${\tilde{\chi }}^j$ and ${\tilde{\beta }}^j$) for the next step $(i+1)$ from the current step $i$ as follows:

$$\begin{array}{rl}{\mathbf{X}}_{i+1}& ={\mathbf{X}}_i+h{\mathbf{V}}_i+\frac{h^2}{2}{\mathbf{a}}_i\\ {\mathbf{V}}_{i+1}& ={\mathbf{V}}_i+\frac{h}{2}({\mathbf{a}}_i+{\mathbf{a}}_{i+1})\end{array}$$

Where ${\mathbf{a}}_i=\mathbf{a}({\mathbf{X}}_i,{\mathbf{V}}_i)$ and $\mathbf{a}$ is defined (in our case) by:

$$\mathbf{a}({\mathbf{X}}_i,{\mathbf{V}}_i)=\sin ({\mathbf{X}}_i){\mathbf{V}}_i$$

This can be implemented by the following (simplified) Python code:

def leapfrog(fun, x0, v0, t_span, step_size=1E-4):
  t0, t_end = t_span
  N = int((t_end - t0)/step_size) # num steps
  t = np.linspace(t0, t_end, N)
  x = np.zeros((N, len(x0)))
  v = np.zeros((N, len(v0)))
  x[0] = x0
  v[0] = v0
  h = step_size
  for i in range(N-1):
    # calculate a[i]
    a_current_step = fun(t[i], x[i], v[i])
    x[i+1] = x[i] + v[i]*h + 1/2 * a_current_step*h**2
    # calculate a[i+1]
    if fun_args is not None:
      a_next_step = fun(t[i], x[i+1], v[i], **fun_args)
    else:
      a_next_step = fun(t[i], x[i+1], v[i])
    v[i+1] = v[i] + (a_current_step + a_next_step)/2 * h
  return {
      "t": t,
      "x": x,
      "v": v
  }

Note: The code for the numerical solver can be found in the Jupyter notebook Simulations/UndulatorDesignSimulations.ipynb in the laser research repository. We have found that it has yielded the best results of any numerical method for computing beam trajectories of any of the methods discussed.

Validation of numerical solution with SciPy

To check that the solution of the leapfrog integrator is correct, an alternative (though very-similar) formulation of the ODEs can be solved via solve_ivp’s ODE solver. To start, we have:

$$\begin{array}{rl}m\ddot{x}& =-q\dot{z}{B}_y\\ m\ddot{z}& =q\dot{x}{B}_y\end{array}$$

Or alternatively, after dividing both sides by $c$:

$$\begin{array}{rl}{\dot{\beta }}_x& =-\frac{q}{m}{\beta }_z{B}_y\\ {\dot{\beta }}_z& =\frac{q}{m}{\beta }_x{B}_y\end{array}$$

Where $B_y=B_0\sin (k_{u}z)$. Now after defining $Z=z/c$, $X=x/c$ as well as ${\omega }_0=q{B}_0/m$ we have:

$$\begin{array}{rl}{\dot{\beta }}_x& =-{\omega }_0{\beta }_z\sin (\omega Z)\\ {\dot{\beta }}_z& ={\omega }_0{\beta }_x\sin (\omega Z)\end{array}$$

Where $\omega =k_{u}c$ is a precomputed constant (of the order $10^{-1}-10^{-2}$ in units of THz), and ${\omega }_0$ is of the same order of magnitude (in units of THz). We thus solve for the scaled positions $X,Z$ and the scaled velocities ${\beta }_x,{\beta }_z$, which has the advantage that the latter are direct derivatives of the former, without needing to multiply by a factor of $c$. Note also that $e/m$, the electron charge-to-mass ratio, has a well-established CODATA value that can be used when computing ${\omega }_0$. However, since the oscillation frequencies are still so high (in the THz range), we perform a transformation of variables $\xi ={\omega }_{0}t$ and $({\Lambda }_x,{\Lambda }_z)=(\omega X,\omega Z)=\frac{\omega }{c}(x,z)$. Thus:

$$\begin{array}{rl}{\beta }_x& ={\omega }_0\frac{dX}{d\xi }=\frac{{\omega }_0}{\omega }\frac{d{\Lambda }_x}{d\xi }\\ {\dot{\beta }}_x& ={\omega }_0\frac{d{\beta }_x}{d\xi }=\frac{{\omega }_0^2}{\omega }\frac{d^2{\Lambda }_x}{d{\xi }^2}\\ {\beta }_z& ={\omega }_0\frac{dZ}{d\xi }=\frac{{\omega }_0}{\omega }\frac{d{\Lambda }_z}{d\xi },\\ {\dot{\beta }}_z& ={\omega }_0\frac{d{\beta }_z}{d\xi }=\frac{{\omega }_0^2}{\omega }\frac{d^2{\Lambda }_z}{d{\xi }^2}\end{array}$$

Therefore, after substitution, we have transformed the system into two second-order quasilinear differential equations for ${\Lambda }_x(\xi )$ and ${\Lambda }_z(\xi )$:

$$\begin{array}{rl}\frac{d^2{\Lambda }_x}{d{\xi }^2}& =-\frac{d{\Lambda }_z}{d\xi }\sin ({\Lambda }_z)\\ \frac{d^2{\Lambda }_z}{d{\xi }^2}& =\frac{d{\Lambda }_x}{d\xi }\sin ({\Lambda }_z)\end{array}$$

Or equivalently, in matrix-vector form:

$$\frac{d^2}{d{\xi }^2}\left(\begin{array}{c}{\Lambda }_x\\ {\Lambda }_z\end{array}\right)=\left(\begin{array}{cc}0& -\sin {\Lambda }_z\\ \sin {\Lambda }_z& 0\end{array}\right)\frac{d}{d\xi }\left(\begin{array}{c}{\Lambda }_x\\ {\Lambda }_z\end{array}\right)$$

This completely removes dependence on the frequencies, which are now entirely contained within the definitions of the new variables ${\Lambda }_x,{\Lambda }_z,\xi$. Note that in our initial conditions, we must use the following transforms:

$$\begin{array}{rl}\frac{d{\Lambda }_x}{d\xi }{|}_{{\xi }_{initial}}& =\left(\frac{\omega }{{\omega }_0}\right){\beta }_x{|}_{t_{initial}}\\ \frac{d{\Lambda }_z}{d\xi }{|}_{{\xi }_{initial}}& =\left(\frac{\omega }{{\omega }_0}\right){\beta }_z{|}_{t_{initial}}\end{array}$$

Note: The code for the SciPy solver can be found in the Jupyter notebook Simulations/UndulatorDesignSimulations.ipynb in the laser research repository. However, the results are still work-in-progress and yield inconclusive results at the present time.

Approximate solution with OSCARS

We can also use the dedicated OSCARS code to numerically compute the results. A table of the results for different beam energies are shown below; all tests are done with $N=10$ and $L_u=\text{2 cm}$.

Magnetic field	Beam energy	Current*	Orientation	Success?
0.1 T	1 MeV	1.0 A	$\theta =0$	Yes
0.05	520 keV	1.0 A	$\theta =0$	Yes
0.03	515 keV	1.0 A	$\theta =0$	Yes
0.03	510 keV	1.0 A	$\theta =0$	No
0.01	510 keV	1.0 A	$\theta =0$	No
0.01	512 keV	1.0 A	$\theta =45^{\circ }$	No
0.01	512 keV	1.0 A	$\theta =15^{\circ }$	No
0.01	512 keV	1.0 A	$\theta =0^{\circ }$	Yes
0.005	512 keV	1.0 A	$\theta =0$	Yes
0.005	510 keV	1.0 A	$\theta =0$	No
0.005	510 keV	1.0 A	$\theta =45^{\circ }$	No

*: Here $\vec{d}=(\sin \theta ,0,\cos \theta {)}^T$ (in $x,y,z$ coordinates) is the initial direction of the particles, with $\theta$ being the angle to the $z$-axis (optical axis).

Note: The code for the numerical solver can be found in Oscars/OscarsBasicSimulation.ipynb in the laser research repository.

The results indicate a cutoff beam energy of 512 keV would be required, which is indeed monstrously high. This, however, may be due to the fact that OSCARS was designed for simulating free-electron lasers with electron beam energies in the GeV range, not free-electron masers with beam energies in the keV range, which is a difference of 6 orders of magnitude.

Approximate solution with SciPy (archived)

Note: This section is an older section, and is preserved here only for archival purposes.

We’ll now cover an alternative method for numerically calculating the beam trajectories using the solve_ivp solver from the SciPy library.

Typically, ODE solvers expect coupled ODEs to be written in vector ODE form $\dot{x}=f(x,z),\dot{z}=f(x,z)$. This means we must write the system in its phase-space form. Let us define $\mathbf{U}=(x,z,\dot{x},\dot{z})=(U_0,U_1,U_2,U_3)$ such that we have:

$$\dot{\mathbf{U}}=f(\mathbf{U}),f(\mathbf{U})=\left(\begin{array}{c}{v}_x\\ v_z\\ -{\omega }_0\sin (k_{u}z)\dot{z}\\ {\omega }_0\sin (k_{u}z)\dot{x}\end{array}\right)=\left(\begin{array}{c}{U}_2\\ U_3\\ -{\omega }_0\sin (k_u{U}_1)U_3\\ {\omega }_0\sin (k_u{U}_1)U_2\end{array}\right)$$

This can then be numerically implemented with the following code:

# right-hand side of differential equation
def f_vector(t, y_vec, omega0=1E9, ku=2*np.pi/5):
	# decompose the phase-space vector y into its components
	x, z, vx, vz = y_vec
	# prefactor
	fac = omega0 * sin(ku*z)
	return [vx, vz, -fac*vz, fac*vx]

Then, supplying f_vector as the first argument to solve_ivp yields the solution. See the SciPy solve_ivp docs for further information.

Fully relativistic solution with SciPy

Our previous analysis assumed that we could set $\gamma \approx 1$, which holds in the non-relativistic and mildly-relativistic regime. However, if we want our solution to hold in the relativistic regime, we cannot make this assumption, and must solve with $\gamma$ explicitly in the differential equations.

Unfortunately, since the Lorentz factor itself also depends on $\dot{x}$ and $\dot{z}$, it cannot simply be factored out of the derivative; instead, we must use the expanded form, given by:

$$\dot{\gamma }\frac{dz}{dt}-\gamma \frac{d^{2}z}{d{t}^2}=\frac{q{B}_0}{m}\dot{x}\sin k_{0}x$$

Which can be rearranged to yield:

$$\frac{d^{2}z}{d{t}^2}=\frac{1}{\gamma }\left[\dot{\gamma }\frac{dz}{dt}-\frac{q{B}_0}{m}\dot{x}\sin k_{0}x\right]$$

However, since we now have second-order time derivatives, we must perform a reduction of order. We can do this by converting the two differential equations into four first-order differential equations:

$$\begin{array}{rl}\frac{dx}{dt}& =v_0,\frac{dz}{dt}=v_z,\\ \frac{d{v}_x}{dt}& =0,\frac{d{v}_z}{dt}=\frac{1}{\gamma }\left[\dot{\gamma }{v}_z-\frac{q{B}_0}{m}{v}_x\sin k_{0}x\right]\end{array}$$

Note that $\dot{\gamma }$, the time derivative of the Lorentz factor, is given by:

$$\begin{array}{rl}\frac{d\gamma }{dt}& =\frac{d}{dt}\left(\frac{1}{\sqrt{1-v^2/c^2}}\right)\\ & =\frac{v}{c^2(1-v^2/c^2{)}^{3/2}}\frac{dv}{dt}\\ & =\frac{v}{c^2(1-v^2/c^2{)}^{3/2}}|\mathbf{a}|,v=|\mathbf{v}|=\sqrt{{\dot{x}}^2+{\dot{z}}^2},|\mathbf{a}|=\sqrt{{\dot{v}}_x^2+{\dot{v}}_z^2}\end{array}$$

Where $v=\sqrt{{\dot{x}}^2+{\dot{z}}^2}$ and $|\mathbf{a}|=\sqrt{{\ddot{x}}^2+{\ddot{z}}^2}=|{\dot{v}}_z|$ (since we have $\ddot{x}={\dot{v}}_x=0$ from our system of differential equations). Thus, substituting, we have:

$$\frac{d\gamma }{dt}=\frac{v}{c^2(1-v^2/c^2{)}^{3/2}}\frac{d{v}_z}{dt}$$

We note that since $\gamma =\frac{1}{(1-v^2/c^2{)}^{1/2}}$ we can also write our result as follows:

$$\frac{d\gamma }{dt}=\frac{v}{c^2}{\gamma }^3\frac{d{v}_z}{dt}$$

However, our system is still not in the correct form, since we have a term involving time derivatives of $v_z$ on both the left- and right-hand sides for the 4th ODE in the system (the one for $\frac{d{v}_z}{dt}$). To see why, we can expand the 4th ODE using our definition of $\dot{\gamma }$, giving us:

$$\frac{d{v}_z}{dt}=\frac{v}{c^2}{v}_z{\gamma }^2\frac{d{v}_z}{dt}-\frac{1}{\gamma }\frac{q{B}_0}{m}{v}_x\sin k_{0}x$$

By slight rearrangement, we can resolve our issue by moving the $\frac{d{v}_z}{dt}$ term to the left-hand side, giving us (after multiplying again by $\gamma$):

$$\begin{array}{c}\frac{d{v}_z}{dt}\left[\frac{v}{c^2}{v}_z{\gamma }^2-1\right]=\frac{1}{\gamma }\frac{q{B}_0}{m}{v}_x\sin k_{0}x\\ \gamma \frac{d{v}_z}{dt}\left[\frac{v}{c^2}{v}_z{\gamma }^2-1\right]=\frac{q{B}_0}{m}{v}_x\sin k_{0}x\\ \frac{d{v}_z}{dt}=\frac{v_x}{v{\gamma }^3{v}_z/c^2-\gamma }\frac{q{B}_0}{m}\sin k_{0}x\end{array}$$

Thus, our system of ODEs becomes:

$$\begin{array}{rl}\frac{dx}{dt}& =v_0,\frac{dz}{dt}=v_z,\\ \frac{d{v}_x}{dt}& =0,\frac{d{v}_z}{dt}=\frac{v_x}{v{\gamma }^3{v}_z/c^2-\gamma }\frac{q{B}_0}{m}\sin k_{0}x\end{array}$$

Which can be written in the matrix form ${\mathbf{y}}^{\prime }=f(\mathbf{y},t)$ as follows:

$$\mathbf{y}=\left(\begin{array}{c}x\\ z\\ v_x\\ v_z\end{array}\right),f(\mathbf{y},t)=\left(\begin{array}{c}{v}_0\\ v_z\\ 0\\ \frac{v_x}{v{\gamma }^3{v}_z/c^2-\gamma }\frac{q{B}_0}{m}\sin k_{0}x\end{array}\right)$$

We can translate this system to Python code as follows:

import numpy as np
 
c = 3E8 # m/s
electron_charge = -1.602E-19 # coulomb
B0 = 1.2 # tesla
electron_mass = 9.109E-31 # kilograms
 
# calculate the speed (scalar) from
# a velocity vector
def speed(vx, vz):
	return np.linalg.norm([vx, vz])
 
# calculate the magnitude of the acceleration
# from an acceleration vector
def accel(ax, az):
	return np.linalg.norm([ax, az])
 
# calculate the lorentz factor from the velocity
def lorentz(v, c=c):
	return 1/np.sqrt(1 - v**2 / c**2)
 
# right-hand side of differential equation
def f_vector(t, y_vec, 
		c=c, v0=0.6*c, q=electron_charge, B0=B0, m=electron_mass):
	# decompose the phase-space vector y into its components
	x, z, vx, vz = y_vec
	v = speed(vx, vz)
	# calculate the Lorentz factor
	gamma = lorentz(v, c)
	# now supply the derivatives
	dx_dt = v0
	dz_dt = vz
	dvx_dt = 0
	dvz_dt = vx/(v * vz * gamma**3/c**2 - gamma) * q*B0/m * sin(k0 * x)
	return [dx_dt, dz_dt, dvx_dt, dvz_dt]

Note: A real electron beam is not a perfect curve and has a nonzero beam cross-section, which is caused by electrons having different energies (and thus different velocities) and some electrons travelling off-centerline. We can simulate this by numerically integrating with different starting positions $(x_0,z_0)$ and with different initial velocities $v_0$. This gives us multiple trajectories that we can overlay together to get a better picture of what a realistic beam would look like. An alternative approach is detailed when we get to electron gun modelling (which will come later).

Modelling of permanent magnets for undulators

The magnetic field can be found using the magnetostatic equations of motion:

$$\{\begin{array}{l}\nabla \cdot (\mu (-\nabla {\psi }_m)+\mathbf{M})=0\\ \mathbf{B}=\mu \mathbf{H}=-\mu (-\nabla {\psi }_m+\mathbf{M})\end{array}$$

Note: Here $\mu$ is the permeability of the magnets, which can be found on tables of material properties and is also usually provided by the manufacturer that made the magnets. For our purposes we can assume $\mu =\text{const.}$

Analytical and numerical solution

For the undulator, we have so far used a simple analytical expression for the magnetic field $B_y$, which comes from Jackson. It is given by:

$$B_y=B_0\sin k_{u}z$$

The reason this expression for the field works reasonably well is that we expect the dominant mode in the Fourier expansion of the field to be the fundamental mode. To describe a realistic beam, however, will involve a Fourier series and will depend on $x$ and $z$, that is:

$$B_y=B_0\sum _{k_u^{(n)}}{c}_n\sin (k_u^{(n)}z)e^{-c_n{x}^2}$$

Where $k_u^{(n)}$ is the wavevector associated with the $n$-th mode, and where we assume a Gaussian cross-section of the field. It is possible to find the coefficients $c_n$ by either fitting to an experimentally-measured beam or approximating it with some theoretical assumptions, but here we will not do so to avoid needless complexity and verbosity. To be able to find more accurate solutions for fields within the undulator, please see Numerical simulations for magnetostatics.

Electron gun modelling

Understanding electron beams also involves being able to model electron guns, which are what produce electron beams in the first place. Electron guns come in a wide variety of designs, but here we will focus on electrostatic electron guns, and discuss how they can be modelled and simulated. In an electrostatic electron gun, a potential is maintained between a cathode (negatively-charged plate) and anode (positively-charged plate). This potential difference is what “pushes” free electrons between the plates, which can then be collimated to form a beam (see Vacuum tube Lorentz factor calculation for more details). We will explain how to model an (electrostatic) electron gun below.

Modelling electron trajectories through electron gun

The key aim in electron gun modelling is beam tracing: tracing the paths of electrons as they exit the cathode surface and are propelled and collimated into a beam. This is significant because the whole point of an electron gun is to create a collimated beam of free electrons. Since power is proportional to both current and voltage, the more free electrons are in the beam, the more efficient the electron gun will be.

In an electron gun, several key physical mechanisms are at play, which must all be accounted for in a realistic model of an electron gun. First, electrons “boil off” the cathode due to primarily thermionic emission, at which point the electrons have random velocities and are moving in random directions². To be able to focus the electrons into a beam, a potential is applied. The applied potential difference $\Delta V_{ac}$ between the cathode and anode is a constant value $V_0$, that is:

$$\Delta V_{ac}=V_0$$

Note: Be careful to note that while the potential difference is constant, the electric potential itself isn’t! Indeed, if it was, we would have an equipotential surface, meaning that there would be no electrostatic acceleration at all!

Considering the problem in 1 dimension, we may solve for the potential via the boundary conditions³ $V_{ac}(0)=0$ and $V_{ac}(d)=V_0$, where $d$ is the distance between the cathode and anode. Thus, by solving Laplace’s equation for electrostatics (that is, ${\nabla }^{2}V=0$) in one dimension (in which case it becomes $V^{\prime \prime }(z)=0$ where $z$ is the beam’s axis of propagation), one finds that:

$$V_{ac}(x)=\alpha z$$

Where $\alpha \equiv V_0/d$ is the linear potential density (potential per unit length). The applied electric field is thus given by:

$${\mathbf{E}}_{ac}=-\nabla V=-\alpha \hat{\mathbf{z}}$$

However, the applied electric field is not the only field that affects the motion of the electrons. First, we must consider the effects of the electron self-potential, caused by the electrostatic repulsion of the electrons within the beam. This causes the beam to diverge, which requires collimating magnetic lenses to correct for the divergence and collimate the beam. Second, we must model the collimating magnets themselves, in which case we must consider the magnetic term in the Lorentz force law as well as the electric term.

Given the complexity of the dynamics involved, it is easiest to model the magnetic field the electron beam divergence caused by the self-potential as two separate problems. Of course, the underlying assumption is that they are indeed able to be treated as separate systems; mathematically, this assumption is valid if the self-potential-caused beam divergence is minimal over the length scales where the magnetic field is strong. We will assume this for the rest of our discussion.

We will now consider the behavior of the electrons after they boil off the cathode (we’ll discuss what happens prior in a section later down), excluding the electrostatic acceleration from the anode-cathode-potential. The velocity distribution of the electrons are thus encoded in the initial conditions, which is decoupled from their equations of motion. To compute the electron trajectories, we will once again need to use the Lorentz force law, whose most general form is given by:

$$\frac{d}{dt}(\gamma m\mathbf{v})=q(\mathbf{E}+\mathbf{v}\times \mathbf{B})$$

In the non-relativistic regime, we have:

$$\frac{d}{dt}(m\mathbf{v})=q(\mathbf{E}+\mathbf{v}\times \mathbf{B})$$

This is a reasonable approximation if the transverse acceleration of the electrons is generally non-relativistic or mildly-relativistic, which is the case for most electron guns used in free-electron masers. We now consider the electron self-potential $V_e$ caused by electron-electron repulsion in the beam. We may state it as follows:

$$V_e=\frac{1}{2}\sum _n\sum _m\frac{-e}{4\pi {\epsilon }_0}\frac{1}{|r_m-r_n|}$$

By tracing the paths of the electrons chosen from different positions across the plate, one may find the general shape of the electron beam. As for their velocities, since the magnitude of the electron velocities must be equal or greater than the work function, the kinetic energy of the electrons is thus (non-relativistically):

$$K=\frac{1}{2}m{v}_x^2\ge W$$

Assuming the minimal kinetic energy $K=W$ to overcome the material’s work function and thus allow the electrons within the material to escape, we must have:

$$v_0=\sqrt{2W/m}$$

The problem can then be solved numerically with suitable assumptions for the initial conditions (for instance, $v_{0x}\approx v_0$ and $v_{0y}=v_{0z}\approx 0$). There are different algorithms for this purpose, such as the Barnes-Hut algorithm in electrostatics and the multigrid method for Poisson’s equation. A promising alternative is the particle-in-cell method, commonly used in plasma physics. However, a discussion of all of these algorithms would be far, far too lengthy here. Instead, we will simply show the analytical solution to the problem.

Analytical solution for electron trajectories with self-potential

In the absence of a magnetic field, the only interaction that shapes the electron beam is the inter-electron repulsion, which leads to beam divergence. The analytical solution of the problem, which is valid for both the non-relativistic and relativistic regimes, is given in implicit form by⁴:

$$z(\chi )=\frac{R_{0}F(\chi )}{\sqrt{2K}}$$

Where $z$ is the axis of propagation for the beam, $R_0$ is the beam waist (the initial width of the beam at $z=0$), $\chi =x/R_0$ is the ratio of the beam radius $x$ to the beam waist, $K$ is a constant, and $K$ is a dimensionless variable known as the generalized beam perveance, given by⁵:

$$K=\frac{e{I}_0}{2\pi {\epsilon }_0{m}_e(v_z\gamma {)}^3}$$

Where $e$ is the magnitude of the elementary charge, ${\epsilon }_0$ is the permittivity of free space, $m_e$ is the electron rest mass, $\gamma$ is the Lorentz factor, $I_0$ is the beam current, and $v_z$ is the electron velocity in the $z$ direction. Note that the Lorentz factor $\gamma$ can be computed from the beam energy (see Vacuum tube Lorentz factor calculation for more details), and if we assume that the divergence is fairly small ($|v|\approx v_z$) then:

$$v_z=\beta c=c\left(1-\frac{1}{{\gamma }^2}\right)$$

Finally, $F(\chi )$ is defined by the following logarithmic integral:

$$F(\chi )={\int }_1^{\chi }\frac{dy}{\sqrt{\ln y}}$$

Note: it is also possible to write the solution in the more familiar (implicit) form $z=z(x)$ with $z(x)=\frac{R_{0}F(x/R_0)}{\sqrt{2K}}$. However, when the beam divergence is very small (on the order of micrometers), this becomes difficult to evaluate numerically, so it is preferred to use the dimensionless variable $\chi$ rather than $x$.

Note that as long as the distance between the cathode and anode is short (i.e. $z$ is small), the solution can be approximated via:

$$z(\chi )\approx R_0\sqrt{\frac{2\ln \left(\chi \right)}{K}}$$

In which case it is possible to write out the solution explicitly as a function of $z$, as follows:

$$\begin{array}{c}{z}^2=R_0^2\frac{2\ln \chi }{K}\\ \\ \frac{K}{2R_0^2}{z}^2=\ln \chi \\ \\ \chi (z)=e^{K{z}^2/2R_{0^2}}\end{array}$$

This comes from the leading-order term in the series expansion for $F(\chi )$, which is a Puisseux series in the form:

$${\int }_1^{\chi }\frac{dy}{\sqrt{\ln y}}=2\ln {\left(\chi \right)}^{\frac{1}{2}}+\frac{2}{3}\ln {\left(\chi \right)}^{\frac{3}{2}}+\frac{1}{5}\ln {\left(\chi \right)}^{\frac{5}{2}}+\dots$$

One may derive this with the substitution $u=\ln y$⁶, $du=\frac{1}{y}dy=e^{-u}dy$. This then gives us:

$${\int }_1^{\chi }\frac{dy}{\sqrt{\ln y}}={\int }_0^b{u}^{-1/2}{e}^{u}du$$

Where $b\equiv \ln \chi$. This integral can be evaluated by means of the following identity:

$$\int x^n{e}^{ax}dx=\frac{(-1{)}^n}{a^{n+1}}\Gamma (1+n,-ax)$$

Where $\Gamma (s,x)$ is the incomplete upper gamma function, given by:

$$\Gamma (s,x)={\int }_x^{\infty }{t}^{s-1}{e}^{-t}dt$$

In our case, we have $n=-1/2$ and $a=1$, therefore:

$${\int }_0^b{u}^{-1/2}{e}^{u}du=-i\Gamma \left(\frac{1}{2},-x\right){|}_0^b$$

Note that by the properties of the incomplete gamma function, $\Gamma \left(\frac{1}{2},x\right)=\sqrt{\pi }\mathrm{erfc}(\sqrt{x})$ where $\mathrm{erf}(x)$ is the error function and $\mathrm{erfc}(x)=1-\mathrm{erf}(x)$ is the complementary error function Thus the above reduces to:

$$-i\Gamma \left(\frac{1}{2},-x\right){|}_0^b=-i\sqrt{\pi }\mathrm{erfc}(\sqrt{-x}){|}_0^b=-i\sqrt{\pi }(\mathrm{erfc}(\sqrt{-b})-1)=i\sqrt{\pi }\mathrm{erf}(\sqrt{-b})$$

The Maclaurin series of the error function is given by:

$$\begin{array}{rl}\mathrm{erf}(z)& =\frac{2}{\sqrt{\pi }}\sum _{n=0}^{\infty }\frac{(-1{)}^n{z}^{2n+1}}{n!(2n+1)}\\ & =\frac{2}{\sqrt{\pi }}\left(z-\frac{1}{3}{z}^3+\frac{1}{10}{z}^5-\frac{1}{42}{z}^7+\frac{1}{216}{z}^9+\dots \right)\end{array}$$

Thus we have:

$$\mathrm{erf}(\sqrt{-b})=\frac{2i}{\sqrt{\pi }}\left(b^{1/2}+\frac{1}{3}{b}^{3/2}+\frac{1}{10}{b}^{5/2}+\frac{1}{42}{b}^{7/2}+\frac{1}{216}{b}^{9/2}+\dots \right)$$

And thus:

$$i\sqrt{\pi }\mathrm{erf}(\sqrt{-b})=2\left(b^{1/2}+\frac{1}{3}{b}^{3/2}+\frac{1}{10}{b}^{5/2}+\frac{1}{42}{b}^{7/2}+\frac{1}{216}{b}^{9/2}+\dots \right)$$

Thus, the solution of the integral can be written in series form as:

$${\int }_0^{\ln (\chi )}{u}^{-1/2}{e}^{u}du=2\ln (\chi {)}^{1/2}+\frac{2}{3}\ln (\chi {)}^{3/2}+\frac{2}{10}\ln (\chi {)}^{5/2}+\frac{2}{42}\ln (\chi {)}^{7/2}+\frac{2}{216}\ln (\chi {)}^{9/2}+\dots$$

And thus, to leading-order, we obtain the following approximation, which we may then substitute into $z(\chi )$:

$${\int }_1^{\chi }\frac{dy}{\sqrt{\ln y}}\approx 2\sqrt{\ln \left(\chi \right)}$$

Magnetic lens modelling

Now, as the electrons pass through the magnetic lens comprised of 2 sextupole magnets, we also need to factor in the effect of the magnetic field. Here it is easiest to work with the (magnetic) vector potential $\mathbf{A}$ rather than the magnetic field $\mathbf{B}$; they are related by $\mathbf{B}=\nabla \times \mathbf{A}$. In our case, the vector potential satisfies the differential equation:

$${\nabla }^2\mathbf{A}=0$$

This is simply Laplace’s equation for magnetostatics, and it can be solved via the multipole expansion. In polar coordinates one finds that the vector potential takes the form:

$$\begin{array}{r}{A}_x={\left(\frac{A}{C\cos 3\theta }\right)}^{1/3}\cos \theta \\ A_y={\left(\frac{A}{C\cos 3\theta }\right)}^{1/3}\sin \theta \end{array}$$

Where $A,C$ are some constants determined by the boundary conditions. This polar form is useful because the sextupole magnet is rotationally-symmetric. We will ignore the $A_z$ component here since it turns out it is unimportant.⁷

Note: See a visualization of the magnetic vector potential at https://www.math3d.org/DSWoli41Gv. Note that the magnetic field $\mathbf{B}$ is always perpendicular to the vector potential $\mathbf{A}$.

The equations of motion are easiest to find from the electromagnetic Hamiltonian, which may be written in terms of the position $\mathbf{R}$ and canonical momentum $\mathbf{P}$ as follows:

$$\hat{H}(\mathbf{R},\mathbf{P})=\frac{(\mathbf{P}-q\mathbf{A})}{2m}+qV$$

This gives you Hamilton’s equations, a system of two first-order differential equations, which are respectively given by:

$$\begin{array}{rl}\frac{d\mathbf{R}}{dt}& =\frac{\mathbf{P}-q\mathbf{A}}{m}=\mathbf{v}\\ \frac{d\mathbf{P}}{dt}& =q\left(-\nabla V-\frac{\partial \mathbf{A}}{\partial t}+\mathbf{v}\times (\nabla \times \mathbf{A})\right)\end{array}$$

Excluding the electric potential $V$ (since we are considering only the effects of the magnetic lens) and noting that $\mathbf{A}$ is time-independent (since we are considering only magnetostatics), we have:

$$\begin{array}{rl}\frac{d\mathbf{R}}{dt}& =\mathbf{v}\\ \frac{d\mathbf{P}}{dt}& =q\mathbf{v}\times (\nabla \times \mathbf{A})\end{array}$$

Note that since $\mathbf{B}=\nabla \times \mathbf{A}$, the first-order equations obtained from the Hamiltonian approach end up exactly equivalent to the Lorentz force law. Now expanding with the formulas for the curl in polar coordinates⁸, we obtain a system of differential equations for $(\mathbf{x},\mathbf{P})$, given by:

$$\begin{array}{rl}\dot{x}& =P_x/m\\ \dot{y}& =P_y/m\\ {\dot{P}}_x& =-\frac{q{P}_y}{m}\frac{f(\theta )}{r}\\ {\dot{P}}_y& =\frac{q{P}_x}{m}\frac{f(\theta )}{r}\end{array}$$

Where $r=\sqrt{x^2+y^2}$, $\theta ={\tan }^{-1}(y/x)$, and:

$$f(\theta )={\left(\frac{A}{C\cos 3\theta }\right)}^{1/3}\tan \theta$$

It should be more than obvious that this system is impossible to solve analytically without a slew of approximations. However, it can be solved numerically using the same techniques as described previously for the undulator beam.

Note: See notebooks/magnetic-lens-calculations.ipynb for the derivation of the equations of motion in explicit form.

Thermocathode modelling

The last part, which we have left out so far, is how the electrons are ejected from the first place from the cathode prior to being collimated into a beam. This is a complex topic, so it is out of scope for this page, but please see Thermocathode design for electron gun for more details.

Software implementation of numerical simulations

For our simulations of the undulator, we use OSCARS (an open-source specialized electron beam/undulator modelling software from Brookhaven). OSCARS has an advanced beam-tracing algorithm that calculates the magnetic fields and beam trajectories inside an undulator. However, we also use our own custom-written codes. The simulations are all open-source and can be found here here.

References

• OSCARS reference examples
  • Magnetic field modelling (also includes information on how to load magnetic fields from external softwares)
  • Undulator electron trajectory simulations
  • Particle beam modelling (this also includes multi-particle beam modelling, which is very important for our use case)

Footnotes

1. We substituted in $z\approx v_{0}t$ into $B_y=B_0\sin k_{u}z$ and then substituted in $B_y$ to obtain this result. ↩

2. Note that electrons that “boil off” the cathode are initially randomly-spread out and have a wide distribution of velocities (and thus kinetic energies) that is roughly a Gaussian distribution (see Guide to modelling lasers and masers for more details), although more accurately it is a Fermi-Dirac distribution given that electrons are fermions. ↩

3. Technically speaking, since the cathode is negatively-charged and the anode positively-charged (with equal charges on both), we would have $V_{ac}(0)=-V_0/2$ and $V_{ac}(d)=V_0/2$. However, since the electric potential is arbitrary and only the potential difference can be measured, we simply add the constant $V_0/2$ to make the problem nicer to work with (without changing the physics). ↩

4. This result comes from the following COMSOL tutorial, which itself references the book Charged Particle Beams by Stanley Humphries. Also, the Wikipedia article on charged particle beams may be useful reading. ↩

5. Some other texts (including this document from CERN’s accelerator school, see eq. 24) write out the formula with $v_z=c^3{\beta }^3$, where $\beta \equiv |v|/c$. This form is reasonable if the beam is paraxial (i.e. does not diverge too much) such that $v\approx v_z$. ↩

6. As suggested by this Physics forum article. ↩

7. Discussed on pg. 23 of the lecture notes from the US particle accelerator school ↩

8. See this Wikipedia article, which lists the vector differential operators in Cartesian, cylindrical, and spherical coordinates. ↩