RF cavity simulation with aperture

An RF cavity is the equivalent of a laser optical cavity, but hosting EM waves in the microwave range rather than visible light. It is an essential component of a maser, and can be solved both analytically and numerically, so we detail how to do so below.

The motivation of using RF cavities

The reason why we need to use a RF cavity in a free-electron maser (indeed, among masers in general) is that the cavity is a resonant cavity, meaning that it causes amplification of EM waves. The precise mechanism behind why this amplification occurs is hard to rigorously explain without solving Maxwell’s equations, but the general effect is that an RF cavity is the crucial component necessary to create a maser beam.

In a resonant cavity that is typical for a laser/maser, two mirrors are placed some distance apart, and joined with a tube made of a optically- or microwave-transparent material (depending on the operating wavelength of the laser/maser). This allows EM waves to be reflected between the metal plates, but not reflect off the walls of the tube, which creates a highly-directional and focused beam.¹ This page is not focused on the precise mechanics or engineering of a resonant (optical or RF) cavity; please see Guide to modelling lasers and masers for more information on the topic. Rather, it is focused on describing how to mathematically model and simulate the electromagnetic field inside and outside an RF cavity. Being able to do so will help us design better resonant cavities for high-performance masers.

A basic 2D model

Note: Our coordinate conventions are to use $z$ to denote the optical axis (direction of propagation of EM waves), and $r$ to denote the radial distance away from the optical axis, with 2D coordinate points labelled by $(r,z)$. When working in 3D, we use cylindrical coordinates $(r,\varphi ,z)$, where $r,z$ are still defined the same way, and $\varphi$ denotes the azimuthal (angular) coordinate. These are chosen for consistency with our maser coordinates, as well as the typical coordinates used for describing Gaussian beams.

We consider a RF cavity of total width $w$, aperture (hole) width $a_0$, and total length $L$, as shown in the below diagram:

Note: The above diagram is a 2D cross-section. The full 3D cavity (which has a cylindrical shape) is not shown.

To make things simpler, to start, we restrict our analysis to two dimensions only. The (scalar) Helmholtz equation for the $E_r$ component of the electric field $\mathbf{E}=E_r\hat{r}+E_z\hat{z}$ is given by:

$$({\nabla }^2+k^2)E_r=0,E_z=0$$

Note: The reason we can treat the electric field as a scalar problem is that we assume that there are no longitudinal modes of the electric field. That is to say, the electric field is always transverse (perpendicular) to the optical axis, or in mathematical terms, $E_{\parallel }=E_z=0$, leaving only the $E_r$ component left.

In our chosen coordinates and in two dimensions, the Helmholtz equation takes the form:

$$\frac{{\partial }^2{E}_r}{\partial z^2}+\frac{{\partial }^2{E}_r}{\partial r^2}+k^2{E}_r=0$$

Using analytical and numerical means, we’ll tackle this problem in detail.

Analytical solution

Note: While the overall solution is believed to be generally correct and mathematically sound, and the quantitative result of ${\lambda }_c=2L$ is well-known, the intermediary steps described here may not be fully correct. Researchers are encouraged to add citations to the literature (for instance, Pozar’s Microwave Engineering) and cross-check to verify that the derivation presented is fully correct.

In two dimensions, it is possible to solve the problem analytically in a fairly straightforward fashion using separation of variables, with the following assumptions:

• The left and right walls of the RF cavity (located at $z=0$ and $z=L$ respectively) are perfect conductors, meaning they reflect all incident EM waves and neither leak nor absorb energy
• The top and bottom walls of the RF cavity (located at $r=-w/2$ and $r=w/2$ respectively) are perfectly transmissive, meaning they allow EM waves to pass through freely
• The aperture is so small as to be negligible in size compared to the width of the RF cavity (that is, $a_0/w\ll 1$), so no energy leaks out the cavity

Note: While it may appear that energy would be needlessly dissipated over the top and bottom walls of the RF cavity, this actually does not happen, since the electromagnetic field quickly becomes zero for $r>R$.

With all of these idealized assumptions, together with the fact that we are ignoring the third dimension, it suffices to say that the results we will obtain will be highly approximate. However, it will prove to be useful before we tackle the much more complex problem of the 3D case, with a non-negligible aperture and non-trivial energy leakage and absorption.

Here, we will not go through the entire process of separation of variables, since it is quite lengthy and is more of an exercise in mathematics than physics or engineering. It suffices to know that after performing the separation of variables procedure, the general solution of the Helmholtz equation in our 2D Cartesian-like coordinates $(r,z)$ is given by²:

$$E_r(r,z)=\sum _{m=1}^{\infty }(A_m{e}^{mr}+B_m{e}^{-mr})[C_m\sin \left(\sqrt{m^2+k^2}z\right)+D_m\cos \left(\sqrt{m^2+k^2}z\right)]$$

We now impose boundary conditions to find a specific solution. We want our solution to apply for all $r$ and all $z\ge 0$, and therefore $A_m$ must be zero, or otherwise the solution blows up at $r\to \infty$. Meanwhile, we impose the requirement that $E_r$ must vanish at the surface of each of the reflective boundary walls (that is, at $z=0$ and $z=L$); this comes from the fact that mirrors are typically made of conductive materials (typically metals), and the tangential components of the electric field at the surface of any perfect conductor is zero. Here, $E_r$ is the tangential component, and therefore $E_r=0$ at both ends of the cavity. To satisfy the $z=0$ boundary condition, $D_m$ must be zero. Now, if we define $K_m=B_m{C}_m$, we have:

$$E_r(r,z)=\sum _{m=1}^{\infty }{K}_m{e}^{-mr}\sin \left(\sqrt{m^2+k^2}z\right)$$

Note: Here, $K_m$ describes the amplitude of the wave as it propagates in space, and depends on the value of $m$. The precise value of $K_m$ can be found through Fourier analysis, but we will not cover this for now.

Meanwhile, to satisfy the $z=L$ boundary condition, we must have $\sin (\sqrt{m^2+k^2}L)=0$. Solving this equation gives us $\sqrt{m^2+k^2}L=n\pi$, or alternatively:

$$k^2=\frac{n^2{\pi }^2}{L^2}-m^2,n=1,2,3,\dots$$

Thus, the modes that the cavity can support are directly related to the integer values of $m$ and $n$, so they are quantized. Now solving for $k$, we obtain:

$$k=\sqrt{\frac{n^2{\pi }^2}{L^2}-m^2},n=1,2,3,\dots$$

We observe that $k$ is only real-valued if $n^2{\pi }^2/L^2\ge m^2$. This means that the minimum value of $m$ is given by the case where $n=1$, that is:

$$m=\frac{\pi }{L}$$

From dimensional analysis we can tell that $m$ has the same units as $k$, and takes the interpretation of a cutoff, prohibiting modes above a certain wavelength³. Indeed, we may define a cutoff wavelength ${\lambda }_c=2\pi /m$, which is given by:

$${\lambda }_c=\frac{2\pi }{(\pi /L)}=2L$$

This is the longest wavelength that the RF cavity permits; it is impossible for any EM wave of a longer wavelength to be present within the cavity. Thus, the RF cavity restricts the possible modes of the electric (and magnetic) fields, allowing us to preferentially select for specific modes. Indeed, if an RF cavity is constructed very carefully to minimize attenuation/losses, then the fundamental mode dominates, and the RF cavity supports EM waves of (almost) purely one wavelength, that being $\lambda ={\lambda }_c$. We can convert this to frequency with $f=c/\lambda$, and thus $f_c=c/{\lambda }_c$. For instance, a $L=10\mathrm{cm}$ RF cavity would therefore have (in our 2D simplification) a fundamental mode with a wavelength of $\lambda =20\mathrm{cm}$, and thus a frequency of around $f=1.5\mathrm{GHz}$⁴.

Additional information: In purely 2D space, there are only TE (transverse electric) modes in the plane of the cavity; a transverse magnetic field must point out of the plane (into the 3rd dimension) and is thus out of our consideration.

A more advanced 2D model

In our basic model of a maser RF cavity, we assumed that the aperture was negligible in relation to the rest of the cavity. This, of course, is an over-simplification. We will now consider what happens if we do consider the effects of the aperture.

The physical effect of the aperture is that it allows some portion of the electromagnetic field inside the RF cavity to pass outside the cavity. This is a complex process that is not easily modelled analytically. However, we can model it with the following boundary condition (which we may term a transmissive boundary condition):

$$E_r(r,L)=T{E}_r(r,0)$$

Where $T$ is the transmission coefficient, and is a dimensionless number between zero and one that describes what percent of the energy (or more precisely, the power density) of the electromagnetic field is allowed to pass outside the RF cavity⁵. We will not care right now about how exactly one can derive $T$ from first principles, just that it is some known value, and is expressed as the ratio between the electric field’s amplitude at $z=L$ and at $z=0$. That is:

$$T=\frac{|E_r(r,L)|}{|E_r(0,L)|}$$

Where we require that $T$ be real-valued, and thus we must take the absolute value, since $E_r$ can be complex-valued. This definition is useless for actually calculating what $T$ is from first principles, since it requires knowing $E_r$ at the left and right ends of the cavity, which (by definition) we don’t know. However, It can be experimentally-determined by performing actual measurements, or numerically-determined by performing a curve fit on the results of finite-element simulations (more on that later)⁶.

In addition, we can recognize some basic qualitative features of $T$. In the limit as $a_0\to 0$, then $T\to 0$, meaning that 0% of the microwaves are transmitted (equivalently, 100% of the microwaves are reflected). This would make for excellent amplification, but would also make it impossible to extract any energy out of the interior of the cavity. Meanwhile, in the limit as $a_0\to w$, meaning that 100% of the microwaves are transmitted (equivalently, none of the microwaves are reflected). This would make the RF cavity equally useless, since all the energy is allowed to pass out of the cavity, and therefore there is no amplification that can happen. But if $T$ is somewhere in between - for instance, if $T$ is between 5-10% - it can allow for the RF cavity to still effectively amplify the EM field within, while allowing a reasonable amount of power out of the cavity. We will now solve the problem analytically to show this effect, and solve for the wavelengths (and frequencies) permitted by our waveguide.

Analytical solution

With our new transmissive boundary condition, let us start at our general solution from before, which was given by:

$$E_r(r,z)=\sum _{m=1}^{\infty }(A_m{e}^{mr}+B_m{e}^{-mr})[C_m\sin \left(\sqrt{m^2+k^2}z\right)+D_m\cos \left(\sqrt{m^2+k^2}z\right)]$$

Once again, we want the solution to be value for all $r$, and therefore $A_m$ must be zero for the solution to not blow up as $r\to \infty$. This gives us:

$$E_r(r,z)=\sum _{m=1}^{\infty }{B}_m{e}^{-mr}[C_m\sin \left(\sqrt{m^2+k^2}z\right)+D_m\cos \left(\sqrt{m^2+k^2}z\right)]$$

This time, we don’t use the perfectly reflective boundary conditions on the left and right wall of the cavity (i.e. setting $E_r\to 0$ at $z=0$ and $z=L$). Instead, we apply our new transmissive boundary condition, giving us (after dropping out common factors):

$$\begin{array}{rl}{C}_m\sin \left(\sqrt{m^2+k^2}L\right)& +D_m\cos \left(\sqrt{m^2+k^2}L\right)\\ & =T\cancel{C_m\sin \left(\sqrt{m^2+k^2}\cdot 0\right)}+T{D}_m\underset{1}{\underbrace{\cos \left(\sqrt{m^2+k^2}\cdot 0\right)}}\end{array}$$

After simplifying, we obtain:

$$C_m\sin \left(\sqrt{m^2+k^2}L\right)+D_m\cos \left(\sqrt{m^2+k^2}L\right)=T{D}_m$$

We can remove the (direct) dependency on $C_m$ and $D_m$ by defining a new constant $\rho \equiv \frac{C_m}{D_m}$. Thus, after dividing all sides by $D_m$, we have:

$$\rho \sin \left(\sqrt{m^2+k^2}L\right)+\cos \left(\sqrt{m^2+k^2}L\right)=T$$

Now solving for $\rho$, we have:

$$\rho =\frac{T-\cos (\sqrt{m^2+k^2}L)}{\sin (\sqrt{m^2+k^{2}L})}$$

We may rewrite $\rho$ in simpler fashion by defining a new dimensionless quantity $\phi$, such that:

$$\rho (\phi )=\frac{T-\cos (\phi )}{\sin (\phi )},\phi =\sqrt{m^2+k^2}L$$

What we see is that $\rho$ is negative until $\phi =\pi /2$, after which the sign is switched and we have positive $\rho$. Meanwhile, the first root of $\rho$ is present at $\phi ={\cos }^{-1}$. In addition, at $\phi =n\pi$, the value of $\rho$ asymptotes and therefore $E_r$ is ill-defined. This tells us that it is impossible for $\sqrt{m^2+k^2}L=n\pi$. Finally, if we still assume that the left wall of the cavity (not the right wall!) is sufficiently reflective enough to be near-perfectly-reflective, then it must be the case that at $z=0$ we have $E_r(r,0)=0$, meaning that $D_m=0$ (or otherwise it would be impossible for this to be the case).

All of these conditions together allow us to find the possible modes. By the requirement that $D_m=0$, it necessarily require that $\rho =0$, since $\rho =\frac{C_m}{D_m}$, and thus $\rho D_m=C_m$. This reduces the problem into finding the roots of $\rho (\phi )$ such that for these values of $\phi$, we have $\rho =0$ - equivalently, to find the roots ${\phi }_n^{\ast }$ of $\rho (\phi )$. Now recalling that $\phi \equiv \sqrt{m^2+k^2}L$, we can rearrange to find that:

$$k^2=\frac{({\phi }_n^{\ast }{)}^2}{L^2}-m^2,n=1,2,3,\dots$$

Or equivalently, we have:

$$k=\sqrt{{\left(\frac{{\phi }_n^{\ast }}{L}\right)}^2-m^2},n=1,2,3,\dots$$

Where ${\phi }_n^{\ast }$ (that is the roots of $\rho (\phi )$) can be analytically solved for by solving the transcendental equation $\rho ({\phi }_n^{\ast })=0$, and are given by:

$${\phi }_n^{\ast }=\{\begin{array}{ll}n\pi -{\cos }^{-1}(T),& n=2,4,6,8,\dots \\ (n-1)\pi +{\cos }^{-1}(T),& n=1,3,5,7,\dots \end{array}$$

Now again by the requirement that the square root be real, it must be the case that:

$$m^2\le {\left(\frac{{\phi }_n^{\ast }}{L}\right)}^2,{\lambda }_n=\frac{2\pi }{m}=\frac{2\pi L}{{\phi }_n^{\ast }}$$

Therefore, the fundamental wavelength ${\lambda }_c=2\pi /m$ is given by:

$${\lambda }_c=\frac{2\pi L}{{\phi }_1^{\ast }}$$

For $T\approx 0.1$ (that is, if 10% of the microwaves were allowed to pass out of the cavity) and $L=10\mathrm{cm}$, the first few roots of $\rho (\phi )$ are numerically given by the following table:

Mode	Root (${\phi }_n^{\ast }$)	Fundamental wavelength	Frequency
1st (fundamental)	1.47063	42.7 cm	701.7 MHz
2nd	4.81256	13.1 cm	2.3 GHz
3rd	7.75381	8.1 cm	3.7 GHz

Finally, as a consistency check, let us note that occurs in the limit $T\to 0$ (where the aperture disappears and both the left and right walls of the cavity become perfectly reflective). In this case, we have:

$$\underset{T\to 0}{\lim }\rho (\phi )=\frac{\sin (\phi )}{\cos (\phi )}=\tan \phi$$

Now setting the requirement that $D_m\to 0$ for both the left and right walls of the cavity (that is, at $z=0$ and $z=L$), as with before, this requires that $\rho =\tan \phi =0$. The roots of $\tan \phi$ are $\phi =n\pi$, and thus we have:

$$k=\sqrt{\frac{({\phi }_n^{\ast }{)}^2}{L^2}-m^2}=\sqrt{\frac{(n\pi {)}^2}{L^2}-m^2},n=1,2,3,\dots$$

Therefore, upon taking the requirement that the square root be real-valued, we have:

$$m=\frac{n\pi }{L}$$

Upon which taking the fundamental ($n=1$) mode gives us:

$${\lambda }_c=\frac{2\pi }{(\pi /L)}=2L$$

This is exactly the result we obtained previously, and shows that our solution is indeed consistent as it reproduces the result of our first derivation in the perfectly reflective, zero-aperture limit.

Numerical solution

To solve the Helmholtz equation numerically, we use the finite element method with the open-source FreeFEM software. This is essential to model physically-realistic cavities that have energy losses. However, in this case, we will not consider such losses in a precise fashion, and find a numerical solution purely based off our simpler model, where we abstract away the precise physics based on a generic transmission coefficient $T$. We interpret $T$ as combining the two sources of power losses:

1. Transmission due to the aperture, which is designed to let out a portion of the microwaves
2. Transmission due to power leakage from the cavity, which is undesirable but nonetheless present in any realistic RF cavity

We will also leave the 3D cylindrical case and other more complex effects for later, and consider the basic 2D Cartesian case only. While this is certainly a very crude model, this numerical simulation is intended to be for familiarization with the concepts of the finite-element method; we will eventually discuss how to perform a numerical simulation that is physically-meaningful and can yield realistic results.

A crash course in the finite element method

To speak in broad terms, the finite element method approximates a solution to a PDE as a sum of lots of simple piecewise functions called basis functions⁷. If we let ${\psi }_i(\mathbf{x})$ be the $i$-th basis function, and $c_i$ be a constant coefficient to multiply the basis function by, then the solution $E(\mathbf{x})$ can be approximated using the sum of basis functions as:

$$E(\mathbf{x})\approx \sum _i{c}_i{\psi }_i(\mathbf{x})$$

The goal of the finite element is to find the unknown coefficients $c_i$ that best approximate the true solution $E(\mathbf{x})$ of our PDE (the Helmholtz equation). The first step in doing so is to convert the PDE into the variational form (also known as the “weak form”) - essentially, an integral equation that contains the same information as the PDE. We’ll do this in several steps. The first step is to write the PDE in standard form, meaning that we need to move all terms to the left-hand side such that the right-hand side is equal to zero. Luckily, the Helmholtz equation is already in this form:

$${\nabla }^{2}E+k^{2}E=0$$

The next step is to multiply every term of the Helmholtz equation with a test function $\phi$, and integrate the entire left-hand side. A test function is a simple function that has a well-known form, such as a polynomial. This results in:

$${\int }_{\Omega }\phi {\nabla }^{2}E+{\int }_{\Omega }{k}^2\phi E=0$$

Where $\Omega$ denotes the domain in which we’re solving the Helmholtz equation over. Now, the crucial part of the finite element method is that we have to write this integral equation as two terms: one integrated over the domain $\Omega$, and one integrated over the boundary of the domain, which we denote as $\partial \Omega$. To do this, we use integration by parts. Recall that the integration by parts formula in one variable says that:

$${\int }_a^{b}udv=uv{|}_a^b-{\int }_a^{b}vdu$$

The multivariable version of integration by parts takes the form:

$${\int }_{\Omega }udv={\int }_{\partial \Omega }uv-{\int }_{\Omega }vdu$$

Here, we apply integration by parts to our term ${\int }_{\Omega }\phi {\nabla }^{2}E$. Let $u=\phi$ and $dv={\nabla }^{2}E$. Then, we have $du=\nabla \phi$ and $v={\int }_{\Omega }{\nabla }^{2}E=\nabla E$. So using the formula, we get:

$${\int }_{\Omega }\phi {\nabla }^{2}E={\int }_{\partial \Omega }\phi \nabla E-{\int }_{\Omega }\nabla \phi \cdot \nabla E$$

Substituting this result back into the integral equation we got previously, we have:

$${\int }_{\partial \Omega }\phi \nabla E-{\int }_{\Omega }\nabla \phi \cdot \nabla E+{\int }_{\Omega }{k}^2\phi E=0$$

We can combine the two integrals over $\Omega$ and multiply by $-1$ so that we have the simplified result in a specific order, where bilinear (that is, nonlinear) terms come first, and linear terms come after:

$${\int }_{\Omega }\nabla \phi \cdot \nabla E-{\int }_{\Omega }{k}^2\phi E-{\int }_{\partial \Omega }\phi \nabla E=0$$

This is the most general variational form of the Helmholtz equation. The boundary term (the integral over $\partial \Omega$) is the one we’re most interested in, because it depends on our choice of boundary conditions (more on that in the next section). Note that in the finite element world it is common to use the following condensed notation⁸:

$$(\nabla \phi ,\nabla E)-(k^2\phi ,E)-⟨\phi ,\hat{n}\cdot \nabla E⟩=0$$

This comes from the mathematical notation for functional analysis, where the inner product of two functions $f(\mathbf{x}),g(\mathbf{x})$ is denoted $(f,g)$ or $⟨f,g⟩$, and expands to the following integrals:

$$(f,g)={\int }_{\Omega }f(\mathbf{x})g(\mathbf{x}),⟨f,g⟩={\int }_{\partial \Omega }f(\mathbf{x})g(\mathbf{x})\cdot \hat{n}$$

The reason we need to do all this work to get the variational form is that finite element software (like FreeFEM, which we use) can convert this to a matrix equation to solve for $c_i$ using the methods of numerical linear algebra. And since integrals are well-defined even over discontinuous or non-smooth domains (unlike partial derivatives), the variational form used by the finite element method means that you can solve PDEs on very complex geometries, including those with sharp corners, holes, and edges, which are commonly encountered in real-life objects; in such cases, the partial derivatives in the so-called strong form (original differential form) are ill-defined.

Setting up boundary conditions

In our simplified model, we aim to solve for the region $z\in [0,L]$ and $r\in [-R,R]$, where $R>w/2$. To start, we have the transmissive boundary condition from before, where $T$ is a constant:

$$E_r(r,L)=T{E}_r(r,0)$$

In practice, it is easier to set $E_r(r,0)$ to a fixed value and then set $E_r$ based on the transmissive boundary condition. We can set an order-of-magnitude estimate by assuming that the interior electric field is sourced from some spatially and temporally-varying potential source $V(z,t)=V_0{e}^{ikz-i\omega t}$, where $\omega =kc=2\pi c/\lambda$. Therefore, we have:

$$\mathbf{E}=-\nabla V=-\frac{\partial V}{\partial z}=ik{V}_0{e}^{ikz-i\omega t}=V_{0}k{e}^{ikz-i\omega t+i\pi /2}$$

Therefore, the amplitude of the electric field $E_0$ at $z=0$ is given by:

$$E_0=|\mathbf{E}|=V_{0}k=\frac{2\pi V_0}{\lambda }=\frac{2\pi V_{0}f}{c}$$

Making a basic order-of-magnitude approximation for the frequency/wavelength (which can always be refined later), we arrive at a figure for about $42\mathrm{V}/\mathrm{m}$ for $f=1\mathrm{GHz}$ and a $2\mathrm{V}$ power source.

Having placed boundary conditions for the left and right walls, we now need to consider the top and bottom sides of the RF cavity. This is a more complex situation since these regions are open boundaries that allow microwaves to pass unimpeded and radiate away to $r\to \pm \infty$; it is in fact this particular property that allows our RF cavity to be any use at all for use in a maser.

To accurately model outgoing waves that radiate to infinity, we must use the Sommerfeld radiation condition, which, in $d$ dimensions, takes the form:

$$\underset{r\to \infty }{\lim }{r}^{\frac{d-1}{2}}\left(\frac{\partial }{\partial |\mathbf{r}|}-ik\right)E(\mathbf{r})=0$$

Where $r\equiv |\mathbf{r}|$ is the magnitude of the radial vector pointing from the origin of the coordinate system (and is always positive). For instance, the two- and three-dimensional variants of the Sommerfeld radiation condition take the form:

$$\begin{array}{rlr}\underset{r\to \infty }{\lim }\sqrt{r}\left(\frac{\partial }{\partial |\mathbf{r}|}-ik\right)E(\mathbf{r})=0,& & \text{(2D)}\\ \underset{r\to \infty }{\lim }r\left(\frac{\partial }{\partial |\mathbf{r}|}-ik\right)E(\mathbf{r})=0,& & \text{(3D)}\end{array}$$

The requirement for this to be exact, of course, is that the domain is infinitely-large. Of course, we know that simulating infinite domains computationally is impractical. However, we can approximate this boundary condition by surrounding our RF cavity in a large circular region of radius $R$ (shown below), where $R\gg L$ (the circular region is much larger than the RF cavity). Thus, we can say that $R\approx \infty$; it may seem preposterous, but our assumption that $R\gg L$ means that the circular domain can be effectively treated as infinite.

Note: the same general idea can easily be extended into 3D: the only difference is that the circular boundary will need to be changed into a spherical boundary.

We note that since only the top and bottom sides of the RF cavity permit the unimpeded passage of outgoing waves, that is, along the $\pm r$ axes, we must modify the 2D condition to selectively choose the $r$ axis, rather than the general distance $|\mathbf{r}|$. Thus, we have:

$$\underset{r\to \infty }{\lim }\sqrt{r}\left(\frac{\partial }{\partial r}-ik\right)E(r,z)=0$$

If the domain is very large, it is possible to approximate the radiative boundary condition with simpler periodic boundary conditions. See the visualization here for a numerical simulation using periodic boundary conditions (this simulation assumes $T=0$). Indeed, one can approximate the 2D Sommerfeld radiation condition up to first-order with⁹:

$$\frac{\partial E}{\partial r}-ikE=0$$

The solutions to this differential equation over slices of constant $z$ are sinusoids, by inspection; thus this approximation results in plane-waves propagating far away from the sides of the cavity. We will now numerically demonstrate how to compute the electric field. First, recall that our weak form was given by:

$${\int }_{\Omega }\nabla \phi \cdot \nabla E-{\int }_{\Omega }{k}^2\phi E-{\int }_{\partial \Omega }\phi \nabla E=0$$

We now apply our boundary conditions on the third term, which describes an integral performed over the boundaries of the domain the problem is to be solved over. We incorporated the transmissive boundary conditions by setting $E_r(r,0)$ and $E_r(r,L)$ to fixed values with $E_r(r,L)=T{E}_r(r,0)$ and $|E_r(r,0)|=E_0$. These do not affect the boundary term since by nature of being constant values (formally, by being a Dirichlet boundary condition) the gradient of the electric field is zero. However, our open boundaries at $r=-w/2$ and $r=w/2$ do factor in, since they specify a non-vanishing derivative of the electric field at the aforementioned boundaries. Note that since we are considering a 2D domain $\Omega$, the boundary of our 2D domain is a 1D curve, so integrals over $\Omega$ are surface integrals while integrals over $\partial \Omega$ are line integrals. Thus adding the appropriate integral notations, we have:

$$\begin{array}{rl}{\oint }_{\partial \Omega }(\phi \nabla E\cdot \hat{n})d\mathbf{r}& =\underset{r=-w/2}{\oint }\phi \frac{\partial E}{\partial r}d\mathbf{r}+\underset{r=w/2}{\oint }\phi \frac{\partial E}{\partial r}d\mathbf{r}\\ & =\underset{r=-w/2}{\oint }\phi E\left(-\frac{w}{2},z\right)d\mathbf{r}+\underset{r=w/2}{\oint }\phi E\left(\frac{w}{2},z\right)d\mathbf{r}\end{array}$$

Therefore, the weak form of the Helmholtz equation for our chosen boundary conditions (and with the integrals fully written-out) is thus given by:

$${\iint }_{\Omega }(\nabla \phi \cdot \nabla E)dA-{\iint }_{\Omega }{k}^2\phi EdA-\underset{r=-w/2}{\oint }\phi E\left(-\frac{w}{2},z\right)d\mathbf{r}-\underset{r=w/2}{\oint }\phi E\left(\frac{w}{2},z\right)d\mathbf{r}=0$$

Where here, we simulate only the interior field (that is, the field inside the RF cavity itself); we may also simulate both the interior and exterior fields with the circular domain method described previously, but we will omit it for simplicity.

Simulation results

The FreeFEM code for the simulation can be found in freefem/dev/lasers/solve-Efield.edp. The result is shown below:

We see some notable features from the simulation. First, the radiation pattern rapidly decays away to infinity and doesn’t stay collimated. This suggests that it would be better to add waveguide to the aperture instead of a bare opening as the output coupler of a maser. Second, as we would expect, the field is zero almost everywhere but the interior of the RF cavity and right in front and around the aperture. Indeed, this is not very different from visible light in a dark room with a small hole in it (see camera obscura).

Note: The code does have some issues right now, including the fact that it is necessary to hard-code a value of the electric field at the aperture (here it is set to $E=E_0=1$) when in theory there should be no physical requirement for this.

The simulation also gives reference values of the electric field strengths for us to compare experimental results against, and thus improve our models. A physically-sound but impractical idea is to put a square plate in front of the maser, with lots of tiny wires crisscrossing the plate. Then, we can measure the voltage across each of the wires at regular intervals to get an approximate idea of the voltage across the entire plate (or we can interpolate numerically). Once we know the voltage $V(x,y)$ we can then figure out the power density (intensity) of the EM waves and compare this against the simulation. The much more practical method is to place a sensitive antenna (or series of antennas) in front of the maser. The antenna(s) can then measure the microwave beam from the maser, much like a radio telescope, allowing us to have an accurate “picture” of how the beam spreads.

A realistic 3D model

Note: for the analytical solution, reference Pozar, _Excitations of Waveguides - Aperture Coupling (Ch. 4.8). The numerical solution will use the actual geometry with a PML. See https://pyoomph.readthedocs.io/en/latest/tutorial/spatial/helmholtz.html

Having considered a simplified 2D model, we now turn our attention to modelling a realistic cylindrical waveguide in 3D. This problem is quite complex and some aspects of it can only be solved numerically, but we will find that analytical methods can get us a long way.

Before we start, we should give some motivation for why solving this problem is necessary at all. Why, for instance, do we need to solve for a cylindrical RF cavity, and not a rectangular one? The reason is that if we want to form a microwave beam by carefully letting power out of an RF cavity, we must have azimuthal symmetry. Azimuthal symmetry means that cross-sections of the beam are always symmetric about the azimuthal angle $\varphi$; otherwise, the “beam” would not be a beam, so to speak, but a propagating plane wave (with cross-sections of squares or rectangles). Thus, cylindrical RF cavities are essential for making functional masers, which is our main motivation for specifically analyzing a cylindrical RF cavity.

Second, why doesn’t a 2D approximation work? The reason is that the Helmholtz equation contains a Laplacian term (that is, ${\nabla }^2$); the Laplacian takes different forms in different coordinate systems, and thus solutions to the Helmholtz equation in different coordinates are distinct. While our 2D solution might suffice as a rough approximation to a 2D cross-section of the 3D fields, it is a poor approximation in general. This means that we have to start from scratch instead of tweaking our 2D solution.

Analytical solution

We will tackle the problem in several steps - analytically to begin with, numerically afterwards:

1. Solve the basic problem analytically for perfectly-reflecting walls and no power losses (that is, with no aperture and no power leakage)
2. Apply a basic transmissive boundary condition, that is, $E_r(r,\varphi ,L)=T{E}_r(r,\varphi ,0)$
3. Perturbatively incorporate the effects of lossy materials
4. Find a numerical solution using finite elements

We will solve for just the TE modes, since microwave beams are approximated by Gaussian beams, which are written in terms of the electric field, not the magnetic field.

To solve, we will first expand out the vector-valued Helmholtz equation $({\nabla }^2+k^2)\mathbf{E}=0$ for each of the nonzero components of the electric field. Since electric fields are always transverse to the direction of wave propagation, there is no $E_z$ field, so we have:

$$\begin{array}{r}{\nabla }^2{E}_r+k^2{E}_r=0\\ {\nabla }^2{E}_{\varphi }+k^2{E}_{\varphi }=0\end{array}$$

Since we are solving in a problem with azimuthal symmetry, we will write out the Laplacian in cylindrical coordinates, where it reads:

$$\begin{array}{rl}{\nabla }^2& =\frac{1}{r}\left(r\frac{\partial }{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\varphi }^2}+\frac{{\partial }^2}{\partial z^2}\\ & =\frac{1}{r}\left(\frac{\partial }{\partial r}+r\frac{{\partial }^2}{\partial r^2}\right)+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\varphi }^2}+\frac{{\partial }^2}{\partial z^2}\end{array}$$

We will now perform the classic separation of variables procedure. Let $E_j(r,\varphi ,z)=u_j(r,\varphi )Z_j(z)$. We may therefore compute the derivatives as follows:

$$\begin{array}{rl}\frac{\partial E_j}{\partial r}& =\frac{\partial u_j}{\partial r}{Z}_j(z)\\ \frac{{\partial }^2{E}_j}{\partial r^2}& =\frac{{\partial }^2{u}_j}{\partial r^2}{Z}_j(z)\\ \frac{{\partial }^2{E}_j}{\partial {\varphi }^2}& =\frac{{\partial }^2{u}_j}{\partial {\varphi }^2}{Z}_j(z)\\ \frac{{\partial }^2{E}_j}{\partial z^2}& =u_j\frac{d^2{Z}_j}{d{z}^2}\end{array}$$

Thus if we substitute our generic form of $E$ into the Helmholtz equation in cylindrical coordinates, we obtain:

$$\begin{array}{rl}{\nabla }^2{E}_j+k^2{E}_j& =\left[\frac{1}{r}\left(\frac{\partial }{\partial r}+r\frac{{\partial }^2}{\partial r^2}\right)+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\varphi }^2}+\frac{{\partial }^2}{\partial z^2}\right]E_j+k^2{E}_j\\ & =\frac{1}{r}\frac{\partial u_j}{\partial r}{Z}_j+\frac{{\partial }^2{u}_j}{\partial r^2}{Z}_j+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\varphi }^2}Z+u_j\frac{d^{2}Z}{d{z}^2}+k^2{u}_j{Z}_j\\ & =\left(\frac{1}{r}\frac{\partial u_j}{\partial r}+\frac{{\partial }^2{u}_j}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2{u}_j}{\partial {\varphi }^2}\right)Z_j+u_j{Z}_j^{\prime \prime }+k^2{u}_j{Z}_j\\ & =0\end{array}$$

Now dividing by $u_j{Z}_j$ from all sides one obtains:

$$\begin{array}{c}\frac{1}{u_j{Z}_j}\left[\left(\frac{1}{r}\frac{\partial u_j}{\partial r}+\frac{{\partial }^2{u}_j}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2{u}_j}{\partial {\varphi }^2}\right)Z_j+u_j{Z}_j^{\prime \prime }+k^2{u}_j{Z}_j\right]=0\\ \frac{1}{u_j}\left(\frac{1}{r}\frac{\partial u_j}{\partial r}+\frac{{\partial }^2{u}_j}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2{u}_j}{\partial {\varphi }^2}\right)+\frac{Z_j^{\prime \prime }}{Z_j}=-k^2\\ \frac{1}{u_j}\left(\frac{1}{r}\frac{\partial u_j}{\partial r}+\frac{{\partial }^2{u}_j}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2{u}_j}{\partial {\varphi }^2}\right)=-\frac{Z_j^{\prime \prime }}{Z_j}-k^2=-m^2\end{array}$$

Therefore, one now obtains two differential equations, as follows:

$$\begin{array}{c}{Z}_j^{\prime \prime }=-(k^2-m^2)Z_j\\ \left(\frac{1}{r}\frac{\partial u_j}{\partial r}+\frac{{\partial }^2{u}_j}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2{u}_j}{\partial {\varphi }^2}\right)u_j=-m^2{u}_j\end{array}$$

The first differential equation is straightforward to solve; the solution is given by:

$$Z_j(z)=C_1{e}^{i(\sqrt{k^2-m^2})z}+C_2{e}^{-i(\sqrt{k^2-m^2})z}$$

It is often more convenient to use the real-valued form:

$$Z_j(z)=C_1\cos (\sqrt{k^2-m^2}z)+C_2\sin (\sqrt{k^2-m^2}z)$$

For the second differential equation, we may solve by performing another separation of variables. First, we can multiply through by $r^2$ and rearrange, giving us:

$$r^2\frac{{\partial }^{2}u}{\partial r^2}+r\frac{\partial u}{\partial r}+\frac{{\partial }^{2}u}{\partial {\varphi }^2}=-r^2{m}^2{u}_j$$

We may solve this differential equation by separation of variables. Let $u=R(r)\Phi (\varphi )$. Then the derivatives of $u$ are thus given by:

$$\frac{\partial u}{\partial r}=\Phi \frac{dR}{dr},\frac{{\partial }^{2}u}{\partial r^2}=\Phi \frac{d^{2}R}{d{r}^2},\frac{{\partial }^{2}u}{\partial {\varphi }^2}=R\frac{d^2\Phi }{d{\varphi }^2}$$

Upon substitution, we have:

$$r^2\Phi \frac{d^{2}R}{d{r}^2}+r\Phi \frac{dR}{dr}+R\frac{d^2\Phi }{d{\varphi }^2}=-r^2{m}^{2}R\Phi$$

Now dividing by $R(r)\Phi (\varphi )$ from both sides, we have:

$$\begin{array}{c}\frac{r^2}{R}\frac{d^{2}R}{d{r}^2}+\frac{r}{R}\frac{dR}{dr}+\frac{1}{\Phi }\frac{d^2\Phi }{d{\varphi }^2}=-r^2{m}^2\\ \Rightarrow \frac{r^2}{R}\frac{d^{2}R}{d{r}^2}+\frac{r}{R}\frac{dR}{dr}+m^2{r}^2=-\frac{1}{\Phi }\frac{d^2\Phi }{d{\varphi }^2}=n^2\end{array}$$

Where $n^2$ is a separation constant. This gives us two differential equations to solve:

$$\begin{array}{rl}\frac{d^2\Phi }{d{\varphi }^2}& =-n^2\Phi \\ r^2\frac{d^{2}R}{d{r}^2}& +r\frac{dR}{dr}+(m^2{r}^2-n^2)R=0\end{array}$$

The first differential equation is trivial to solve; its solution is, by inspection:

$$\Phi (\varphi )=A\cos (n\varphi )+B\sin (n\varphi )$$

Meanwhile, we identify the second differential equation as Bessel’s differential equation, whose solution are given in terms of the Bessel functions:

$$R(r)=C{J}_n(m|r|)+D{Y}_n(m|r|)$$

Where $J_n$ is the $n$-th Bessel function of the first kind, and $Y_n$ is the $n$-th Bessel function of the second kind. Putting it all together, we thus have:

$$\begin{array}{rl}{E}_j(r,\varphi ,z)& =R(r)\Phi (\varphi )Z(z)\\ & =[A\cos (n\varphi )+B\sin (n\varphi )][C{J}_n(m|r|)+D{Y}_n(m|r|)]\\ & \times [F\cos (\sqrt{k^2-m^2}z)+G\sin (\sqrt{k^2-m^2}z)]\end{array}$$

Where $A,B,C,D,F,G$ are all undetermined coefficients. The most general solution is obtained by summing over all possible $m$ and $n$, giving us:

$$\begin{array}{rl}{E}_j(r,\varphi ,z)=\sum _{n=0}^{\infty }\sum _{m=0}^{\infty }[& A_n\cos (n\varphi )+B_n\sin (n\varphi )]\\ & \times [C_{nm}{J}_n(m|r|)+D_{nm}{Y}_n(m|r|)]\\ & \times [F_m\cos (\sqrt{k^2-m^2}z)+G_m\sin (\sqrt{k^2-m^2}z)]\end{array}$$

We may now apply specific boundary conditions to obtain the particular solutions for $E_r$ and $E_{\varphi }$. To begin, our requirement of azimuthal symmetry requires that we impose the periodic boundary condition $E_j(r,\varphi ,z)=E_j(r,\varphi +2\pi ,z)$. This is satisfied by default as long as $n=0,1,2,3,\dots$ is some positive integer value and if either one of $A_m,B_m$ is zero; if we specifically specify that $E_j(r,0,z)=0$ then . Meanwhile, we also impose the requirement that the solution does not blow up for any value of $r$. Because of how $Y_n$ (the Bessel functions of the second kind) are constructed, they will lead to asymptotic behavior at $r=0$, meaning that we must have $D_{mn}=0$ to preserve a physical solution. This simplifies our particular solution to:

$$\begin{array}{rl}{E}_j(r,\varphi ,z)& =\sum _{n=0}^{\infty }\sum _{m=0}^{\infty }[A_n\cos (n\varphi )+B_n\sin (n\varphi )]C_{nm}{J}_n(m|r|)\\ & \times [F_m\cos (\sqrt{k^2-m^2}z)+G_m\sin (\sqrt{k^2-m^2}z)]\end{array}$$

Note: In this particular solution, a time-dependent factor of $e^{-i\omega t}$ is implied; it is not very physically important in RF analysis however, since the oscillation frequencies are in the GHz range (that is, billions of times per second) and thus we usually care only about the spatial variation of the EM field, not its temporal variation.

Note that this particular solution automatically satisfies the open boundary conditions at $r=\pm w/2$ since the solution smoothly transitions over the sides, meaning that physically, EM waves are allowed to freely pass through. It is also physically consistent since ${\lim }_{r\to \infty }{J}_n(m|r|)=0$, i.e. the solution is well-defined for all values of $r$ (and it is not hard to show that it is also well-defined for all values of $z\in [0,L]$ and $\varphi \in [0,2\pi ]$ within our domain).

We must now impose specific boundary conditions that are different for $E_r$ and $E_{\varphi }$, based on the physics of the problem. First, we have the transmission boundary condition $E_r(r,\varphi ,L)=T{E}_r(r,\varphi ,0)$. For this to be the case, we must have:

$$G_m\sin (\sqrt{k^2-m^2}L)+F_m\cos (\sqrt{k^2-m^2}L)=T{F}_m$$

Dividing by $F_m$ on all sides and defining $\rho \equiv G_m/F_m$ we obtain:

$$\rho \sin (\sqrt{k^2-m^2}L)+\cos (\sqrt{k^2-m^2}L)=T$$

Now defining $\overline{\phi }=\sqrt{k^2-m^2}L$, we have:

$$\rho =\frac{T-\cos \overline{\phi }}{\sin \overline{\phi }}$$

This is, in fact, nearly the exact same equation as our previous 2D result, with the only difference being that we use $\overline{\phi }$ rather than $\phi$. Therefore, the roots ${\overline{\phi }}_n^{\ast }$ are also almost exactly the same, and are given by:

$${\overline{\phi }}_n^{\ast }=\{\begin{array}{ll}n\pi -{\cos }^{-1}(T),& n=2,4,6,8,\dots \\ (n-1)\pi +{\cos }^{-1}(T),& n=1,3,5,7,\dots \end{array}$$

Solving for the values of $k$ that satisfy ${\overline{\phi }}_n^{\ast }=\sqrt{k^2-m^2}L$, we find that:

$$k=\sqrt{{\left(\frac{{\overline{\phi }}_n^{\ast }}{L}\right)}^2+m^2},n,m=1,2,3,\dots$$

This is rather surprising in that in this case, we do not have a restriction on $m$ from the square root (as we saw before); rather, the values of $k$ are restricted based on the fixed values of $m$.

More on electromagnetic boundary conditions

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 ${ϵ}_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 ${ϵ}_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 ${ϵ}_r^{\prime }$. 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 },{ϵ}_r{E}_{\mathrm{interior}}^{\perp }={ϵ}_r^{\prime }{E}_{\mathrm{mirror}}^{\perp }$$

Fourth, we know that the beam has azimuthal symmetry, and therefore we naturally have the periodic boundary condition $E_{\varphi }(r,\varphi ,z)=E_{\varphi }(r,\varphi +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_{\varphi }$ and $E_r$ (the radial and azimuthal components of the field); moreover, $E_{\mathrm{interior}}^{\parallel }=E_{\mathrm{mirror}}^{\parallel }=0$.

Now, upon performing separation of variables in cylindrical coordinates, we previously obtained a particular solution that could be written as a sum of Bessel functions $J_m(x)$ (more precisely, the Bessel functions of the first kind)¹⁰. That is, we have found a solution in the form:

$$E(r,\varphi ,z)=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{A}_{mn}{J}_n(r\sqrt{k^2+m^2})\cos (n\varphi )e^{-i\sqrt{k^2-m^2}z}$$

Which can be decomposed into modes $E_n(r,\varphi ,z)$ as follows:

$$\begin{array}{rl}E(r,\varphi ,z)& =\sum _{n=0}^{\infty }{E}_n(r,\varphi ,z)\\ E_n(r,\varphi ,z)& =\sum _{m=0}^{\infty }{A}_{mn}{J}_n(r\sqrt{k^2+m^2})\cos (n\varphi )e^{-i\sqrt{k^2-m^2}z}\end{array}$$

To replicate the Gaussian beam, we replace the amplitude coefficients $A_{mn}$ by a continuous function $A(r,z)$ (which does not depend on $\Phi$ due to azimuthal symmetry). We may now proceed using a combination of a heuristic line of reasoning and mathematical analysis. Here, we will describe both.

In the heuristic approach, note that for the solution to become zero at $z=0$ and $z=L$, as well as at $r=a$ (for the rightward mirror), it would make sense for $A(r,z)$ to be the product of a Gaussian and a sinusoid (since $\sin (0)=0$ and for a resonant cavity it is always the case that $\sin (kL)\to 0$ in the case of zero transmission). More precisely, we expect the Gaussian to have a radius (or more precisely, RMS width) of $a$. Thus we assume an ansatz of the form:

$$A(r,z)=A(r)=E_0{e}^{-r^2/a^2}$$

Where we may write $A(r,z)$ as simply $A(r)$ since our ansatz does not depend explicitly on $z$. Now, to obtain a Gaussian beam solution, we must replace the infinite sum in the solution with a continuous integral. The reason is because just as we have replaced the constant coefficients $A_{mn}$ with a continuous function $A(r)$, so we must also transform a discrete sum into a continuous integral for the solution to be mathematically consistent. By substituting in our ansatz for $A(r,z)$ and swapping the sum for an integral, we obtain:

$$\begin{array}{rl}{E}_n(r,\varphi ,z)& =\int dkA(r)J_n(r\sqrt{k^2+m^2})\cos (n\varphi )e^{-i\sqrt{k^2-m^2}z}\\ & =\cos (n\varphi )e^{-i\sqrt{k^2-m^2}z}\int dkA(r)J_n(r\sqrt{k^2+m^2})\end{array}$$

The integral we see is in the form of an inverse Hankel transform. If the transform is performed, we find that the fundamental mode ($n=0$) for our solution is just a (simplified) Gaussian beam!

Note: Our solution is not precisely a Gaussian beam because we assumed in our calculation that the width of the beam is always constant. This cannot physically occur for $z>L$ since the diffraction limit tells us that all electromagnetic waves passing through an aperture must eventually diverge.

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.

Note: Add the effects of impedance to describe a lossy conductor. See Ch. 1.7 in Pozar, Microwave Engineering for more information.

Numerical solution of the realistic model

Note: Finite element simulations should be added, with streamline plots and/or vector field plots of both the electric field and magnetic field, as well as the radiation pattern plot (magnitude of the Poynting vector)

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,\varphi ,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 $({\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{array}{r}\frac{{\partial }^2{E}_{\perp }}{\partial r^2}+\frac{{\partial }^2{E}_{\perp }}{\partial z^2}+k^2{E}_{\perp }=0\\ \frac{{\partial }^2{E}_{\parallel }}{\partial r^2}+\frac{{\partial }^2{E}_{\parallel }}{\partial z^2}+k^2{E}_{\parallel }=0\end{array}$$

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¹ 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¹⁰, 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\frac{\partial u}{\partial n}+c_3=0$$

Substituting our solution $u(\mathbf{x})$ in, we have:

$$\begin{array}{c}{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}}+r{e}^{-ik\cdot \mathbf{x}})+ik{c}_2(e^{ik\cdot \mathbf{x}}-r{e}^{-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=-\frac{c_3+(c_1-ik{c}_2)e^{ik\cdot \mathbf{x}}}{(c_1+ik{c}_2)e^{-ik\cdot \mathbf{x}}}\end{array}$$

Thus, we have:

$$\begin{array}{rl}R& =|r{|}^2\\ & ={\left|\frac{c_3+(c_1-ik{c}_2)e^{ik\cdot \mathbf{x}}}{(c_1+ik{c}_2)e^{-ik\cdot \mathbf{x}}}\right|}_{(\mathbf{x}=\partial \Omega =L)}^2\\ & ={\left|1+\frac{c_3}{(c_1+ik{c}_2)e^{-ikL}}\right|}^2\end{array}$$

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+\frac{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{array}{rl}R& ={\left|1+\frac{\alpha c_1}{c_1}{e}^{ikL}\right|}^2\\ & ={\left|1+\alpha e^{ikL}\right|}^2\\ & =1+{\alpha }^2+2\alpha \cos (kL)\end{array}$$

The solution can be found using the quadratic formula:

$$\begin{array}{c}{\alpha }^2+2\alpha \cos (kL)+(1-R)=0\\ \Rightarrow \alpha =-\cos (kL)-\sqrt{{\cos }^2(kL)-(1-R)}\end{array}$$

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{array}{c}{c}_{1}u(\mathbf{x})+c_2\frac{\partial u}{\partial n}+c_3=0\\ \Rightarrow c_{1}u(\mathbf{x})+\alpha c_1=0\\ \Rightarrow u(L)=-\alpha \end{array}$$

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 =633\mathrm{nm}$ light, substituting in $k=2\pi /\lambda$ and $L=30\mathrm{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)+\underset{\text{only real for }{\cos }^2(kL)=1}{\underbrace{\sqrt{{\cos }^2(kL)-1}}}=\cos (kL)$$

Which describes outgoing plane waves $u(x)=\cos (kx)=\mathrm{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 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(\frac{\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{array}{c}{E}_{\perp }(r,0)=0,E_{\perp }(r,L)=\eta E_0\\ \frac{\partial E_{\perp }(z,R)}{\partial r}-ik{E}_{\perp }(z,R)=0,\\ \frac{\partial E_{\perp }(z,-R)}{\partial r}-ik{E}_{\perp }(z,-R)=0,\\ r\in [0,R],z\in [0,L]\end{array}$$

• The aperture has a Sommerfeld-like radiative boundary condition
• The two “mirrors” of the RF cavity (really, they’re just metal plates and don’t have to be polished) have a reflective (homogeneous Dirichlet) boundary condition
• The rest of the RF cavity (including the “tube” connecting the two mirrors) also has a Sommerfeld-like radiative boundary condition; this represents a material transparent to microwaves

In this case, it is easier to use cylindrical coordinates $(r,\Phi ,z)$, where $z$ is the optical axis (direction of propagation) and $(r,\Phi )$ are the radial and polar angles respectively. We again assume our cavity to have length $L$. Our boundary conditions for the Helmholtz equation thus correspond mathematically to:

$$\begin{array}{rlr}& E(r,\Phi ,0)=0,& 0\le r\le R,\\ & E(r,\Phi ,L)=0,& a<r\le R,\\ & \underset{r\to {\infty }^{+}}{\lim }E(r,\Phi ,z)=0,& r>R,\\ & \underset{z\to {\infty }^{+}}{\lim }E(r,\Phi ,z)=0,& 0\le r\le a,\\ & E(r,\Phi ,z)=E(r,\Phi +2\pi ,z),& 0\le \Phi \le 2\pi \end{array}$$

Where our domain is given by $z\in [0,z_{\text{max}}],r\in [0,r_{\text{max}}],\Phi \in [0,2\pi ]$, and we want to pick a suitably large simulation domain (where $z_{\text{max}}\ge 5L$ and $r_{\text{max}}\ge 5R$, roughly) to be able to see both the near-field and the far-field behavior of the beam. Our first two boundary conditions encode the fact that the electromagnetic waves in the cavity must reflect off the two mirrors at each end of the cavity. Meanwhile, our last two boundary conditions encode the fact that the electromagnetic waves must pass freely through the tube connecting the two mirrors, as well as the aperture (output coupler) at the end of the cavity. Our last boundary condition enforces azimuthal symmetry since we expect our solution to be symmetric with respect to $\Phi$. While it may appear that energy would be needlessly dissipated over the sides of the maser, this actually does not happen, since the electromagnetic field quickly becomes zero for $r>R$.

Note: This means that the resultant electric field will be complex-valued; we take only the real part for a physically-meaningful solution.

Maser gain analysis

Lastly, we don’t just want to see what the beam’s wave profile; we also want to know how well the maser amplifies light¹¹. The corresponding quantity of interest to us in the numerical simulation, which describes this amplification, is called the gain of the maser. The gain compares the output power of the maser to that of an isotropic (point-source) radiator. It is usually measured in decibels (dB), which are a logarithmic unit, and it can be written as:

$$G_{\mathrm{dB}}(\Phi ,\varphi )={\log }_{10}\left(\frac{P_{\text{maser}}(\Phi ,\varphi )}{P_{\text{iso}}(\Phi ,\varphi )}\right)={\log }_{10}\left(\frac{P_{\text{maser}}(\Phi ,\varphi )4\pi |\mathbf{r}{|}^2}{P_{\text{input}}}\right)$$

Where $P_{\text{iso}}=\frac{P_{\text{input}}}{4\pi |\mathbf{r}{|}^2}$ (in cylindrical coordinates we have $|\mathbf{r}|=\sqrt{r^2+z^2}$), $P_{\text{input}}$ is the input power used to drive the maser, the power of the maser is given by $P=|\mathbf{S}|$, and $\mathbf{S}$ is just the Poynting vector. The input power $P_{\text{input}}=(I\Delta V{)}_{\text{input}}$, where $I$ and $\Delta V$ are the driving current and voltage used to power the maser respectively. If we substitute our analytical solution, we find that:

$$P_{\text{maser}}\sim I_0{e}^{-2r^2/a^2},I_0\sim E_0^2$$

Where here, $I_0$ is the maximum intensity of the beam. Thus the gain is approximately:

$$G_{\text{dB}}\sim {\log }_{10}\left(\frac{4\pi I_0|\mathbf{r}{|}^2{e}^{-2r^2/a^2}}{P_{\text{input}}}\right)$$

If we assume that the initial intensity $I_0$ is equal to the input power $P_{\text{input}}$ (that is, we assume that the maser is 100% efficient mechanically-wise), this gives us:

$$G_{\text{dB}}\sim {\log }_{10}(4\pi |\mathbf{r}{|}^2{e}^{-2r^2/a^2})$$

Note that our result is independent of angle and also independent of $z$. However, this is not entirely accurate, as again we are using a simplified solution as opposed to the physically-accurate solution (the Gaussian beam, where $w(z)\ne a$). For a Gaussian beam, the gain becomes:

$$G_{\text{dB}}\sim {\log }_{10}\left(4\pi \frac{(|\mathbf{r}|w_0{)}^2}{w(z{)}^2}{e}^{-2r^2/w(z{)}^2}\right),w_0=a,w(z)\approx \frac{w_{0}z}{L}$$

Where we used the fact that $w(z)\sim z$ for $z\gg L$. This result can be simplified and written explicitly, up to some constants, to:

$$G_{\text{dB}}\sim {\log }_{10}\left(\frac{r^2+z^2}{z^2}{e}^{-2r^2/z^2}\right)\sim {\log }_{10}\left[\left(1+\frac{r^2}{z^2}\right)e^{-2r^2/z^2}\right]$$

It is easy to see that in the limit as $z\to 0$, we have $G_{\text{dB}}\to \infty$. Thus, at its source, a (perfect) Gaussian beam would have effectively infinite gain, which corresponds to infinite amplification of light (once again, here “amplification” means “concentration of energy” not “increasing energy”). Of course, as $z$ increases, the gain of the beam decreases; at some finite distance, $G_{\text{dB}}=0$ (meaning that the Gaussian beam has diverged to the point that it is effectively isotropic and no longer performs useful amplification of light). A minimally-divergent Gaussian beam (whose profile does not change much as $z$ increases) would thus have very high gain; a highly-divergent Gaussian beam would have very low gain. Therefore, decreasing the divergence of a Gaussian beam is absolutely fundamental to improving its amplification and thus its efficiency.

Footnotes

1. This is due to the transverse and longitudinal reflections interfering with each other, leading to the buildup of undesired non-Gaussian modes in the RF cavity. In simplified terms, any off-optical-axis reflections will interfere with the desirable on-axis reflections, causing power to be dispersed and the beam to lose directionality. ↩ ↩²

2. From Wolfram Mathworld’s article article on solving the Helmholtz equation in cartesian coordinates. Owing to our coordinate conventions, we make the substitutions $x\to r,y\to z$. ↩

3. This follows the derivation in Pozar, Microwave Engineering (4th ed.) in Chapter 3.2. ↩

4. For more information, see the Wikipedia article on RF waveguides ↩

5. This was inspired by the boundary conditions described in Freund & Antonssen’s Principles of Free Electron Lasers (4th ed.), ch. 9.1-9.4 as well as the discussion on the transmission coefficient $T_x$ in hole-coupled free-electron lasers from Varro et. al., Free Electron Lasers ,Ch. 3, pg. 68 ↩

6. Directly based off Varro et. al., Free Electron Lasers ,Ch. 3, pg. 68 ↩

7. This is only a basic overview of the finite element method - more details and a much more comprehensive guide to the finite element method can be found on the COMSOL software blog. ↩

8. This comes from the form shown in the documentation for the pyoomph library. Note that we swap the order of the inner products, but this does not affect the results since the inner product is commutative. ↩

9. Taken from an oomph library example for solving the Helmholtz equation with the finite element method. ↩

10. This is a simplified form of the solution given on Wolfram Mathworld’s article on solving the Helmholtz equation in cylindrical coordinates. In particular, we choose $B_{mn}=F_n=0$ so that the solution does not blow up (become infinite) at $z\to \infty$ and at $z\to 0$. ↩ ↩²

11. Technically-speaking, a maser (or laser) doesn’t amplify light; at least, not in the sense of producing more (radiant) energy than the energy used to power it. Rather, “amplification” is a rough term that describes the fact that a maser concentrates light so that it is amplified in one particular direction, while reduced in all other directions. ↩