Guide to modelling lasers and masers

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,\varphi ,z)=E_0\frac{w_0}{w(z)}\exp \left(-\frac{r^2}{w(z{)}^2}\right)\exp \left(-i\left(kz+\frac{k{r}^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) 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:

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¹:

$$\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)u+2ik\frac{\partial u}{\partial z}=0$$

Note that if we define ${\nabla }_{(2D)}^2=\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 }_{(2D)}^{2}u+2ik\frac{\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,\varphi ,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\frac{\partial }{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\varphi }^2}+2ik\frac{\partial u}{\partial z}=0$$

Now, by assuming that the beam has azimuthal symmetry (that is, it doesn’t depend on the azimuthal angle $\varphi$), the derivative with respect to $\varphi$ drops out, giving us:

$$\begin{array}{c}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial }{\partial r}\right)+2ik\frac{\partial u}{\partial z}=0\\ \Rightarrow \frac{1}{r}\frac{\partial u}{\partial r}+\frac{{\partial }^{2}u}{\partial r^2}+2ik\frac{\partial u}{\partial z}=0\end{array}$$

To derive the Gaussian beam solution in particular, we assume that $u(\mathbf{r})$ takes the form of a Gaussian profile²:

$$u(\mathbf{r})=f(\mathbf{r})\exp (-\alpha r^2),\alpha =\frac{1}{w(z{)}^2}+\frac{ik}{2R(z)}$$

Where $w(z)$ describes the beam width, and $R(z)\sim z\left(1+\frac{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 of the first kind. The general solution reads as follows:

$$u(r,\varphi ,z)=\sum _m\sum _n{A}_{mn}{J}_n(kr)e^{im\varphi }{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 $\varphi$) and assume continuous coefficients (such that the sum becomes an integral), we have:

$$u(r,\varphi ,z)=\int dkA(r)J_n(kr)e^{im\varphi }{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,\varphi ,z)=e^{-k_{z}z}\int dkA(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,\varphi ,0)=C{e}^{-a{r}^2}$$

Where $C$ and $a$ are some undetermined constants. This gives us:

$$u(r,\varphi ,z)=C{e}^{-k_{z}z}\int kdk{e}^{-a{r}^2}{J}_0(kr)$$

Here, the extra factor of $k$ in $kdk$ comes from the fact that the integral $\int kdk{e}^{-\alpha r}{J}_0(kr)$ is an (inverse) Hankel transform. This particular transform has a well-known solution in the form $e^{-b{r}^2}$, where $b$ is a constant that is related to $a$. This gives us:

$$u(r,\varphi ,z)=C{e}^{-b{r}^2}{e}^{-k_{z}z}=C{e}^{-(b{r}^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.³

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, 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,\varphi )=E(r,\varphi +2\pi )$) with $E(r)\to 0$ as $r\to \infty$. Since the electric field and its derivative are non-vanishing on the partially-transmissive 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.

Modelling realistic laser resonant cavities

As we discussed prior, 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. The exact details for solving this boundary-value problem are given in RF cavity simulation with aperture.

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) 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, which describes the interaction of an atom with a quantized electromagnetic field. It is based on the Jaynes-Cummings Hamiltonian, given by⁴:

$$\hat{H}=\hslash {\omega }_c{\hat{a}}^{†}\hat{a}+\hslash {\omega }_a\frac{{\hat{\sigma }}_z}{2}+\frac{1}{2}\hslash \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=\hslash {\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=\frac{\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) for the type(s) of atom(s) that the laser’s gain medium is filled with.

Continuing on, $\Omega$ is the 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}}^{†}),E_{ZPF}=\sqrt{\frac{\hslash {\omega }_c}{2{\epsilon }_{0}V}}$$

Where $\hat{a},{\hat{a}}^{†}$ are the creation and annihilation operators that model the emission and absorption of photons, $E_{ZPF}$ is the energy associated with zero-point fluctuations 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=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)$$

Phew, that’s a lot! Luckily, there exists free and open-source software like QuTiP 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 can be found at this YouTube link.

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. (Note however that most conventional lasers are modelled using a mixture of classical and quantum methods). More specifically, they use clever techniques to amplify 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{array}{rl}& \{\begin{array}{l}\mathbf{B}={\mathbf{B}}_u(\mathbf{r})+{\mathbf{B}}_q(\mathbf{r},{\mathbf{r}}^{\prime })+{\mathbf{B}}_r(\mathbf{r})\\ \mathbf{E}={\mathbf{E}}_q(\mathbf{r},{\mathbf{r}}^{\prime })+{\mathbf{E}}_r(\mathbf{r})\\ \frac{d}{dt}(\gamma m{\mathbf{v}}^{\prime })=q(\mathbf{E}+{\mathbf{v}}^{\prime }\times \mathbf{B})\end{array}\\ & \{\begin{array}{l}{\nabla }^2{\mathbf{B}}_r=-k^2\mathbf{B}\\ {\nabla }^2{\mathbf{E}}_r=-k^2\mathbf{E}\end{array}\\ & \{\begin{array}{l}\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{array}\end{array}$$

Where ${\mathbf{v}}^{\prime }={\dot{\mathbf{r}}}^{\prime }$ 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,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_{0}t,z(t)=\frac{q{B}_0}{m{\omega }_0{k}_0}\sin {\omega }_{0}t,{\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.5c$, 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 \frac{n\pi x}{L}\cos \Omega t,\Omega =\frac{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 \frac{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(\frac{n\pi x}{L}\right)\cos \left(\frac{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\frac{d^2{\mathbf{r}}_c}{d{t}^2}=q\mathbf{E}=-e{E}_0\hat{\mathbf{x}}\sin \frac{n\pi x}{L}$$

The electric field is only along $x$, so this simplifies to:

$$\frac{d^2{x}_c(t)}{dt}=-\frac{e}{m}{E}_0\sin \frac{n\pi x}{L}$$

This differential equation is nonlinear, so we must linearize it with $\sin \varphi \approx \varphi$. Thus we have:

$$\frac{d^2{x}_c(t)}{dt}\approx -\frac{e}{m}{E}_0\frac{n\pi }{L}x$$

Which gives us the solution:

$$x_c(t)=-\frac{e{E}_0}{mL}n\pi \sin ({\omega }_{c}t),{\omega }_c=\sqrt{\frac{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_{0}t-\frac{e\pi E_0}{mL}\sin ({\omega }_{c}t),z(t)=-\frac{e{B}_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(\frac{n\pi x}{L}\right)\cos \left({\omega }_{c}t\right),\sqrt{\frac{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 from Brookhaven National Laboratory with extensive examples 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).

Footnotes

1. From the Wikipedia article on the paraxial Helmholtz equation ↩

2. See Section 2.4, Paraxial Wave Equation and Gaussian Beams from the MIT OpenCourseWare course Fundamentals of Photonics: Quantum Electronics. Despite the name, note that the derivation is entirely classical. ↩

3. Also see this Physics SE article for a more rigorous semi-derivation ↩

4. See the Wikipedia article on the Jaynes-Cummings model and also the article on detuning. ↩