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. ![[undulator-trajectory.excalidraw.svg|400]] The electron beam trajectory $\mathbf{r}(t)$ follows the relativistic Lorentz force equation for magnetic fields: $ \dfrac{d\mathbf{p}}{dt} = \dfrac{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. ### 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{eB_0 L_u}{2\pi mc^2}$. This definition of $K$ is clearly in Gaussian units due to his extra factor of $c$, since in SI we have $K = \frac{eB_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' \equiv \omega_0/\gamma$ we can write this more elegantly as $a = \frac{\omega_0'}{k_u^2 c}$. Thus (in our notation): $ x \approx a \sin (k_u z) = \frac{\omega_0'}{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'$. Since $\omega_0' \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 c t$ 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, \quad 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{align*} \dfrac{d}{dt}(\gamma m \dot x) = q(\cancel{v_y B_z}^0 - v_z B_y) \\ \dfrac{d}{dt}(\gamma m \dot y) = q(\cancel{v_z B_x}^0 - \cancel{v_x B_z}) \\ \dfrac{d}{dt}(\gamma m \dot z) = q(v_x B_y - \cancel{v_y B_x}^0) \end{align*} $ Thus, we find that the electrons move purely in the $xz$-plane, and their equations of motion are given by: $ \begin{align*} \dfrac{d}{dt}(\gamma m \dot x) &= -q v_z B_y \\ \dfrac{d}{dt}(\gamma m \dot z) &= qv_x B_y \end{align*} $ 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: $ \dfrac{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[^4]: $ \dot \gamma \dfrac{dx}{dt} - \gamma \dfrac{d^2 x}{dt^2} = \omega_0 v_0 \sin(k_u v_0 t), \quad \omega_0 \equiv qB_0/m $ Where we once again define $\omega_0 \equiv qB_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: $ \dfrac{d^2 x}{dt^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' = \omega_0/\gamma$ as the "relativistic" form of $\omega_0$ gives us: $ x = -\frac{\omega_0'}{v_0 k_u^2} \sin k_u z = a \sin k_uz, \quad a =-\frac{\omega_0'}{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'}{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{align*} \dfrac{d}{dt}(\gamma m \dot x) &= -q \dot z B_y \\ \dfrac{d}{dt}(\gamma m \dot z) &= q\dot x B_y \end{align*} $ 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{align} \ddot x &= - \dot z qB_y/m \\ \ddot z &= \dot x qB_y/m \end{align} $ Or in vector form: $ \frac{d}{dt} \begin{pmatrix} \dot x \\ \dot z \end{pmatrix} = \frac{q B_y}{m}\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \dot x \\ \dot z \end{pmatrix} $ 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} \begin{pmatrix} \dot x \\ \dot z \end{pmatrix} = \omega_0\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \dot x \\ \dot z \end{pmatrix} $ 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: ```python # 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{align*} 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{align*} $ 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: $ \begin{pmatrix} \ddot x \\ \ddot z \end{pmatrix} = \sin(k_u z)\begin{pmatrix} 0 & -\omega_0 \\ \omega_0 & 0 \end{pmatrix} \begin{pmatrix} \dot x \\ \dot z \end{pmatrix} $ The system can then be rewritten in the following (phase-space) form, where $\mathbf{y} = (x, z, v_x, v_z)$: $ \dot{\mathbf{y}} = \begin{pmatrix} \dot x \\ \dot z \\ \dot v_x \\ \dot v_z \end{pmatrix} = \begin{pmatrix} 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{pmatrix} \begin{pmatrix} x \\ z \\ v_x \\ v_z \end{pmatrix} $ 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^* = \frac{\pi}{2k_u}$. Now expanding about this steady state is equivalent to evaluating the Jacobian at $z^* = \frac{\pi}{2k_u}$. The Jacobian in this case is given (in general) by: $ J(x, z, v_x, v_z) = \begin{pmatrix} 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{pmatrix} $ Evaluating at $z^* = \pi/2k_u$ gives us: $ J(z^*) = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & -\omega_0 \\ 0 & 0 & \omega_0 & 0 \end{pmatrix} $ The linearized version of the differential equation is then given by $\dot{\mathbf{y}} = J(z^*)\mathbf{y}$. The solution can be found by first finding the eigenvalues and eigenvectors of the Jacobian $J(z^*)$ (which we'll notate as $J^*$ 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: ```python 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: $ \begin{bmatrix} x(t) \\ z(t) \\ v_x(t) \\ v_z(t) \end{bmatrix} = \begin{bmatrix} C_1 \\ C_2 \\ 0 \\ 0 \end{bmatrix} + C_3\begin{bmatrix} 1 \\ i \\ -i\omega_0 \\ \omega_0 \end{bmatrix} e^{-i\omega_0 t} + C_4 \begin{bmatrix} 1 \\ -i \\ i\omega_0 \\ \omega_0 \end{bmatrix} e^{i\omega_0 t} $ Simplifying with some new constants $\tilde C_3, \tilde C_4$ gives us: $ \begin{align} x(t) = C_1 + \tilde C_3\cos \omega_0 t \\ z(t) = C_2 + \tilde C_4 \sin \omega_0 t \end{align} $ 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](http://egun-igun.com/), 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{align} \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{align} $ 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{align} \ddot x &= - \dot z qB_y/m \\ \ddot z &= \dot x qB_y/m \end{align} $ Or in vector form: $ \frac{d}{dt} \begin{pmatrix} \dot x \\ \dot z \end{pmatrix} = \frac{q B_y}{m}\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \dot x \\ \dot z \end{pmatrix} $ 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](https://en.wikipedia.org/wiki/Nondimensionalization). To begin, the (non-relativistic) Lorentz force law with zero electric field can be written in the form: $ \ddot x^i = \dfrac{qB_\perp}{m} \dot x^j = \omega_0 f(x^k) \dot x^j, \quad 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: $ \dfrac{d}{dt} = \dfrac{d\tau}{dt}\frac{d}{d\tau} = \frac{1}{t_c}\frac{d}{d\tau}, \quad \dfrac{dx^j}{dt} = \frac{x_c}{t_c} \dfrac{d\chi^j}{d\tau} $ Likewise: $ \dfrac{d^2 x^i}{dt^2} = \dfrac{d}{dt} \left(\dfrac{dx^i}{dt}\right) = \frac{x_c}{t_c^2} \dfrac{d^2\chi^i}{d\tau^2} $ Substituting back into the Lorentz force law gives us: $ \frac{x_c}{t_c^2} \dfrac{d^2\chi^i}{d\tau^2} = \omega_0 f(x_c \chi^k) \frac{x_c}{t_c} \dfrac{d\chi^j}{d\tau} $ Simplifying and rearranging, we have: $ x_c\dfrac{d^2\chi^i}{d\tau^2} = t_c\omega_0 x_c f(x_c \chi^k) \dfrac{d\chi^j}{d\tau} $ Now, let us define $\beta^k = \frac{1}{c} \frac{d\chi^j}{d\tau}$. We can thus write the above as: $ x_c\dfrac{d^2\chi^i}{d\tau^2} = t_c\omega_0 x_c f(x_c \chi^k) \dfrac{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 in theory somewhat arbitrary since it is divided from both sides. However, note that for large but not highly-relativistic velocities $0.01c \leq v \leq 0.5c$, there is another choice for $t_c$ that is more useful for numerical work. If we define $t_c = 1/c\omega_0$, then: $ \dfrac{d^2\chi^i}{d\tau^2} = \frac{1}{c} f(x_c \chi^k) \dfrac{d\chi^j}{d\tau} $ In doing so, we can express the velocities in terms of the dimensionless parameter $\beta = v/c$ (as is commonly used in special relativity), giving us: $ \dfrac{d^2\chi^i}{d\tau^2} = f(x_c\chi^k) \dfrac{d\beta^j}{d\tau} $ Nothing qualitatively changes about the differential equation here; the only difference is that our velocities are now expressed as a percent of the speed of light. This has the advantage of keeping our velocities within the range $0 \leq \beta < 1$, instead of being massive numbers on the order of $v \sim 10^8$! With all of these substitutions, our system may be written as: $ \dfrac{d^2}{d\tau^2} \begin{pmatrix} \tilde x \\ \tilde z \end{pmatrix} = \sin(x_ck_u\tilde z) \begin{pmatrix} -\tilde \beta_z \\ \tilde \beta_x \end{pmatrix} $ Where $\chi^j = (\tilde x, \tilde z)$ and $t = \frac{1}{\omega_0 c}\tau$. We can choose $x_c = 1/k_u$ to remove all dimensions from the differential equation, so that it takes a particularly-simple form: $ \dfrac{d^2}{d\tau^2} \begin{pmatrix} \tilde x \\ \tilde z \end{pmatrix} = \sin(\tilde z) \begin{pmatrix} -\tilde \beta_z \\ \tilde \beta_x \end{pmatrix} $ However, it turns out that this choice actually works rather poorly for numerical solving for lower velocities. It is better to leave $x_c = 1$, giving us: $ \dfrac{d^2}{d\tau^2} \begin{pmatrix} \tilde x \\ \tilde z \end{pmatrix} = \sin(k_u\tilde z) \begin{pmatrix} -\tilde \beta_z \\ \tilde \beta_x \end{pmatrix} $ At the cost of needing one non-dimensionless parameter, this variant of the system is much better-suited to solving. 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. > **Note:** The code for the numerical solver can be found in `Oscars/OscarsBasicSimulation.ipynb` #### Approximate solution with SciPy 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}),\quad f(\mathbf{U}) =\begin{pmatrix} v_x \\ v_z \\ -\omega_0 \sin(k_u z)\dot z \\ \omega_0 \sin(k_u z) \dot x \end{pmatrix} = \begin{pmatrix} U_2 \\ U_3 \\ -\omega_0\sin(k_u U_1) U_3 \\ \omega_0 \sin(k_u U_1) U_2 \end{pmatrix} $ #### 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 \dfrac{dz}{dt} - \gamma \dfrac{d^2 z}{dt^2} = \dfrac{qB_0}{m} \dot x \sin k_0 x $ Which can be rearranged to yield: $ \dfrac{d^2 z}{dt^2} = \dfrac{1}{\gamma}\left[\dot \gamma \dfrac{dz}{dt} - \dfrac{qB_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{align*} \dfrac{dx}{dt} &= v_0, \quad \dfrac{dz}{dt} = v_z, \\ \dfrac{dv_x}{dt} &= 0, \quad \dfrac{dv_z}{dt} = \dfrac{1}{\gamma}\left[\dot \gamma v_z - \dfrac{qB_0}{m} v_x \sin k_0 x\right] \end{align*} $ Note that $\dot \gamma$, the time derivative of the Lorentz factor, is given by: $ \begin{align*} \dfrac{d\gamma}{dt} &= \dfrac{d}{dt}\left(\dfrac{1}{\sqrt{1 - v^2/c^2}}\right) \\ &= \dfrac{v}{c^2(1 - v^2/c^2)^{3/2}} \dfrac{dv}{dt} \\ &= \dfrac{v}{c^2(1 - v^2/c^2)^{3/2}} |\mathbf{a}|, \quad v=|\mathbf{v}| = \sqrt{\dot x^2 + \dot z^2}, \quad |\mathbf{a}| = \sqrt{\dot v_x^2 + \dot v_z^2} \end{align*} $ 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: $ \dfrac{d\gamma}{dt} = \dfrac{v}{c^2(1 - v^2/c^2)^{3/2}}\dfrac{dv_z}{dt} $ We note that since $\gamma = \dfrac{1}{(1 - v^2/c^2)^{1/2}}$ we can also write our result as follows: $ \dfrac{d\gamma}{dt} = \dfrac{v}{c^2} \gamma^3 \dfrac{dv_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 $\dfrac{dv_z}{dt}$). To see why, we can expand the 4th ODE using our definition of $\dot \gamma$, giving us: $ \dfrac{dv_z}{dt} = \dfrac{v}{c^2} v_z\gamma^2 \dfrac{dv_z}{dt} - \dfrac{1}{\gamma} \dfrac{qB_0}{m} v_x \sin k_0 x $ By slight rearrangement, we can resolve our issue by moving the $\dfrac{dv_z}{dt}$ term to the left-hand side, giving us (after multiplying again by $\gamma$): $ \begin{gather*} \dfrac{dv_z}{dt}\left[\dfrac{v}{c^2} v_z\gamma^2 - 1\right] = \dfrac{1}{\gamma} \dfrac{qB_0}{m} v_x \sin k_0 x \\ \gamma\dfrac{dv_z}{dt}\left[\dfrac{v}{c^2} v_z\gamma^2 - 1\right] = \dfrac{qB_0}{m} v_x \sin k_0 x \\ \dfrac{dv_z}{dt} = \dfrac{v_x}{v\gamma^3 v_z /c^2 - \gamma}\dfrac{qB_0}{m} \sin k_0 x \end{gather*} $ Thus, our system of ODEs becomes: $ \begin{align*} \dfrac{dx}{dt} &= v_0, \quad \dfrac{dz}{dt} = v_z, \\ \dfrac{dv_x}{dt} &= 0, \quad \dfrac{dv_z}{dt} = \dfrac{v_x}{v\gamma^3 v_z /c^2 - \gamma}\dfrac{qB_0}{m} \sin k_0 x \end{align*} $ Which can be written in the matrix form $\mathbf{y}' = f(\mathbf{y}, t)$ as follows: $ \mathbf{y} = \begin{pmatrix} x \\ z \\ v_x \\ v_z\end{pmatrix}, \quad f(\mathbf{y}, t) = \begin{pmatrix} v_0 \\ v_z \\ 0 \\ \dfrac{v_x}{v\gamma^3 v_z /c^2 - \gamma}\dfrac{qB_0}{m} \sin k_0 x \end{pmatrix} $ We can translate this system to Python code as follows: ```python 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{cases} \nabla \cdot (\mu(-\nabla \psi_m) + \mathbf{M}) = 0 \\ \mathbf{B} = \mu\mathbf{H} = -\mu(-\nabla\psi_m + \mathbf{M}) \end{cases} $ > **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[^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''(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: $ \dfrac{d}{dt}(\gamma m \mathbf{v}) = q(\mathbf{E} + \mathbf{v} \times \mathbf{B}) $ In the non-relativistic regime, we have: $ \dfrac{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 \varepsilon_0}\dfrac{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} mv_x^2 \geq 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](https://arborjs.org/docs/barnes-hut) in electrostatics and the [multigrid method for Poisson's equation](https://vuduc.org/teaching/cse8803-pna-sp08/slides/cse8803-pna-sp08-12.pdf). A promising alternative is the [particle-in-cell method](https://www.particleincell.com/2011/velocity-integration/), 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 \varepsilon_{0} m_{e}(v_{z} \gamma)^3} $ Where $e$ is the magnitude of the elementary charge, $\varepsilon_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) = \dfrac{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{gather} 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{gather} $ This comes from the leading-order term in the series expansion for $F(\chi)$, which is a [Puisseux series](https://en.wikipedia.org/wiki/Puiseux_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](https://en.wikipedia.org/wiki/Incomplete_gamma_function), given by: $ \Gamma(s, x) = \int_{x}^{\infty} t^{s-1}e^{-t}d t $ 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)\Bigg|_{0}^b $ Note that by the [properties of the incomplete gamma function](https://en.wikipedia.org/wiki/Incomplete_gamma_function#Special_values), $\Gamma\left( \frac{1}{2}, x \right) = \sqrt{ \pi } \operatorname{erfc}(\sqrt{ x })$ where $\operatorname{erf}(x)$ is the [error function](https://en.wikipedia.org/wiki/Error_function) and $\operatorname{erfc}(x) = 1 - \operatorname{erf}(x)$ is the [complementary error function](https://mathworld.wolfram.com/Erfc.html) Thus the above reduces to: $ -i \,\Gamma \left( \frac{1}{2}, -x \right)\Bigg|_{0}^b = -i\sqrt{ \pi } \operatorname{erfc}(\sqrt{ -x })\Bigg|_{0}^b = -i\sqrt{ \pi }(\operatorname{erfc}(\sqrt{ -b }) - 1) = i\sqrt{ \pi }\operatorname{erf}(\sqrt{ -b }) $ The Maclaurin series of the error function is given by: $ \begin{align*} \operatorname{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{align*} $ Thus we have: $ \operatorname{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 }\operatorname{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{align*} 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{align*} $ 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} + q V $ This gives you Hamilton's equations, a system of two first-order differential equations, which are respectively given by: $ \begin{align*} \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{align*} $ 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{align*} \frac{d\mathbf{R}}{dt} &= \mathbf{v} \\ \frac{d \mathbf{P}}{dt} &= q\mathbf{v} \times (\nabla \times \mathbf{A}) \end{align*} $ 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{align} \dot x &= P_x/m \\ \dot y &= P_y/m \\ \dot P_x &= -\frac{qP_y}{m}\frac{f(\theta)}{r} \\ \dot P_y &= \frac{qP_x}{m}\frac{f(\theta)}{r} \end{align} $ 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](https://oscars.bnl.gov/) (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](https://codeberg.org/elaraproject/laser-research/src/branch/main/Oscars). ## References - OSCARS reference examples - [Magnetic field modelling](https://oscars.bnl.gov/examples/AllAbout_MagneticFields.html) (also includes information on how to load magnetic fields from external softwares) - [Undulator electron trajectory simulations](https://oscars.bnl.gov/examples/Example_002_UndulatorTrajectory.html) - [Particle beam modelling](https://oscars.bnl.gov/examples/AllAbout_ParticleBeams.html) (this also includes multi-particle beam modelling, which is very important for our use case) [^1]: From the [FreeFEM tutorial on simulating a microwave oven](https://doc.freefem.org/tutorials/complexNumbers.html). Note that they use $\theta$ instead of $T$ as the symbol for the temperature, and write the Laplacian as $\Delta$ rather than $\nabla^2$. The quartic $\mu(T^4 - T^4_\text{out})$ term comes from the [Stefan-Boltzmann law](https://en.wikipedia.org/wiki/Heat_equation#Accounting_for_radiative_loss) and models heat loss due to thermal radiation. [^2]: Note that $\sigma$ has units of $\pu{C^2*K*J^{-2}}$ (Coulomb-squared-Kelvin per Joule-squared). [^3]: No relationship between $\sigma$ and the electrical conductivity/charge density, which use the same symbol [^4]: 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. [^5]: 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](https://en.wikipedia.org/wiki/Fermi%E2%80%93Dirac_statistics#Distribution_of_particles_over_energy) given that electrons are fermions. [^6]: 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). [^7]: This result comes from [the following COMSOL tutorial](https://www.comsol.com/model/relativistic-diverging-electron-beam-17065), which itself references the book _[Charged Particle Beams](https://books.google.com/books/about/Charged_Particle_Beams.html?id=1GjCAgAAQBAJ&source=kp_book_description)_ by Stanley Humphries. Also, the [Wikipedia article on charged particle beams](https://en.wikipedia.org/wiki/Charged_particle_beam) may be useful reading. [^8]: As suggested by [this Physics forum article](https://www.physicsforums.com/threads/antiderivative-of-1-sqrt-lnx-c.83520/). [^9]: Some other texts (including [this document from CERN's accelerator school](https://cas.web.cern.ch/sites/default/files/lectures/senec-2012/chauvin.pdf), 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$. [^10]: Taken from the [lecture notes from the US particle accelerator school](https://uspas.fnal.gov/materials/14Knoxville/Lecture01.pdf) - the equation is from there (pg. 39 and 40) [^11]: See [this Wikipedia article](https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates#Del_formula), which lists the vector differential operators in Cartesian, cylindrical, and spherical coordinates. [^12]: Discussed on pg. 23 of the [lecture notes from the US particle accelerator school](https://uspas.fnal.gov/materials/14Knoxville/Lecture01.pdf)