Guide to modelling 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. This guide is a comprehensive overview of how we can model electron beams mathematically and computationally.

Note: This is a general overview on modelling electron beams, with primary emphasis on analytic methods. For specific information about using numerical methods for electron beam modelling, please see Numerical modelling of electron beams.

Modelling of undulator electron beams

Here, we will use slightly different coordinate conventions compared to Free-electron maser physics, in an attempt to match the standard literature convention. 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 system of differential equations, one we will need to solve (either through analytical or numerical means) to determine the trajectory of the electron beam through undulator. It is pretty obvious that an exact analytical solution is near-impossible to obtain, and we will need to use numerical methods.

Why do we care about the electron beam trajectories? Well, because a well-behaved electron beam is essential to high-efficiency operation of the maser. If the electron beam loses collimation or diverges, the power output of the maser will drop, and energy would be wasted; indeed, in some cases, the maser may not be able to function at all. This is why it is so important to calculate (or simulate) electron beam trajectories, prior to the start of construction.

Fundamental parameters controlling electron beams

From the Lorentz force equations, we can clearly see that the trajectories of charged particles in a magnetic field are dependent on their velocities, their charge, and the magnetic field strength (and orientation). In our case, the beam is made of one type of particle (electrons) so the charge is fixed; however, the electron velocity can be varied by changing the accelerator voltage (see Vacuum tube Lorentz factor calculation) and the magnetic field strength can be modified by raising and lowering the magnets. These two will be the key parameters used to control the electron beam.

A heuristic derivation of the beam energy threshold

The electron beam’s behavior is dependent on many variables, but a particularly important one is the beam energy: the average kinetic energy of electrons in the beam. A beam energy that is too low will lead to a beam that is deflected straight back by the undulator magnets. A beam energy that is too high is very inefficient. Finding the right beam energy is essential, but very difficult to calculate by analytical methods.

A good way to try to estimate the beam energy necessary is to consider just one magnet, which we assume produces a constant magnetic field of $\mathbf{B}=B_0\hat{z}$. Note that this field is 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 avoid ending up in a circular orbit, and successfully “escape” to reach the next magnet in the undulator. Note that we can rearrange our previous equation in terms of the normalized velocity $\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_{mass}}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}$$

For a field strength of $B_0=0.1\mathrm{T}$, 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$. However, using weaker field strengths, the beam energy requirements are much lower; for instance, at a field strength of $5\mathrm{mT}$ and $L_u=10\mathrm{cm}$ the required velocity is less than a percent of the speed of light.

Undulator modelling

Now that we have a rough benchmark, we can attempt to analytically solve for the electron beam trajectories. This will only allow us to get an approximate result, for reasons that will become clear shortly; but it does provide us with a useful reference result.

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 very idealized conditions, 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 methods for undulator simulations

Information about numerical methods can be found in Numerical modelling of electron beams.

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[^5]. 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[^6] $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[^7]:

$$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[^9]:

$$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$[^8], $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.[^12]

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[^11], 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.

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. ↩