Numerical modelling of electron beams

In Guide to modelling electron beams, we discussed (primarily) analytical methods for modelling electron beams, a key component in our free-electron maser. However, in the vast majority of cases, it is impossible to find analytical solutions and we must resort to numerical methods. Here, we give a description of 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).

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 our code samples, we will use Python, but any other programming language can be used to implement the numerical methods described here.

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 loss of significance. To avoid this we will perform nondimensionalization. First, note that 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).

The Hamiltonian approach

Lastly, we’ll discuss one approach that is more advanced, but especially powerful when used alongside symplectic integrators. This is the Hamiltonian approach. We start with the classical Hamiltonian of an electron in a magnetic field, which is given by:

$$H(\mathbf{r},\mathbf{p},t)=\frac{1}{2m}(\mathbf{p}+e\mathbf{A}(\mathbf{r}){)}^2$$

Where $\mathbf{A}$ is the magnetic vector potential, and $\mathbf{B}=\nabla \times \mathbf{A}$. The equations of motion are then naturally expressed in terms of the Hamiltonian, as follows:

$$\frac{d\mathbf{r}}{dt}=\frac{\partial H}{\partial \mathbf{p}},\frac{d\mathbf{p}}{dt}=-\frac{\partial H}{\partial \mathbf{r}}$$

Now, we know that the (idealized) magnetic field of the system is given by:

$$\mathbf{B}=B_y\hat{y}=B_0\sin (k_{u}z)\hat{y}$$

Therefore by working backwards, we have (up to the addition of an irrotational field):

$$\mathbf{A}=-A_0\cos (k_{u}z)\hat{x},A_0\equiv \frac{B_0}{k_u}$$

Taking the curl it is easy to see that $\nabla \times \mathbf{A}=B_0\sin (k_{u}z)\hat{y}$, as we started with. This gives us the following Hamiltonian:

$$H(\mathbf{r},\mathbf{p},t)=\frac{1}{2m}(\mathbf{p}-e{A}_0\cos (k_{u}z)\hat{x}{)}^2$$

Note: since $H$ does not depend on $x,y$, Hamilton’s equations immediately tell us that $\ddot{x},\ddot{y}=0$.

We can explicitly evaluate this Hamiltonian by substituting in $\mathbf{r}$ and $\mathbf{p}$, giving us a regular system of first-order ODEs that can be solved by any ODE solver. However, there are also specialized libraries like pyhamsys that can directly numerically-integrate a given Hamiltonian with symplectic integrators, ensuring energy stability throughout the simulation.

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.}$

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.

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)