This page gives a general description of how lasers and masers are described and simulated for Project Elara. > **Note:** For specific modelling of free-electron laser components, including electron guns and undulators, please see [[Numerical modelling of electron beams]]. ## Physical foundations of lasers All lasers, including free-electron lasers (and masers), produce electromagnetic waves that are tightly collimated into a beam. This beam is mathematically described by the **Gaussian beam**, which is a solution to Maxwell's equations that is given (in cylindrical coordinates by: $ \mathbf{E}(r, \phi, z) = E_0 \frac{w_0}{w(z)} \exp\left(-\frac{r^2}{w(z)^2}\right) \exp \left(-i\left(kz + \frac{kr^2}{2R(z)}\right)\right) \hat{\mathbf{r}} $ The Gaussian beam represents the beam produced by an idealized laser, and can be derived from some basic physical assumptions along with physically-imposed boundary conditions (here whenever we say "lasers" we mean "lasers and masers"). This holds for (nearly) all lasers, regardless of their mode of operation, since it is a product of the specific *boundary conditions* imposed by the design of **laser cavities**, which have a large number of commonalities across different designs. In particular, the goal of a laser cavity is to achieve **resonance**: the amplification of electromagnetic energy along a specific axis to create a powerful, highly-collimated beam. For lasers that operate in the visible light range, this is also known as **optical resonance**. This is typically accomplished by means of two mirrors, a fully-reflective (also called _silvered_) mirror, and a partially-reflective (also called _semi-silvered_) mirror. The two mirrors are connected by a tube that allows light (or microwaves for masers) to pass through unimpeded. Within the tube, it is common (though not always) to find a material that fills the interior of the tube: this is called the _gain medium_, and it can be everything from a gas (carbon dioxide, helium, neon, and hydrogen are all possibilities) to a semiconductor (in [laser diodes](https://en.wikipedia.org/wiki/Laser_diode)) to some type of plastic (for certain types of masers). The gain medium acts as both a source that *generates* the electromagnetic field and a material that *amplifies* the electromagnetic field within the laser cavity: it is said to be _pumped_ when it receives (is "pumped with") energy from an external source, like an electric arc, flash lamp, or even another laser. A general diagram of a laser cavity is shown below: ![[laser-cavity.excalidraw.svg]] The partially reflective mirror allows some of the electromagnetic waves in the cavity to escape while reflecting the rest of the rest back into the cavity, ensuring that energy is continually built up within the cavity to further amplify the field. The escaped waves are what we observe as a characteristic laser beam, and gives lasers their distinctive properties.  _An image of a He-Ne laser, a common laser type. The walls of its optical cavity are glass (transparent), while two mirrors at the end achieve optical resonance necessary for amplification._ To be able to model lasers (and masers) effectively, it is necessary to have an in-depth understanding of the Gaussian beam solution, and how it arises from the boundary conditions imposed by the physical construction and materials used within the laser cavity. We will describe this, and more, in depth throughout the guide. ### Standard derivation of the Gaussian beam The most basic derivation of the Gaussian beam starts by assuming that the electric field is transverse (perpendicular) to the axis of propagation, and takes the form: $ E(\mathbf{r}) = E_0 \hat{\mathbf{x}}~u(\mathbf{r}) e^{ikz} $ Where $z$ is the direction the beam propagates (the *optical axis* or *transmission axis*). By substituting this into the Helmholtz equation, we get the **paraxial Helmholtz equation**[^1]: $ \left(\dfrac{\partial^2}{\partial x^2} + \dfrac{\partial^2}{\partial y^2}\right)u + 2ik \dfrac{\partial u}{\partial z} = 0 $ Note that if we define $\nabla^2_{(2D)} = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ the paraxial Helmholtz equation can also be written in a more familiar form: $ \nabla^2_{(2D)}u + 2ik \dfrac{\partial u}{\partial z} = 0 $ We can now solve the paraxial Helmholtz equation. Since laser cavities are typically cylindrical, it is conventional to use cylindrical coordinates $(r, \phi, z)$ instead of Cartesian coordinates $(x, y, z)$. If we transform the PDE to cylindrical coordinates, the paraxial Helmholtz equation reads: $ \frac{1}{r}\frac{\partial}{\partial r}\left(r \dfrac{\partial}{\partial r}\right) + \frac{1}{r^2} \dfrac{\partial^2 u}{\partial \phi^2} + 2ik \dfrac{\partial u}{\partial z} = 0 $ Now, by assuming that the beam has azimuthal symmetry (that is, it doesn't depend on the azimuthal angle $\phi$), the derivative with respect to $\phi$ drops out, giving us: $ \begin{gather*} \frac{1}{r}\frac{\partial}{\partial r}\left(r \dfrac{\partial}{\partial r}\right) + 2ik \dfrac{\partial u}{\partial z} = 0 \\ \Rightarrow \frac{1}{r}\dfrac{\partial u}{\partial r} + \frac{\partial^2 u}{\partial r^2} + 2ik \dfrac{\partial u}{\partial z} = 0 \end{gather*} $ To derive the Gaussian beam solution in particular, we assume that $u(\mathbf{r})$ takes the form of a Gaussian profile[^2]: $ u(\mathbf{r})= f(\mathbf{r}) \exp(-\alpha r^2), \quad \alpha =\dfrac{1}{w(z)^2} + \dfrac{ik}{2R(z)} $ Where $w(z)$ describes the *beam width*, and $R(z) \sim z\left(1 + \dfrac{1}{z^2}\right)$ describes the *beam curvature*. It is possible to prove *why* the beam must take this form, though it is not exactly a simple derivation. A heuristic is that since we know that any physical EM wave must go to zero as $r \to \infty$ and $z \to \infty$, we must pick an appropriate function that decays smoothly to infinity, so a Gaussian makes sense. Meanwhile, $\alpha$ must be a function of $z$ to satisfy the requirement that a beam must diverge over space as per the predictions of electromagnetic theory. We can then work backwards to find the correct $\alpha(z)$ such that the beam diverges *exactly* at the diffraction limit, giving us the Gaussian beam solution. ### Alternative derivation of the Gaussian beam If we don't start with any approximation and use the standard Helmholtz equation $(\nabla^2 + k^2)u = 0$, we can also derive the Gaussian beam, albeit through a much more convoluted process. The idea is to first do separation of variables in cylindrical coordinates. Assuming some basic boundary conditions (namely that the solution doesn't explode at zero or infinite distance), this gives you the general solution in terms of $J_n(x)$, which are the [Bessel functions](https://en.wikipedia.org/wiki/Bessel_function) of the first kind. The general solution reads as follows: $ u(r, \phi, z) = \sum_m \sum_n A_{mn} J_n(kr)e^{im\phi} e^{-k_z z} $ Here, the $A_{mn}$ coefficients are unknown and must be determined. A more rigorous and step-by-step derivation can be found in [[Solving for the fields in a hollow waveguide]]. Now, if we assume azimuthal symmetry (symmetry about $\phi$) and assume continuous coefficients (such that the sum becomes an integral), we have: $ u(r, \phi, z) = \int dk~A(r) J_n(kr)e^{im\phi} e^{-k_z z} $ Where we choose the $m = 0, n=0$ mode since it is the lowest-order (and thus the dominant) mode. Factoring out everything that does not explicitly depend on $k$ ($k_z$ is independent of $k$) we have: $ u(r, \phi, z) = e^{-k_z z}\int dk~A(r) J_0(kr) $ Now, to get the Gaussian beam solution, we assume that the beam is symmetric and its energy is concentrated along the optical axis, so it can be modelled by a Gaussian function in the form: $ u(r, \phi, 0) = C e^{-ar^2} $ Where $C$ and $a$ are some undetermined constants. This gives us: $ u(r, \phi, z) = C e^{-k_z z}\int k\, dk~e^{-ar^2} J_0(kr) $ Here, the extra factor of $k$ in $k dk$ comes from the fact that the integral $\displaystyle \int k\,dk~e^{-\alpha r} J_0(kr)$ is an (inverse) **Hankel transform**. This particular transform has a well-known solution in the form $e^{-br^2}$, where $b$ is a constant that is related to $a$. This gives us: $ u(r, \phi, z) = C e^{-b r^2}e^{-k_z z} = Ce^{-(br^2 + k_z z)} $ This is not _quite_ the correct solution but it is pretty close, and if we choose our Gaussian _ansatz_ more carefully, then we can derive the exact solution.[^7] ### Applicability of the Gaussian beam solution Ideally, the electric field of any laser - including a free-electron laser - would be both *internally* and *externally* described by Gaussian beam. A more accurate representation that accounts for the imperfections in a Gaussian beam is that the electric field is described by a superposition of [Hermite-Gaussian modes](https://en.wikipedia.org/wiki/Gaussian_beam#Hermite-Gaussian_modes), which are each labelled $\text{TEM}_{lm}$, where $l, m$ are integers. The lowest mode ($\text{TEM}_{00}$) is the ideal case and corresponds to the Gaussian beam, while higher-order modes describe small correction to describe realistic laser/maser beams. However, in practice, *numerically* obtaining this result is not as simple. First, we must impose periodic boundary conditions azimuthally (that is, $E(r, \phi) = E(r, \phi + 2\pi)$) with $E(r) \to 0$ as $r \to \infty$. Since the electric field and its derivative are non-vanishing on the semi-silvered side of the laser cavity, neither Dirichlet or Neumann boundary conditions would suffice; rather, it would be necessary to use a complicated *partially-absorbent layer* to describe the fact that some light can escape from the semi-silvered side. This requires a much more complicated model, which we'll discuss next. ## Boundary-value problem for a laser cavity To solve for an *exact* solution for the electromagnetic field both inside and outside a laser cavity, we cannot just _a priori_ assume the Gaussian beam solution. Rather, we must solve the full **boundary-value problem** of the Helmholtz equation; crucially, we will not use the paraxial approximation here, as we want to find an exact solution. To begin, we assume a laser cavity that is cylindrical, with radius $R$ and length $L$ - this is the most typical type of laser cavity and used on a variety of lasers and masers. We must now give the boundary conditions of the problem. First, a fully-reflective mirror can be modelled by a perfect conductor, which implies that the tangential component of the electric field along the surface of the mirror is zero. Second, the gain medium can be modelled by a dielectric with dielectric constant $\epsilon_r$ (also known as the _relative permittivity_). Note that we assume a linear and homogeneous gain medium; otherwise the dielectric constant is no longer constant, and is in general a tensor with 9 components (but considering the full tensorial form of $\epsilon_r$ would make finding an analytical solution to the problem near-impossible). Third, a partially-reflective mirror can be modelled by a dielectric with dielectric constant $\epsilon_{r}'$. By the standard electrodynamics boundary conditions, the electric fields $\mathbf{E}_\mathrm{interior}$ and $\mathbf{E}_\mathrm{mirror}$ must satisfy: $ E_\mathrm{interior}^\parallel = E_\mathrm{mirror}^\parallel, \quad \epsilon_{r} E_\mathrm{interior}^\perp = \epsilon_{r}' E_\mathrm{mirror}^\perp $ Fourth, we know that the beam has azimuthal symmetry, and therefore we naturally have the periodic boundary condition $E_\phi(r, \phi, z) = E_\phi(r, \phi + 2\pi n, z)$. Fifth, the electromagnetic field in the $z$ direction (that is, the $E_z$ component of the field) is necessarily zero, since the laser beam points along $z$ (and thus the wave-vector $\mathbf{k}$ must be in the form $\mathbf{k} = |\mathbf{k}|\hat z$), but the electric field is always transverse to the direction of the wave-vector (that is, $\mathbf{E} \times \mathbf{k} = 0$), and thus $E_z$ must be zero in our case. Thus we only need to solve for $E_\phi$ and $E_r$ (the radial and azimuthal components of the field); moreover, $E_\mathrm{interior}^\parallel = E_\mathrm{mirror}^\parallel = 0$. (Note that due to the azimuthal of the beam, the $E_r$ component is of far more interest to us than $E_\phi$). And finally, the tube between the two mirrors is meant to be an open boundary and allow electromagnetic waves to freely pass through. Thus we have a _radiative boundary condition_ (also known as the Sommerfeld boundary condition): $ \lim_{ r \to \infty } r\left( \frac{\partial}{\partial r} - ik \right)E_{r} = 0, \quad 0 < z < L $ Since we are solving the *exact* Helmholtz equation rather than the paraxial Helmholtz equation, the solution will not be a Gaussian beam, although we can apply the paraxial approximation later on to show that the $E_r$ component of the field *reduces* to a Gaussian beam in the paraxial limit. To solve, we will first write out the Helmholtz equation for each of the components of the electric field: $ \begin{align*} \nabla^2 E_{r} + k^2 E_{r} = 0 \\ \nabla^2 E_{\phi} + k^2 E_{\phi} = 0 \end{align*} $ Now, upon performing separation of variables in cylindrical coordinates, we have a general solution in the following form: $ u(r, \phi, z) = (A \cos (m\phi) + B \sin(m\phi))(CJ_{n}(kr) + DY_{n}(kr)) $ We thus have the ansatz (explain): $ E_{r}^{m,n}(r , \phi, z) \sim A_{mn} \cos(m \phi) J_{n}(k_{r}r) e^{-ik_{z} z}, \quad k_{z} = \frac{n\pi z}{L} $ Check by substitution into the Helmholtz equation that it satisfies the Helmholtz and satisfies the BCs. Also solve for resonant frequencies. Note that if we choose the $m = 0$ mode, the cosine vanishes and the field no longer depends on $\phi$: $ E_{r}^{0,n}(r , \phi, z) \sim A_{0n} J_{n}(k_{r}r) e^{-ik_{z} z}, \quad k_{z} = \frac{n\pi z}{L} $ If we then make the approximation that $J_0(k_r r) \approx e^{-r^2}$ (a valid approximation if $k_r$ is small), then our solution strongly resembles the Gaussian beam. ### Numerical solution To start a basic numerical model, we can start by once again assuming azimuthal symmetry implicitly, which reduces the 3D problem to a 2D problem, where we use the coordinates $(r, \phi, z)$ to denote the transverse (perpendicular to optical axis) and longitudinal (parallel to optical axis) directions, so $E_\perp = E_\perp (r, z)$ and $E_\parallel = E_\parallel(r, z)$. The total electric field is then given by $\mathbf{E} = E_\perp \hat r + E_\parallel \hat z$. A typical laser cavity has a 100% reflective mirror on its left end (at $z = 0$) and a partially-reflective mirror with reflectivity $\eta < 1$ at its right end (it is a bit different for masers but we'll discuss that later). The electric field satisfies the [Helmholtz equation](https://en.wikipedia.org/wiki/Helmholtz_equation) $(\nabla^2 + k^2)\mathbf{E} = 0$, which are essentially the time-independent Maxwell equations. When expanded into vector form, this gives us two PDEs to solve: $ \begin{align*} \dfrac{\partial^2 E_\perp}{\partial r^2} + \dfrac{\partial^2 E_\perp}{\partial z^2} + k^2 E_\perp = 0 \\ \dfrac{\partial^2 E_\parallel}{\partial r^2} + \dfrac{\partial^2 E_\parallel}{\partial z^2} + k^2 E_\parallel = 0 \end{align*} $ As mentioned, the longitudinal field runs along the optical axis (that is, $E_\perp = E_z$), while the transverse field runs perpendicular to it (that is, $E_\parallel = E_r$). We assume that the left end (at $z = 0$) of the laser cavity has a perfectly-reflective mirror, such that $E_\perp = 0$ (the tangential component of the electric field is zero, which is true for all perfect conductors). This is equivalent[^5] to the following boundary condition: $ \hat{\mathbf{n}} \times \mathbf{E} = 0 $ Where in this case we have $\hat{\mathbf{n}} = \hat{\mathbf{z}} = (0, 1)^T$. Now, we examine the right side of the laser cavity (at $z = L$). For this side, we have a *partially-transparent* mirror. This means that the mirror will reflect some light and also transmit some light (we assume no absorption here). There are several different approaches to this. The first method, and a fairly naive one[^6], is to start with the general form of a scattered (plane) wave: $ u(\mathbf{x}) = e^{ik\cdot \mathbf{x}} + r e^{-ik\cdot\mathbf{x}} $ Where $r$ is the *reflection coefficient* and $R = |r|^2$ is the *reflectance*, or essentially the percentage of light that is reflected. We now consider a generalized mixed boundary condition: $ c_1 u(\mathbf{x}) + c_2 \dfrac{\partial u}{\partial n} + c_3 = 0 $ Substituting our solution $u(\mathbf{x})$ in, we have: $ \begin{gather*} c_1 u(\mathbf{x}) + c_2 (ik e^{ik \cdot \mathbf{x}} - ikr e^{-ik\cdot \mathbf{x}}) = -c_3 \\ c_1 (e^{ik\cdot \mathbf{x}}+ re^{-ik\cdot \mathbf{x}}) + i k c_2 (e^{ik\cdot \mathbf{x}} - re^{-ik\cdot \mathbf{x}}) = -c_3 \\ (c_1 + ik c_2) e^{ik\cdot \mathbf{x}} + r(c_1 - ik c_2 )e^{-ik\cdot \mathbf{x}} = -c_3 \\ r = -\dfrac{c_3 + (c_1 - ik c_2)e^{ik\cdot \mathbf{x}}}{(c_1 + ik c_2)e^{-ik\cdot \mathbf{x}}} \end{gather*} $ Thus, we have: $ \begin{align*} R &= |r|^2 \\ &= \left|\dfrac{c_3 + (c_1 - ik c_2)e^{ik\cdot \mathbf{x}}}{(c_1 + ik c_2)e^{-ik\cdot \mathbf{x}}}\right|^2_{(\mathbf{x} = \partial \Omega = L)} \\ &= \left|1 + \dfrac{c_3}{(c_1 + ik c_2)e^{-ikL}}\right|^2 \\ \end{align*} $ This result is very interesting since it tells us that $c_2 < 0$ in order for a physically-meaningful reflection coefficient ($R>1$ would violate conservation of energy). Now, the particular values of $c_1, c_2, c_3$ to make this equation true for a certain value of $R$ are actually *arbitrary*. For simplicity, we can set $c_2 = 0$ so that we can reduce one of our constants, giving us (for *different* $c_1$ and $c_3$ than previously): $ R = \left|1 + \dfrac{c_3}{c_1}e^{ikL}\right|^2 $ Now, let us presume that $c_3 = \alpha c_1$, since again the choice of our constants is arbitrary so long as our equation is fulfilled, which also requires that $\alpha < 0$. Thus we have: $ \begin{align*} R &= \left|1 + \dfrac{\alpha c_1}{c_1}e^{ikL}\right|^2 \\ &= \left|1 + \alpha e^{ikL}\right|^2 \\ &= 1 + \alpha^2 + 2 \alpha \cos(kL) \end{align*} $ The solution can be found using the quadratic formula: $ \begin{gather*} \alpha^2 + 2\alpha \cos (kL) + (1 - R) = 0 \\ \Rightarrow \alpha = -\cos (kL) - \sqrt{\cos^2(kL) - (1 - R)} \end{gather*} $ Where we take the negative root since we know that $\alpha < 0$. This gives us an explicit expression for $\alpha$ in terms of $R, k, L$. Our mixed boundary condition is then: $ \begin{gather*} c_1 u(\mathbf{x}) + c_2 \dfrac{\partial u}{\partial n} + c_3 = 0 \\ \Rightarrow c_1 u(\mathbf{x}) + \alpha c_1 = 0 \\ \Rightarrow u(L) = -\alpha \end{gather*} $ Finally, we obtain our boundary condition: $ u(L) = \cos (kL) + \sqrt{\cos^2(kL) - (1 - R)} $ Where the electric field at the boundary $z = L$ is then given by: $ \hat{\mathbf{n}} \times \mathbf{E} = u(L) $ For instance, for a $R = 95\%$ reflective mirror and using $\lambda = \pu{633 nm}$ light, substituting in $k = 2\pi/\lambda$ and $L = \pu{30 cm}$ gives us $|u(L)| = -0.044$. We note that in the special case of $R \to 1$ (that is, all light is reflected), we have: $ u(L) = \cos (kL) + \sqrt{\cos^2(kL)} = \cos(kL) - \cos (kL) = 0 $ This is the Dirichlet boundary condition that we would expect for a perfectly-reflective boundary in one dimension. Meanwhile, in the case of $R \to 0$ (that is, all light is transmitted), we have: $ u(L) = \cos (kL) + \underbrace{\sqrt{\cos^2(kL) - 1}}_{\text{only real for } \cos^2(kL) = 1} = \cos(kL) $ Which describes outgoing plane waves $u(x) = \cos(kx) = \operatorname{Re}(e^{ikx})$, as we would expect. This is a *very crude* method and a realistic simulation would most likely use a **partially-absorbent boundary condition** to model the fact that the mirror at $z = L$ absorbs *some* light. > **Note:** in addition to what has been mentioned, realistic laser cavities typically use parabolic/spherical mirrors instead of flat mirrors, so the domain will be curved at either end, making things much more complicated. Meanwhile, the sides of the laser cavity are usually glass or some other optically-transparent material, making it an open boundary (also called a *radiative boundary*). This can be modelled by a [perfectly-matched layer](https://en.wikipedia.org/wiki/Perfectly_matched_layer) that allows light to freely pass through. While counterintuitive, this is what maintains the strong directionality of the light along the transmission axis ($+z$) and keeps it as a beam rather than dispersing. However, a PML requires a bit of work to implement, so it is often easier to just use the first-order approximation of the *Sommerfeld radiative boundary condition*: $ \left(\dfrac{\partial \mathbf{E}}{\partial r} -ik\mathbf{E} = 0\right)_{|r| = R} $ Putting everything together, we arrive at the following boundary conditions for the electric field *inside* a laser cavity (although it can be generalized to the outside field as well): $ \begin{gather*} E_\perp(r,0) = 0, \quad E_\perp(r, L)= \eta E_0 \\ \dfrac{\partial E_\perp(z, R)}{\partial r} -ikE_\perp(z, R) = 0, \\ \dfrac{\partial E_\perp(z, -R)}{\partial r} -ikE_\perp(z, -R) = 0, \\ r \in [0, R], \quad z \in [0, L] \end{gather*} $ > **Note**: This means that the resultant electric field will be *complex-valued*; we take only the real part for a physically-meaningful solution. ## Conventional laser/maser modelling To model conventional lasers, which are by definition quantum devices (from _stimulated emission_ in the LASER acronym), requires a combination of classical and quantum models. Lasers can vary greatly in terms of pumping mechanism (or lack thereof), lasing wavelength, number of levels, etc. so it is not easy to come up with a model general enough to describe *all* (quantum) laser systems. Accordingly, to be able to create a basic model that can be numerically and analytically-solved, we must use some approximations. First, we assume that even in three-level or four-level systems, the transition that leads to stimulated emission is the *dominant contribution* to the laser's overall operation. This assumption can be justified (to some extent) by the fact that even three- and four-level systems generally have just one transition that produces stimulated emission, whereas the other transitions are in generally either non-radiative or produce spontaneous emission. Second, we assume that we can separately analyze the *macroscopic* and *microscopic* behavior of lasers. In more technical terms, this means that we can analyze lasers at a macroscopic level with *classical* electrodynamics (i.e. with electric and magnetic fields that obey the [Maxwell equation](https://en.wikipedia.org/wiki/Maxwell%27s_equations)) while analyzing the microscopic behavior of laser with *quantum electrodynamics*. This is justified by the fact that while lasers function because of *quantum* phenomena (involving individual atoms and photons), the light they produce at macroscopic scales (where we consider quadrillions of atoms and photons) is accurately described by the Maxwell equations. This also means that we can use the model we already described for modelling basic laser cavities to model the macroscopic-scale behavior of lasers, saving us a lot of work! Next, let's look at modelling the microscopic behavior of lasers. The most common quantum model of lasers is the [Jaynes-Cummings model](https://en.wikipedia.org/wiki/Jaynes%E2%80%93Cummings_model), which describes the interaction of an atom with a quantized electromagnetic field. It is based on the **Jaynes-Cummings Hamiltonian**, given by[^3]: $ \hat H = \hbar \omega_c \hat a^\dagger \hat a + \hbar \omega_a \dfrac{\hat \sigma_z}{2} + \dfrac{1}{2} \hbar \Omega \hat E \hat S $ There is a bit of explaining to do Here, $\omega_c$ is the oscillation frequency of the electromagnetic field (also called the *resonant frequency* of the system), and $\omega_a$ is the transition frequency, where $E_a = \hbar \omega_a$ is the difference in the upper and lower energy levels in the two-level system. $\omega_c$ is based on the geometric characteristics of the laser cavity, which restricts the electromagnetic field to only certain modes, each described by integer $n$. Considering only the **fundamental mode** ($n = 1$), a laser cavity of length $L$ would have $\omega_c = kc =\dfrac{\pi c}{L}$, where $c$ is the speed of light and $k$ is the wavenumber of the electromagnetic field, which is always $k = n\pi/L$ for the $n$-th mode. > **Note:** the Jaynes-Cummings model requires $\omega_a$ to be set as an input parameter, and it doesn't tell us what the value of $\omega_a$ should be. Instead, one usually obtains $\omega_a$ based on consulting a database (such as [NIST's atomic spectra database](https://www.nist.gov/pml/atomic-spectra-database)) for the type(s) of atom(s) that the laser's gain medium is filled with. Continuing on, $\Omega$ is the [Rabi frequency](https://en.wikipedia.org/wiki/Rabi_frequency) and is a function of both $\omega_c$ and $\omega_a$. Lastly, $\hat S$ is the spin operator, which is used as a simplified model of the atom (it doesn't mean the atom actually *has* spin, it is just a model that *approximates* the two dominant energy levels of the atom). Lastly, $\hat E$ is the quantized electric field operator, given by: $ \hat E = E_{ZPF}(\hat a + \hat a^\dagger), \quad E_{ZPF} = \sqrt{\frac{\hbar \omega_c}{2\varepsilon_0 V}} $ Where $\hat a, \hat a^\dagger$ are the [creation and annihilation operators](https://en.wikipedia.org/wiki/Creation_and_annihilation_operators) that model the emission and absorption of photons, $E_{ZPF}$ is the energy associated with [zero-point fluctuations](https://en.wikipedia.org/wiki/Zero-point_energy) of the electromagnetic field in the cavity, and $V$ is the volume of the cavity. $\hat \sigma_z$ is the Pauli $z$-operator, which corresponds to the matrix: $ \hat \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} $ Phew, that's a lot! Luckily, there exists free and open-source software like [QuTiP](https://qutip.org/) that can be used to solve the Jaynes-Cummings model, which (among other things) helps us understand the intensity of the radiation through time, as well as the population of the upper and lower levels. These are all crucial for understanding stimulated emission and the overall performance of the laser. > **Note:** An excellent online video series on quantum optics (including the Jaynes-Cummings model) from [LMU Munich University](https://www.lmu.de/en/) can be found [at this YouTube link](https://www.youtube.com/playlist?list=PL4_zMhS4uvR8dSbM-f_KDivMIAHi3guq2). ## Free-electron laser/maser modelling Free-electron lasers/masers are an exception to the standard methods of modelling lasers because unlike regular lasers, they operate via classical principles. More specifically, they use clever techniques to amplify [synchrotron radiation](https://en.wikipedia.org/wiki/Synchrotron_radiation) from electrons accelerated via a magnetic field and achieve coherence. The basics of free electron masers and a theoretical model has already been discussed in [[Free-electron maser physics]] and [[Numerical modelling of electron beams]]. Accordingly, we use the same coordinate conventions as in [[Free-electron maser physics]], where $x$ is the optical axis, $y$ is the direction parallel to the north-south axis of the magnets, and $z$ the direction perpendicular to the magnets. In that same article, we arrived at a model given by: $ \begin{align*} &\begin{cases} \mathbf{B} = \mathbf{B}_u(\mathbf{r}) + \mathbf{B}_q(\mathbf{r}, \mathbf{r}') + \mathbf{B}_r(\mathbf{r}) \\ \mathbf{E} = \mathbf{E}_q(\mathbf{r}, \mathbf{r}') + \mathbf{E}_r(\mathbf{r}) \\ \dfrac{d}{dt}(\gamma m\mathbf{v}') = q(\mathbf{E} + \mathbf{v}' \times \mathbf{B}) \end{cases} \\ &\begin{cases} \nabla^2 \mathbf{B}_r = -k^2 \mathbf{B} \\ \nabla^2 \mathbf{E}_r = -k^2 \mathbf{E} \end{cases} \\ &\begin{cases} \nabla \cdot (\mu_0(-\nabla \psi_m) + \mathbf{M}) = 0 \\ \mathbf{B}_u = \mu_0\mathbf{H} = -\mu_0(-\nabla\psi_m + \mathbf{M}) \end{cases} \\ \end{align*} $ Where $\mathbf{v}' = \dot{\mathbf{r}}'$ is the velocity of the electrons, $\mathbf{B}_u$ is the magnetostatic field generated by the magnets in the undulator cavity, $\mathbf{E}_q, \mathbf{B}_q$ are the fields of the moving electrons, and $\mathbf{E}_r, \mathbf{B}_r$ are the radiation fields of standing waves in the cavity. This is quite a complicated model so let's first try to solve it analytically with suitable approximations. First, we use the result from [[Free-electron maser physics]], where we derived that $\mathbf{B}_u$ takes the form: $ B_{ux} = B_{uz} = 0, \quad B_{uy} = B_0 \sin k_0 x $ Thus we have $\mathbf{B}_u = B_0 \hat{\mathbf{y}} \sin k_0 x$, where $k_0 = 2\pi/L_u$ and $L_u$ is the spatial period of the magnets. We will also use the approximate solution of the Lorentz force equation from [[Numerical modelling of electron beams]], which tells us that: $ x(t) = v_0t, \quad z(t) = \dfrac{qB_0}{m \omega_0 k_0} \sin \omega_0 t, \quad \omega_0 = k_0 v_0 $ Where we set $q = -e$ for electrons. This solution is true assuming that the electron is not moving beyond $\approx 0.5 c$, that is, moving at slower than 50% of the speed of light, and that the electric field doesn't affect the electrons. The first assumption can be kept for the purposes of our simplified solution, but the second must be altered. This is because the electric field in an undulator *does* affect the electrons. As the synchrotron radiation of the electrons is reflected by the walls of the cavity, it forms standing waves, which corresponds to an electric field given by: $ \mathbf{E}_r(x,t) = E_0 \hat{\mathbf{x}} \sin \dfrac{n\pi x}{L} \cos \Omega t, \quad \Omega = \dfrac{n\pi v}{L} $ Where $\Omega = kv$ is the frequency of the waves' oscillations in time, and $k = n\pi/L$ comes from the boundary conditions of the problem. Now, the electrons do not move with uniform velocity, but rather their velocities follow a roughly Gaussian distribution. This means the probability distribution of electron velocities would follow: $ P(v) \approx \dfrac{1}{\sqrt{\pi}} e^{-(v-v_0)^2} $ Where here, $P(v)$ is the probability (density) for an electron to have velocity $v$ (we assume no upper bound for the max velocity of electrons here, other than $v>0$). This means that the electric field in the cavity is, in general, a sum of all possible $v$: $ \mathbf{E}_r(x,t) \sim \int_0^{v_{max}} dv~ E_0 \hat{\mathbf{x}} \sin \left(\dfrac{n\pi x}{L}\right) \cos \left(\dfrac{n\pi v}{L}t\right) $ This means that while the electric field is spatially monochromatic, it is not monochromatic in time, as the electrons emit radiation of different temporal frequencies $\Omega$ depending on their velocity, leading to interference. However, this changes when we also consider the effects of the standing waves on the electrons, which comes from the electromagnetic field confined within the optical (resonant) cavity. Since the Lorentz force is linear in the non-relativistic regime, we can simply add a linear correction $\mathbf{r}_c(t)$ to compensate for the effects of the standing waves (we ignore the time component for now, since the spatial variation of the electric field is dominant in determining the electron trajectories except for very, very long distances). By substitution of the Lorentz force for an electric field, we have: $ \mathbf{F} = m\dfrac{d^2 \mathbf{r}_c}{dt^2} = q\mathbf{E} = -e E_0 \hat{\mathbf{x}} \sin \dfrac{n\pi x}{L} $ The electric field is only along $x$, so this simplifies to: $ \frac{d^2x_c(t)}{dt} = -\dfrac{e}{m} E_0 \sin\dfrac{n\pi x}{L} $ This differential equation is nonlinear, so we must linearize it with $\sin \phi \approx \phi$. Thus we have: $ \frac{d^2x_c(t)}{dt} \approx -\dfrac{e}{m} E_0 \dfrac{n\pi}{L}x $ Which gives us the solution: $ x_c(t) = -\dfrac{eE_0}{mL}n\pi \sin(\omega_c t), \quad \omega_c = \sqrt{\dfrac{n\pi E_0 e}{mL}} $ The complete solution $x(t)$ is obtained by adding our correction $x_c(t)$. Assuming that the fundamental mode ($n = 1$) is dominant, we have: $ x(t) = v_0t - \dfrac{e\pi E_0}{mL} \sin(\omega_c t), \quad z(t) = -\dfrac{eB_0}{m \omega_0 k_0} \sin \omega_0 t $ Notice how the electrons are no longer moving at a constant speed along $x$, but rather speed up and slow down in a harmonic fashion. This is the phenomenon of **microbunching**, where the electrons longitudinally move forwards and backwards in "bunches" based on their oscillation frequency $\omega_c$, which is restricted to specific values of $n$. This means that each bunch of electrons is **in phase** with one another, leading to constructive interference between the radiation they emit, so the electric field is no longer an integral over $v$; rather, it becomes just a sum: $ \mathbf{E}_r(x,t) = E_0 \hat{\mathbf{x}} \sum_n\sin \left(\dfrac{n\pi x}{L}\right) \cos \left(\omega_c t\right), \quad \sqrt{\dfrac{n\pi E_0 e}{mL}} $ This is a generalized result (Fourier series) for a closed-off undulator with perfectly-reflecting mirrors on both sides in one dimension. Applying the specific boundary conditions we discussed at the start of the guide - as well as applying the paraxial approximation - gives us the Gaussian beam as the particular solution. ### Achieving resonance in free electron lasers For a free-electron laser, a key factor in its performance is how well its radiative frequency $f_{r} = c / \lambda_r$ matches with the resonant frequency $f_{c}$ of the laser cavity. ## Numerical simulations for free electron lasers/masers ### General foundations A free-electron laser (or maser) is quite complicated to simulate, as it is a multiphysics problem: interactions between at least four (if not more) dynamical systems must be considered: - The magnetostatic field generated by the magnets inside the undulator - The electron beam and its trajectory in the static magnetic field of the undulator magnets - The electromagnetic radiation ($\mathbf{E}$ field) caused by the transverse (synchrotron) and longitudinal (pondermotive) acceleration of the electrons - The electromagnetic waves propagating in the laser/maser cavity, which requires solving the (paraxial) Helmholtz equation for the electric/magnetic fields as an eigenvalue problem to obtain the modal decomposition of the field within the cavity While theoretical models (like the one outlined in [[Free-electron maser physics]]) can be solved exactly in simple, idealized cases, they cannot be solved exactly for most cases. It is then necessary (and almost *always* necessary in general) to use a numerical solver. For this, it is possible to write custom code in association with general PDE/ODE solvers like FreeFEM (for the magnetic fields) and SciPy (for numerical integration of beam trajectories). ### Software implementations However, there exist dedicated software specifically designed for computer simulations of particle beams and undulators, such as the open-source [OSCARS software](https://oscars.bnl.gov/index.php) from Brookhaven National Laboratory with extensive [examples](https://oscars.bnl.gov/examples.php) and support for importing CAD models (a must-have for realistic free electron laser/maser designs). It is recommended to use these when possible, instead of custom-written solvers, for improved performance and accuracy. For visualization, an exported `.vtk` file is very useful, as it can be exported to professional visualization software like Paraview. ## Reference table The fundamental parameters of an undulator are: | Name | Symbol | Units | | --------------------------------- | -------------- | ------------------------------------------ | | Electron energy | $E_e$ | keV, MeV, GeV | | Electron beam voltage | $\Delta V$ | kV, MV, GV (for high-energy installations) | | Undulator magnetic field strength | $B_0$ | T, Gauss | | Undulator (spatial) period | $L_u$ | cm | | Number of magnet pairs along axis | $N_\text{mag}$ | dimensionless | > Note that some of these are interrelated. For instance, $E_e = e\Delta V$ (where $e$ is the elementary charge). [^1]: From the Wikipedia article on the [paraxial Helmholtz equation](https://en.wikipedia.org/wiki/Helmholtz_equation#Paraxial_approximation) [^2]: See [Section 2.4, _Paraxial Wave Equation and Gaussian Beams_](https://ocw.mit.edu/courses/6-974-fundamentals-of-photonics-quantum-electronics-spring-2006/871c32e6e4a44cbb1546741ef06f0f2f_parax_wav_eq_gau.pdf) from the MIT OpenCourseWare course [Fundamentals of Photonics: Quantum Electronics](https://ocw.mit.edu/courses/6-974-fundamentals-of-photonics-quantum-electronics-spring-2006/pages/lecture-notes/). Despite the name, note that the derivation is entirely classical. [^3]: See the Wikipedia article on the [Jaynes-Cummings model](https://en.wikipedia.org/wiki/Jaynes%E2%80%93Cummings_model) and also the article on [detuning](https://en.wikipedia.org/wiki/Laser_detuning). [^4]: See [this table](https://en.wikipedia.org/wiki/Electrical_resistivity_and_conductivity#Table). [^5]: From [this COMSOL guide](https://www.comsol.com/model/download/1405071/models.rf.rcs_sphere.pdf) for modelling radars. It refers to something known as a *perfect electrical conductor* boundary condition, which is described more in [this other COMSOL article](https://www.comsol.com/blogs/modeling-metallic-objects-in-wave-electromagnetics-problems/). [^6]: Inspired by [this Physics SE article](https://physics.stackexchange.com/questions/495354/boundary-condition-for-partial-reflection) [^7]: Also see [this Physics SE article](https://physics.stackexchange.com/a/597710) for a more rigorous semi-derivation