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. ## A basic model We consider a RF cavity (maser cavity) of aperture width $a$, total width $w$, and total length $L$. We'll assume the cavity is square and of unit length $L = 1$ for simplicity. The (scalar) Helmholtz equation for the electric field $E(\mathbf{x})$ is given by: $ (\nabla^2 + k^2)E = 0 $ Here, we'll mostly focus on the 2D case, where the Helmholtz equation takes the form: $ \dfrac{\partial^2 E}{\partial x^2} + \dfrac{\partial^2 E}{\partial y^2} + k^2 E = 0 $ Using analytical and numerical means, we'll tackle this problem in detail. ### Analytical solution Surprisingly, the problem actually has an analytical solution if we switch to cylindrical coordinates $(\rho, \phi, x)$. This method uses the [Fourier optics approach](https://en.wikipedia.org/wiki/Fourier_optics) which is rather complicated, so we'll not discuss it in-depth. However, the end-result is that the radiation pattern of the power density of the electromagnetic field is given (up to some constant factors) by[^3]: $ I(\rho, x) \sim \left (\frac{2J_1(a\pi \rho /\lambda x)}{a\pi \rho/\lambda x} \right)^{2} $ Where we choose $x$ as the optical axis on which the microwaves propagate, $a$ once again is the aperture width (width as in _diameter_, not radius), $\lambda$ is the wavelength of the microwaves, and $J_1(\rho)$ is a 1st-order [Bessel function of the first kind](https://en.wikipedia.org/wiki/Bessel_function). Here is an [interactive Desmos plot](https://www.desmos.com/3d/cqtiu6nvym) of this radiation pattern, and a plot with $a = 0.1, \lambda=1$ is shown below: ![[cavity-with-aperture.png]] > **Note:** This was generated with the `rf-cavity-with-aperture/cavity_with_aperture.m` Octave/MATLAB script in the `visualizations/` folder of this repository. The important result is that the wave diverges to infinity and becomes unfocused very quickly after leaving the aperture, which can be seen in the plot, but can also be shown analytically (if we take the limit as $x \to \infty$ we have $I \to 0$ for all $\rho$). ### Numerical solution To solve the Helmholtz equation numerically, we use the **finite element method** with the open-source FreeFEM software. This is essential when the cavity is not a perfect cylindrically-symmetric cavity (as would be the case for any realistic RF cavity design). 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_[^1]. 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^2E + 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 u dv = uv\bigg|_a^b - \int_a^b vdu $ The multivariable version of integration by parts takes the form: $ \int_\Omega u dv = \int_{\partial \Omega} uv - \int_\Omega vdu $ Here, we apply integration by parts to our term $\displaystyle \int_\Omega \phi \nabla^2 E$. Let $u = \phi$ and $dv = \nabla^2 E$. Then, we have $du = \nabla \phi$ and $v = \displaystyle \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$ so that we have the simplified result: $ \int_{\partial \Omega} \phi \nabla E - \int_\Omega \nabla \phi \cdot \nabla E - k^2 \phi 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). 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 typical methods of 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. ### Boundary conditions An RF cavity has reflective walls made of metal (which we can assume to be a perfect conductor), so it can be modelled by reflective (Dirichlet) boundary conditions: $ \begin{gather*} E(x, 0) = E(x, L) = 0 \\ E(0, y) = E(L, y) = 0 \end{gather*} $ For the aperture, unfortunately, things are not so simple. To accurately model outgoing waves that radiate to infinity, we must use the **Sommerfeld radiation condition**, which, in $d$ dimensions, takes the form: $ \lim_{r \to \infty} r^{\frac{d-1}{2}} \left(\dfrac{\partial}{\partial r} - ik\right)u = 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{align*} \lim_{r \to \infty} \sqrt{r} \left(\dfrac{\partial}{\partial r} - ik\right)u = 0, &\quad & \text{(2D)} \\ \lim_{r \to \infty} r \left(\dfrac{\partial}{\partial r} - ik\right)u = 0, &\quad &\text{(3D)} \end{align*} $ One can approximate the 2D Sommerfield radiation condition up to first-order with[^2]: $ \dfrac{\partial u}{\partial r} - iku = 0 $ The requirement for this to be valid, 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. ![[helmholtz-RF-cavity.excalidraw.svg|400]] Note that this idea can easily be extended into 3D: the only difference is that the circular boundary will need to be changed into a spherical boundary. ### Simulation results The FreeFEM code for the simulation can be found in `freefem/dev/lasers/solve-Efield.edp`. The result is shown below: ![[rf-cavity.png|400]] 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 more advanced model In our basic model of a maser RF cavity, we assumed that the cavity had perfectly-reflecting metal walls. This, of course, is an oversimplification. Real masers do not use perfectly-reflecting RF cavities; rather, they use _open_ cavities where two metal plates are placed some distance apart from each other, but the tube joining them is made of a microwave-transparent material. This allows microwaves to be reflected _between_ the metal plates, but *not* reflect off the walls, or otherwise the directionality of the maser would be ruined[^5]. This means that we need to change our boundary conditions to the following: - 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 > **Note:** If the domain is very large, it is possible to approximate the radiative boundary condition with simpler periodic boundary conditions. See [the visualization here](https://visualpde.com/sim/?mini=M93Igvaz) for a numerical simulation using periodic boundary conditions (this simulation assumes no maser output coupler). In this case, it is easier to use cylindrical coordinates $(r, \theta, z)$, where $z$ is the optical axis (direction of propagation) and $(r, \theta)$ 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{align*} &E(r, \theta, 0) = 0, & 0 \leq r \leq R, \\ &E(r, \theta, L) = 0, & a < r \leq R, \\ &\lim_{r \to \infty^+} E(r, \theta, z) = 0, & r > R, \\ &\lim_{z \to \infty^+} E(r, \theta, z) = 0, &0 \leq r \leq a, \\ &E(r, \theta, z) = E(r, \theta + 2\pi, z), & 0 \leq \theta \leq 2\pi \end{align*} $ Where our domain is given by $z \in [0, z_\text{max}], r \in [0, r_\text{max}], \theta \in [0, 2\pi]$, and we want to pick a suitably large simulation domain (where $z_\text{max} \geq 5L$ and $r_\text{max} \geq 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 $\theta$. 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$. ### The analytical solution To solve the boundary-value problem for the Helmholtz equation, we first note that its general solution in cylindrical coordinates can be written as a sum of Bessel functions, that is[^6]: $ E(r, \theta, z) = \sum_m \sum_n A_{mn} J_m(r \sqrt{k^2 + n^2}) \cos (m \theta) e^{-nz} $ If we assume that the dominant mode is given by $m = n = 0$, our solution simplifies to: $ E(r, \theta, z) = A(r, z) J_0(kr) $ Where $A_{mn}$ is now replaced by the function $A(r, z)$ (which does not depend on $\theta$ due to azimuthal symmetry). Now, 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) = 0$). More precisely, we expect the Gaussian to have a radius (or more precisely, [RMS width](https://en.wikipedia.org/wiki/Root_mean_square)) of $a$. Of course, a Gaussian never truly decays to zero, but it decays to zero "fast enough" that it effectively satisfies the boundary conditions. Thus we assume a solution of the form: $ A(r, z) = E_0 e^{-(r^2/a^2 + z^2)} \sin(kz) \quad \Rightarrow \quad E(r, \theta, z) \sim E_0 e^{-(r^2/a^2 + z^2)}J_0(kr)e^{-ikz} $ Where we use complex exponentials instead of real sinusoids when substituting into the expression for the field (so $\sin kz \to e^{-ikz}$). We find that 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**. ### Maser gain analysis For this more advanced analysis, we don't just want to see what the beam's wave profile; we also want to know how well the maser _amplifies_ light[^4]. 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}(\theta, \phi) = \log_{10}\left(\dfrac{P_\text{maser}(\theta, \phi)}{P_\text{iso}(\theta, \phi)}\right) = \log_{10} \left( \dfrac{P_\text{maser}(\theta, \phi)4\pi |\mathbf{r}|^2}{P_\text{input}}\right) $ Where $P_\text{iso} = \dfrac{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}, \quad 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(\dfrac{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) \neq a$). For a Gaussian beam, the gain becomes: $ G_\text{dB} \sim \log_{10}\left(4\pi\dfrac{(|\mathbf{r}|w_0)^2}{w(z)^2} e^{-2r^2/w(z)^2}\right), \quad w_0 = a, \quad w(z) \approx \dfrac{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(\dfrac{r^2 + z^2}{z^2} e^{-2r^2/z^2}\right) \sim \log_{10}\left[\left(1 + \dfrac{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. [^1]: 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](https://www.comsol.com/multiphysics/finite-element-method). [^2]: Taken from an [oomph library example](https://oomph-lib.github.io/oomph-lib/doc/helmholtz/scattering/html/index.html) for solving the Helmholtz equation with the finite element method. [^3]: This comes from [this optics website](https://opticsthewebsite.com/Circular), which additionally goes into more depth about Fourier optics. [^4]: 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. [^5]: 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. [^6]: This is a simplified form of the solution given on Wolfram Mathworld's article on [solving the Helmholtz equation in cylindrical coordinates](https://mathworld.wolfram.com/HelmholtzDifferentialEquationCircularCylindricalCoordinates.html). 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$.