A classical field theory formulation for the numerical solution of time harmonic electromagnetic fields

TL;DR

This paper introduces a classical field theory-based variational formulation for solving time-harmonic electromagnetic fields, simplifying finite element methods by replacing edge elements with nodal elements, and leveraging gauge fixing and differential geometry.

Contribution

It presents a novel variational approach inspired by classical field theory and QED, enabling the use of simpler nodal elements in finite element analysis of Maxwell's equations.

Findings

Replaces edge elements with nodal elements for simplicity.

Provides a robust, symmetry-preserving numerical framework.

Demonstrates effectiveness in time-harmonic electromagnetic problems.

Abstract

Finite element representations of Maxwell's equations pose unusual challenges inherent to the variational representation of the `curl-curl' equation for the fields. We present a variational formulation based on classical field theory. Borrowing from QED, we modify the Lagrangian by adding an implicit gauge-fixing term. Our formulation, in the language of differential geometry, shows that conventional edge elements should be replaced by the simpler nodal elements for time-harmonic problems. We demonstrate how this formulation, adhering to the deeper underlying symmetries of the four-dimensional covariant field description, provides a highly general, robust numerical framework.

Figures

Equations

$$F={\int }_{\Omega }\frac{1}{{\mu }_r}|\nabla \times \vec{E}{|}^2-k_0^2{ϵ}_r|\overrightarrow{E_r}{|}^2+i{k}_0{Z}_0({\vec{E}}^{\ast }\cdot \vec{J}-\vec{E}\cdot {\vec{J}}^{\ast })\mathrm{d}V$$

$$\mathcal{L}=-\frac{1}{2\mu }\mathrm{d}\mathbf{A}\wedge \star \mathrm{d}\mathbf{A}+\mathbf{A}\wedge \mathbf{J}-\frac{1}{2\xi \mu }\mathrm{d}\star \mathbf{A}\wedge \star \mathrm{d}\star \mathbf{A}$$

$$\mathrm{d}\mathbf{F}=0$$

$$\mathrm{d}\star \mathbf{F}=\mathbf{J}$$

$$\mathcal{L}=-\frac{1}{2\mu }\mathrm{d}\mathbf{A}\wedge \star \mathrm{d}\mathbf{A}+\mathbf{A}\wedge \mathbf{J}$$

$${\mathcal{L}}_{\mathrm{GF}}=-\frac{1}{2\mu \xi }\mathrm{d}\star \mathbf{A}\wedge \star \mathrm{d}\star \mathbf{A}$$

$$\begin{gathered} S(\varphi ,\vec{A})=\frac{1}{2}\int ϵ|\nabla \varphi +i\omega \vec{A}{|}^2-\frac{1}{\mu }|\nabla \times \vec{A}{|}^2 \\ -\frac{1}{\mu \xi }|\nabla \cdot \vec{A}-\frac{i\omega }{c^2}\varphi {|}^2-\rho {\varphi }^{\ast }-{\rho }^{\ast }\varphi +\vec{A}\cdot {\vec{J}}^{\ast }+{\vec{A}}^{\ast }\cdot \vec{J}\mathrm{d}V \end{gathered}$$

$${\nabla }^2\varphi +k_0^2\varphi =-\frac{\rho }{ϵ}$$

$${\nabla }^2\vec{A}+k_0^2\vec{A}=-\mu j$$

$$S_{\mathrm{Z}}=\gamma {\int }_{\partial \Omega }{\left|\hat{n}\times \left(-\nabla \varphi -i\omega \vec{A}\right)\right|}^{2}dS$$

$$\alpha =\frac{1}{\xi }$$

$$\alpha n^2{k}_0^2\varphi +{\nabla }^2\varphi -i{k}_0(\alpha -1)\nabla \cdot \vec{A}$$

$$n^2{k}_0^2\vec{A}+{\nabla }^2\vec{A}+(\alpha -1)\nabla (\nabla \cdot \vec{A}+i\frac{\omega }{c^2}\varphi )$$

$$\left(\mathbf{M}{\widetilde{k_0}}^2+\mathbf{C}\widetilde{k_0}+\mathbf{K}\right).\mathbf{a}=0$$

$$\left(\mathbf{M}{\widetilde{k_0}}^2+\mathbf{C}\widetilde{k_0}+\mathbf{K}\right).\mathbf{a}=\mathbf{j}$$

$$\vec{\mathcal{E}}(\vec{r},t)=\mathrm{Re}[\vec{E}(\vec{r})e^{i\omega t+im\theta }]$$

$$\mathcal{L}=-\frac{1}{4\mu }{F}^{\nu \beta }{F}_{\nu \beta }-A_{\nu }{J}^{\nu }-\frac{\alpha }{2}{\left({\partial }_{\nu }{A}^{\nu }\right)}^2$$

$$\begin{gathered} \mathcal{L}=\frac{1}{2}[ϵ(\nabla \varphi +i\omega \vec{A})\cdot (\nabla {\varphi }^{\ast }-i\omega {\vec{A}}^{\ast })-\frac{1}{\mu }(\nabla \times \vec{A})\cdot (\nabla \times {\vec{A}}^{\ast }) \\ -\frac{\alpha }{\mu }(\nabla \cdot \vec{A}+i\frac{\omega }{c^2}\varphi )(\nabla \cdot {\vec{A}}^{\ast }-i\frac{\omega }{c^2}{\varphi }^{\ast })-\rho {\varphi }^{\ast }-{\rho }^{\ast }\varphi +\vec{A}\cdot {\vec{J}}^{\ast }+{\vec{A}}^{\ast }\cdot \vec{J}] \end{gathered}$$

$$\begin{gathered} \mathcal{L}={ϵ}_r(\nabla \tilde{\varphi }+i{k}_0\vec{A})\cdot (\nabla {\tilde{\varphi }}^{\ast }-i{k}_0{\vec{A}}^{\ast })-\frac{1}{{\mu }_r}(\nabla \times \vec{A})\cdot (\nabla \times {\vec{A}}^{\ast }) \\ -\frac{\alpha }{{\mu }_r}(\nabla \cdot \vec{A}+i{k}_0{n}^2\tilde{\varphi })(\nabla \cdot {\vec{A}}^{\ast }-i{k}_0{n}^2{\tilde{\varphi }}^{\ast }) \\ -\tilde{\rho }{\tilde{\varphi }}^{\ast }-{\tilde{\rho }}^{\ast }\tilde{\varphi }+{\mu }_0\vec{A}\cdot {\vec{J}}^{\ast }+{\mu }_0{\vec{A}}^{\ast }\cdot \vec{J} \end{gathered}$$

$$\begin{gathered} S={\int }_{\Omega }{ϵ}_r(\nabla \tilde{\varphi }+i{k}_0\vec{A})\cdot (\nabla {\tilde{\varphi }}^{\ast }-i{k}_0{\vec{A}}^{\ast })-\frac{1}{{\mu }_r}(\nabla \times \vec{A})\cdot (\nabla \times {\vec{A}}^{\ast }) \\ -\frac{\alpha }{{\mu }_r}(\nabla \cdot \vec{A}+i{k}_0{n}^2\tilde{\varphi })(\nabla \cdot {\vec{A}}^{\ast }-i{k}_0{n}^2{\tilde{\varphi }}^{\ast }) \\ -\tilde{\rho }{\tilde{\varphi }}^{\ast }-{\tilde{\rho }}^{\ast }\tilde{\varphi }+{\mu }_0\vec{A}\cdot {\vec{J}}^{\ast }+{\mu }_0{\vec{A}}^{\ast }\cdot \vec{J}\mathrm{d}V+ \\ \gamma {\int }_{\partial {\Gamma }_{\mathrm{Z}}}\left(\hat{n}\times \left(-\nabla \tilde{\varphi }-i{k}_0\vec{A}\right)\right)\cdot \left(\hat{n}\times \left(-\nabla {\tilde{\varphi }}^{\ast }+i{k}_0{\vec{A}}^{\ast }\right)\right)\mathrm{d}S \end{gathered}$$

$$\begin{gathered} {\delta }_{{\tilde{\varphi }}^{\ast }}S={\int }_{\Omega }{ϵ}_r\left(\nabla \tilde{\varphi }+i{k}_0\vec{A}\right)\nabla (\delta {\tilde{\varphi }}^{\ast })-\tilde{\rho }\delta {\tilde{\varphi }}^{\ast } \\ +\frac{i{k}_0{n}^2\alpha }{{\mu }_r}\left(\nabla \cdot \vec{A}+i{k}_0{n}^2\tilde{\varphi }\right)\delta {\tilde{\varphi }}^{\ast }dV \\ -\gamma {\int }_{{\Gamma }_{\mathrm{Z}}}\left(\hat{n}\times \left(-\nabla \tilde{\varphi }-i{k}_0\vec{A}\right)\right)\cdot \left(\hat{n}\times \nabla \delta {\tilde{\varphi }}^{\ast }\right)dS \end{gathered}$$

$$\begin{gathered} ={ϵ}_r{\int }_{\Omega }\left(i{k}_0(\alpha -1)\nabla \cdot \vec{A}-\alpha n^2{k}_0^2\tilde{\varphi }-{\nabla }^2\tilde{\varphi }-\frac{\tilde{\rho }}{{ϵ}_r}\right)\delta {\varphi }^{\ast }dV \\ +{ϵ}_r{\int }_{\partial \Omega }\left(\nabla \tilde{\varphi }+i{k}_0\vec{A}\right)\delta {\tilde{\varphi }}^{\ast }dS \\ +\gamma {\int }_{{\Gamma }_{\mathrm{Z}}}\left[(-\nabla \tilde{\varphi }-i{k}_0\vec{A})-\left((-\nabla \tilde{\varphi }-i{k}_0\vec{A})\cdot \hat{n}\right)\hat{n}\right]\nabla (\delta {\tilde{\varphi }}^{\ast })dS \end{gathered}$$

$$\begin{gathered} {\delta }_{{\tilde{\varphi }}^{\ast }}S={ϵ}_r{\int }_{\Omega }\left(i{k}_0(\alpha -1)\nabla \cdot \vec{A}-\alpha n^2{k}_0^2\tilde{\varphi }-{\nabla }^2\tilde{\varphi }-\frac{\tilde{\rho }}{{ϵ}_r}\right)\delta {\varphi }^{\ast }dV \\ +{ϵ}_r{\int }_{\partial \Omega }\left(\nabla \tilde{\varphi }+i{k}_0\vec{A}\right)\delta {\tilde{\varphi }}^{\ast }dS+\gamma {\int }_{{\Gamma }_{\mathrm{Z}}}\nabla \cdot ({\widetilde{\stackrel{⃗}{E}}}_t^{\prime }\delta {\tilde{\varphi }}^{\ast })-\nabla \cdot ({\widetilde{\stackrel{⃗}{E}}}_t^{\prime })\delta {\tilde{\varphi }}^{\ast }dS \end{gathered}$$

$$\alpha n^2{k}_0^2\varphi +{\nabla }^2\varphi -i{k}_0(\alpha -1)\nabla \cdot \vec{A}=-\frac{\rho }{ϵ}$$

$$\left(-\nabla \varphi -i{k}_0\vec{A}\right)\cdot \hat{n}={\vec{E}}_n=0$$

$${ϵ}_r{E}_n+\gamma \nabla \cdot {\vec{E}}_t^{\prime }=0$$

$$\begin{gathered} {\delta }_{{\vec{A}}^{\ast }}S={\int }_{\Omega }-i{k}_0{ϵ}_r(\nabla \tilde{\varphi }+i{k}_0\vec{A}){\overrightarrow{\delta A}}^{\ast }+{\mu }_0\vec{J}\cdot {\overrightarrow{\delta A}}^{\ast } \\ -\frac{\alpha }{{\mu }_r}(\nabla \cdot \vec{A}+i{k}_0{n}^2\tilde{\varphi })\nabla \cdot {\overrightarrow{\delta A}}^{\ast }-\frac{1}{{\mu }_r}(\nabla \times \vec{A})\cdot (\nabla \times {\overrightarrow{\delta A}}^{\ast })dV \\ +i{k}_0\gamma {\int }_{{\Gamma }_{\mathrm{Z}}}[\hat{n}\times (-\nabla \tilde{\varphi }-i{k}_0\vec{A})]\cdot (\hat{n}\times {\overrightarrow{\delta A}}^{\ast })dS \end{gathered}$$

$$\begin{gathered} ={\int }_{\Omega }[-i{k}_0{ϵ}_r\nabla \tilde{\varphi }+{ϵ}_r{k}_0^2\vec{A}+{\mu }_0\vec{J}+\frac{\alpha }{{\mu }_r}\nabla (\nabla \cdot \vec{A}+i{k}_0{n}^2\tilde{\varphi }) \\ -\frac{1}{{\mu }_r}\nabla \times (\nabla \times \vec{A})]\cdot {\overrightarrow{\delta A}}^{\ast }dV-\frac{1}{{\mu }_r}{\int }_{\partial \Omega }\left({\overrightarrow{\delta A}}^{\ast }\times (\nabla \times \vec{A})\right)\cdot \hat{n}dS \\ +i{k}_0\gamma {\int }_{{\Gamma }_{\mathrm{Z}}}[(-\nabla \tilde{\varphi }-i{k}_0\vec{A})-\left((-\nabla \tilde{\varphi }-i{k}_0\vec{A})\cdot \hat{n}\right)\hat{n}]\cdot {\overrightarrow{\delta A}}^{\ast }dS \\ -\frac{\alpha }{{\mu }_r}{\int }_{\partial \Omega }(\nabla \cdot \vec{A}+i{k}_0{n}^2\tilde{\varphi }){\overrightarrow{\delta A}}^{\ast }\cdot \hat{n}dS \end{gathered}$$

$$\begin{gathered} ={\int }_{\Omega }[i{k}_0{n}^2(\alpha -1)\nabla \tilde{\varphi }+{ϵ}_r{k}_0^2\vec{A} \\ +\frac{1}{{\mu }_r}(\alpha -1)\nabla (\nabla \cdot \vec{A})+\frac{1}{{\mu }_r}{\nabla }^2\vec{A}]\cdot {\overrightarrow{\delta A}}^{\ast }dV \\ +\frac{1}{{\mu }_r}{\int }_{\partial \Omega }(\hat{n}\times (\nabla \times \vec{A})-\alpha (\nabla \cdot \vec{A}+i{k}_0\tilde{\varphi })\hat{n})\cdot {\overrightarrow{\delta A}}^{\ast }dS+ \\ i{k}_0\gamma {\int }_{{\Gamma }_{\mathrm{Z}}}{\vec{E}}_t^{\prime }\cdot {\overrightarrow{\delta A}}^{\ast }dS \end{gathered}$$

$$\alpha n^2{k}_0^2\varphi +{\nabla }^2\varphi -i{k}_0(\alpha -1)\nabla \cdot \vec{A}=-\frac{\rho }{ϵ}$$

$$n^2{k}_0^2\vec{A}+{\nabla }^2\vec{A}+(\alpha -1)\nabla (\nabla \cdot \vec{A}+i\frac{\omega }{c^2}\varphi )=-\mu \vec{J}$$

$$\left(-\nabla \varphi -i{c}_0{k}_0\vec{A}\right)\cdot \hat{n}={\vec{E}}_n=0$$

$$\hat{n}\times (\nabla \times \vec{A})-\alpha \left(\nabla \cdot \vec{A}+i\frac{\omega }{c^2}\varphi \right)\hat{n}=0$$

$${ϵ}_r{E}_n+\gamma \nabla \cdot {\vec{E}}_t^{\prime }=0$$

Full text

I Introduction

As physicists and engineers seek to model increasingly complex electromagnetic systems, from radio-frequency power sources to integrated photonics, the need for efficient and robust full-wave, first-principles numerical field solvers is growing. Finite element (FE) methods, which solve partial differential equations over a discretized problem domain (often a spatial mesh), are a natural solution enjoying widespread use in disciplines from fluid dynamics to structural mechanics.

In electromagnetic problems, the variational formulation driving the FE method presents unique challenges which impede the computational efficiency and accuracy of existing solvers, however. For time harmonic problems where the fields oscillate at angular frequency $\omega$, $\vec{\mathcal{E}}(\vec{r},t)=\mathrm{Re}[\vec{E}(\vec{r})e^{i\omega t}]$, this variational expression is given by eq. 1. Here, $\mu_{r}$ and $\epsilon_{r}$ are the relative permeability and permitivitty, $k_{0}$ is the wave number in free space, $Z_{0}$ is the intrinsic impedance of free space and $\vec{J}$ is the current density. For a more detailed treatment see, for example, ref. Jin, 2015.

[TABLE]

The primary challenge consists in enforcing the divergence constraint associated with Gauss’s law, $\nabla\cdot\vec{E}=0$ (in the source-free case) while ensuring adequate freedom in the basis functions used to expand the approximate solution so as to be able to model discontinuities in the fields Albanese and Rubinacci (1997); Cortes Garcia et al. (2018). Alternatively, when working with the magnetic vector potential, $\vec{A}$, the divergence requirement is necessary to enforce the Coulomb gauge, $\nabla\cdot\vec{A}=0$. These two requirements conflict in standard nodal element based finite element (FE) methods, employed successfully in other fields such as fluid mechanics and structural mechanics.

In problems with charge and current density, $\rho$ and $\vec{J}$, there is the added dilemma of how best to satisfy both the Ampere-Maxwell equation and Gauss’s law. In particle-in-cell codes, used in plasma and accelerator physics, this is critical as discrete charge conservation is not automatically guaranteed. Correction schemes must be applied to either the field calculation or the source deposition to bound the error in $\nabla\cdot\vec{E}$ and avoid any resulting numerical instabilities. Marder (1987); Langdon (1992); CIARLET and LABRUNIE (2009); BARTHELMÉ et al. (2007); Ciarlet et al. (2014) For stationary and low-frequency or broadband problems, such as in electro-quasistatics and integrated circuit design, mixed finite-element solutions combined with tree-cotree splitting of the mesh and/or Lagrange multipliers are commonly applied to account for the contributions from $\rho$ and $\vec{J}$ in the static limit separately Li et al. (2016); Hiptmair et al. (2008); BARTHELMÉ et al. (2007).

Historically, the application of differential geometry in three dimensional (3D) Euclidean space has sucessfully resolved the first issue, providing a theoretical motivation for the use of edge elements in the expansion of the fields. Despite this success, to the authors’ knowledge, the full four dimensional (4D) covariant framework has not been investigated within the context of numerical electromagnetism. We demonstrate herein that such an extension is not only more naturally suited to numerical analysis, the field theory Lagrangian providing a variational form directly applicable to the finite element method, but resolves entirely the two significant issues discussed above.

We propose a formulation for the finite element solution of electromagnetic systems based on the classical field theory Lagrangian with a gauge fixing term adapted from quantum electrodynamics. This is given in the framework of differential geometry by eq. 2, where we adopt the standard terminology given by introductory texts such as ref. Stoker, 1989. As will be introduced in greater detail shortly, the first two terms constitute the classical Lagrangian and the final term, multiplied by the scalar $\frac{1}{2\xi}$, is the gauge fixing addition.

[TABLE]

Here $\mu$ is the magnetic permeability, $\mathbf{A}=A_{\nu}\mathrm{d}x^{\nu}$ is the 4D differential 1-form, $A^{\nu}=(\frac{\phi}{c},\vec{A})$ is the four-potential comprised of the electrostatic potential $\phi$ and magnetic vector potential $\vec{A}$, and $\mathbf{J}$ is the electric current 3-form. In this paper we use the metric signature ($+---$). Instead of taking the variation of eq. 1 to obtain $\vec{E}$, we propose taking the variation of the action $S[\mathbf{A}]=\int\mathcal{L}$. This formulation fully accounts for the charge density, $\rho$ in addition to the current $\vec{J}$ through $\mathbf{J}$ and facilitates a return to the widely used nodal FE framework.

The paper is organized as follows: section II provides the background for this work, including a more in depth discussion of the edge elements and their benefits and challenges. We then move from 3D Euclidean space into 4D Minkowski space, presenting the classical field theory Lagrangian in section III. We demonstrate how this formulation is related to existing $\vec{A}-\phi$ approaches obtained by substituting the potentials into eq. 1 yet differs in a few critical points which ensures a fully 4D formulation and enables the use of nodal elements instead of edge elements. Finally, in section IV and V we introduce a proof of concept implementation and demonstrate the validity of the formulation, benchmarking it against a state of the art edge element field solver in terms of accuracy, numerical robustness and flexibility.

II Background

We motivate this idea by considering how the challenges inherent to eq. 1 are currently resolved. Enforcing the divergence constraint was initially addressed by adding a regularization term of the form $s(\nabla\cdot\vec{E})^{2}$ to eq. 1, with limited successHara et al. (1983); Rahman and Davies (1984); Demkowicz (2017); Adam et al. (1997). While the regularization term eliminates spurious non-solenoidal modes in the solution spectrum, the regularized formulation fails to converge to the correct solution for problem geometries with sharp or re-entrant corners. The explanation for this failing was only recently understood: when nodal basis functions are used in conjunction with the regularization term, the approximate solution space spanned is overly restrictive on non-convex domains Costabel and Dauge (2003). When singularities in the field exist such as at sharp corners, it can be shown that the missing subspace consists of the gradients of solutions to Laplace’s equation on the same domainCostabel and Dauge (2000); Assous et al. (1999). Instead of converging to the correct solution with field singularities, the solution obtained will be the projection of the correct fields on the smooth approximate solution space.

Two approaches exist to resolve this issue. One can either supplement the nodal basis functions with additional singular or non-conforming functions Dhia et al. (1999); Assous et al. (1998); Duan et al. (2014), or relax the regularization term near the singularity Costabel and Dauge (2002); Ciarlet Jr., Patrick et al. (2010); Otin (2010). The former requires computing the coupling between the nodal basis and the additional singular functions and is challenging to extend to three dimensions while the latter is a compromise between enforcing the divergence constraint over the problem domain and not completely restricting the subspace spanned by gradients.

A more robust solution arises by formulating electromagnetism in the language of differential geometry. In 3D Euclidean space, $\vec{E}$ (and $\vec{A}$ in the Coulomb gauge) are both differential 1-forms which should be expanded not on nodes but on edges. This led to the development and widespread adoption of the “edge elements” for electromagnetic problems, as developed separately by Whitney Whitney (1957) and Nédélec Nedelec (1980). By their construction, only tangential continuity is imposed at the faces between elements, resolving the issue of modeling field discontinuities at interfaces and boundaries. However, while the edge elements are divergence-free locally, the discontinuity in the normal field at element interfaces allows for solutions that are not divergence-free globally. The space spanned by edge elements divides into the desired space of weakly divergence-free fields and its co-domain, the kernel of the curl operator (purely gradient functions in topologically trivial domains) Bossavit (1988); Bossavit and Mayergoyz (1998). There are methods to extract the gradient field so as to span only the divergence-free fields, such as the tree-cotree method Albanese and Rubinacci (1988); Trapp et al. (2002); Wang et al. (2010); Manges and Cendes (1997). However, choosing an optimal tree is challenging and poor conditioning of the resulting linear system is a common issue Ticar et al. (2002); Golias and Tsiboukis (1994); Preis et al. (1992); Ahagon and Kameari (2017).

III Covariant 4D Field Theory Formulation

III.1 Insight from Differential Geometry

Compared to the 3D formulation, electromagnetic theory is encoded much more succinctly by differential geometry in 4D Minkowski space, where deeper underlying structure is made explicit, such as gauge and Lorenz invariance. For Lorenzian manifolds, the equations for the electromagnetic field tensor, $\mathbf{F}=\mathrm{d}\mathbf{A}$, are given by eqs. 3-4. As defined in section I, $\mathbf{A}=A_{\nu}\mathrm{d}x^{\nu}$ is the four dimensional (4D) differential 1-form, $A^{\nu}=(\frac{\phi}{c},\vec{A})$ is the four-potential, and $\mathbf{J}=-\rho\mathrm{d}x\wedge\mathrm{d}y\wedge\mathrm{d}z+j_{x}\mathrm{d}t\wedge\mathrm{d}y\wedge\mathrm{d}z+j_{y}\mathrm{d}t\wedge\mathrm{d}z\wedge\mathrm{d}x+j_{z}\mathrm{d}t\wedge\mathrm{d}x\wedge\mathrm{d}y$ is the current 3-form.

$\mathrm{d}\mathbf{F}=0$

(3)

$\mathrm{d}\star\mathbf{F}=\mathbf{J}$

(4)

The classical field theory Lagrangian which encodes these equations, written in terms of $\mathbf{A}$, is given by eq. 5.

[TABLE]

In the covariant treatment it is natural to work with $\mathbf{A}$ as opposed to $\mathbf{E}$ and $\mathbf{B}$, the 3D components of the two-form $\mathbf{F}$. In this case, eq. 3 is automatically satisfied as $\mathrm{d}^{2}=0$. Additionally, the charge density, $\rho$, which does not enter into the conventional variational expression for the fields given by eq. 1, is accounted for in eq. 5 through $\mathbf{J}$.

Applying the variational formulation given by the action of eq. 5 in the finite element method, we first consider the appropriate elements over which to expand the solution, $\mathbf{A}$. As in 3D Euclidean space, where there is a duality between 1-forms and edges, in a 4D mesh the 1-form $\mathbf{A}$ should be expanded using edge elements. This is an intriguing idea to pursue for transient numerical analysis, where one could envision using this formulation on a 4D mesh. However, in the case of time harmonic problems, the focus of this paper, we pursue a different course. The time dimension of the mesh is collapsed and edges of the 4D mesh become points in a 3D mesh. As such, nodal elements should be employed to expand $\mathbf{A}$ rather than edge elements.

III.2 Gauge Invariance and Ill-Conditioned Systems

The linear systems resulting from employing eq. 5 in the finite element method are, unfortunately, highly ill-conditioned. This should come as no surprise when considering that unlike the fields, the four potential is not uniquely defined. The Lagrangian is invariant to gauge transformations of the form $\mathbf{A}\rightarrow\mathbf{A}+d\psi$ where $\psi$ is a scalar 0-form. To resolve this issue, we apply a solution used to address a similar challenge in quantum electrodynamics (QED). In the field discretization in QED, the resulting symbolic matrices are singular due to gauge invariance. One approach to overcome this issue is through the addition of a gauge fixing term to the Lagrangian Schwartz (2014):

[TABLE]

The resulting action integral for time harmonic problems, integrated over the time dimension already, is given by eq. 7. Here we have expanded $\mathbf{A}$ into the conventional three-plus-one notation ($\vec{A}$ + $\phi$) for ease of comparison to existing formulations and to expose some implementation challenges which will be discussed in section IV. Note also that as we are now working in the frequency domain, $\phi$ and $\vec{A}$ as given below are complex quantities. The strong problem corresponding to eq. 7, obtained by taking the variation of this action integral, is derived in the appendix.

[TABLE]

In QED, the gauge fixing term imposes different gauges depending on the value of $\xi$. In our classical context, any $\xi\neq 0$ imposes the Lorenz gauge with residual gauge freedom, $\mathbf{A}\rightarrow\mathbf{A}+\mathrm{d}\psi$ for $\psi$ satisfying the wave equation $\nabla^{2}\psi+k_{0}^{2}\psi=0$. By setting the components of $\mathbf{A}$ or their derivatives explicitly on the boundary, $\psi$ is forced to zero on the boundary, and hence everywhere, and $\mathbf{A}$ will be unique.

Tempting as this may be, doing so in the most straight-forward way (setting $\vec{A}_{t}$ =0 and $\phi=0$ on the boundary for a perfect electric conductor, or $\vec{A}_{n}$ =0 and $\nabla_{n}\phi=0$ for a perfect magnetic conductor) in fact decouples $\phi$ from $\vec{A}$. This reduces the problem to two separate wave equations: eq. 8 amenable to nodal elements and eq. 9 which reduces to the curl-curl equation for $\vec{A}$ and must be solved via edge elements.

[TABLE]

On this note, we have since discovered the work of Boyse and Paulsen, who had started to develop a nodal element based formulation by substituting the potentials into Maxwell’s equations and applying the Lorenz gauge. Paulsen et al. (1992); Boyse et al. (1992). Their weak formulation is missing some of the coupling terms that appear in the Lagrangian formulation with gauge fixing and the resulting strong form of the problem has no direct coupling between $\vec{A}$ and $\phi$ in the volume. However, the main issue in their original formulation is precisely the application of the boundary condition scheme suggested above, decoupling $\vec{A}$ and $\phi$. While not mentioned in the original works, a nodal implementation of the formulation will fail on non-convex domains, as noted in ref. Jog and Nandy, 2014. Nonetheless, it is perfectly acceptable to adopt such an approach, as long as a mixed formulation is employed, with $\vec{A}$ expanded by the edge elements.

This decoupled, mixed formulation approach has been adopted by some low-frequency solvers. In these implementations the Lorenz gauge simplifies to the usual divergence free condition on $\vec{A}$ and the two components are coupled only in the sense that their excitations are related through the charge continuity equation, $\nabla\cdot\vec{J}=i\omega\rho$.Li et al. (2016) In this sense, it is not a fully four dimensional solution.

Instead, we allow the residual gauge freedom to persist by imposing boundary conditions solely through surface integrals instead of explicit Dirichlet boundary conditions. This is not a significant drawback, as surface integrals are commonly used to implement impedance or absorbing boundary conditions in any case. The natural boundary condition in our formulation corresponds to a perfect magnetic boundary. An impedance boundary condition can be imposed through the additional surface integral given by eq. 10. A perfect electric boundary is imposed in the limit where the conductance, directly proportional to $\gamma$ is large. For further discussion of the boundary condition imposed by eq. 10 and for a definition of $\gamma$ in the context of our implementation, please see the appendix.

[TABLE]

In contrast to the mixed formulations, the resulting solution is not a simple superposition of independent solutions for $\vec{A}$ and $\phi$, but a self-consistent solution for $\mathbf{A}$ in its entirety, as is demonstrated in section V. This, in conjunction with the gauge fixing term which regularizes the problem, is what enables the use of nodal elements in the four-potential formulation compared to existing $\vec{A}-\phi$ formulations Li et al. (2016); Duan et al. (2018); Baumanns et al. (2013); Amrouche et al. (1998); Wang et al. (2010). Boyse and Paulsen arrived at a similar conclusion with their Maxwell based $\vec{A}+\phi$ formulation, implementing an impedance boundary condition in a later paper which they then demonstrated working on a 2D wedge geometry. Boyse and Paulsen (1997)

Finally, we conclude this section by introducing a new coefficient for the gauge fixing term:

[TABLE]

This is both for the sake of brevity as the coefficient of the gauge fixing terms is referred to often in the implementation and results section, and also as we are interested in plotting solution properties as a function of $\alpha$ near 0.

III.3 Gauss’ Law

While the ability to use nodal elements is a nice benefit, the primary motivation for our adoption of the Lagrangian formulation is the fully general treatment of the source terms provided. In the appendix, we provide the full derivation of the strong form corresponding to eq. 7 or, in the language of variational calculus, the Euler-Lagrange equations resulting from the variation of the Lagrangian. Equations 12 and 13 give the resulting equations imposed in the volume (there are also surface terms which are provided in the appendix).

[TABLE]

Equation 12 is the result of the variation with respect to $\phi$ while eq. 13 arises through the variation with respect to $\vec{A}$. With the Lorenz gauge implicitly imposed through the gauge fixing term, regardless of the residual gauge these equations reduce to Gauss’ law and Ampere’s law in this gauge. Thus, unlike in eq. 1, both equations are independently satisfied by the solution which minimizes eq. 7, even in the case where discrete charge conservation is not guaranteed.

Not only does the solution explicitly satisfy Gauss’ law, but the use of the nodal elements means it can do so element-wise as well as globally. This is in contrast to the lowest order edge elements, which are divergence-free within each element: any non-zero divergence in the fields arises only through discontinuities in the normal component of the field between elements of the mesh. Furthermore, there is significant flexibility offered by the fact that both $\rho$ and $\vec{J}$ can be used to drive the fields.

Both of these features are beneficial in modeling problems with significant space charge, whether for low-frequency applications where $\rho$ becomes important in the static limit or, as in our motivation for pursuing this approach, in the modeling of high frequency power sources where time harmonic components of the space charge contribute strongly to the fields even at high frequencies.

IV Computational Implementation

As a proof of concept, we have implemented this formulation for 2.5D azimuthally symmetric fields, solving on a 2D mesh and accounting for the azimuthal dependance of the fields, of the form $e^{im\theta}$, a-priori. A few unique challenges arise in the implementation of the Lagrangian finite element formulation. The issue of enforcing boundary conditions through surface integrals instead of having the option of explicitly setting Dirichlet boundary conditions was discussed in section III.2. The other significant difference relative to the curl-curl formulation is the presence of terms linear in $k_{0}$.

In finite element electromagnetic analysis, there are two types of problems which are of interest: eigenmode analysis and driven problems. In driven problems, the driving frequency is known so that only the fields need to be computed. In the eigenmode analysis, the resonant frequencies (eigenvalues) and corresponding four-potentials (eigenvectors) are calculated. The discretized Lagrangian is composed of three finite element matrices, $\mathrm{\mathbf{M}},\mathrm{\mathbf{C}},\mathrm{\mathbf{K}}$ and the resulting matrix equation is a generalized quadratic eigenvalue problem (QEP) where we solve for $\tilde{k_{0}}$ , the approximate resonant frequency, and $\mathbf{a}$, the coefficients of the approximate solution over the discretized space.

[TABLE]

It is the coupling between $\vec{A}$ and $\phi$, appearing in the matrix $\mathrm{\mathbf{C}}$, that results in a quadratic eigenvalue problem instead of the regular generalized eigenvalue problem of the curl-curl equation. QEPs are common in finite element problems, for example in modeling damped structural resonances and a significant body of work exists on the topic, including a comprehensive review paper Tisseur and Meerbergen (2001).

We implemented the sparse non-linear eigenvalue solver, NLFEAST Gavin et al. (2018), a contour integral based solver where we constructed the kernel specifically for our QEP. Our implementation has proven robust, agreeing with the direct solver for small problem sizes where a comparison was possible, and scalable up to matrix sizes on the order of 1E6 (we did not test beyond this as for a 2D mesh, this is a very dense mesh). The condition number of the eigenvalues, as defined in ref. Tisseur and Meerbergen, 2001, are reasonable and uncorrelated to problem size. The conditioning does depend weakly on the gauge fixing term and the need for the gauge fixing term becomes immediately clear from the singularity in the condition number when it is not included, as will be shown in sec. V.

The driven problem employs the same $\mathrm{\mathbf{M}},\mathrm{\mathbf{C}}$ and $\mathrm{\mathbf{K}}$ matrices, but $\tilde{k_{0}}$ is set by the frequency of the driving source terms and $\mathbf{a}$ is determined by solving the resulting linear system.

[TABLE]

For this, we use the Intel Math Kernel Libraries, and in particular, the PARDISO solver. We have tried direct and iterative solvers and found both to be equally effective for the moderate problem sizes we have been working with so far.

Future work will look to scale the implementation to 3D meshes and thus much large matrices. Here, a more advanced solver and the employment of a preconditionner will likely be beneficial. We expect NLFEAST or a similar contour integral solver will still be the optimal choice for the eigenmode analysis.

V Numerical Results and Benchmarking

We have benchmarked the Lagrangian formulation with respect to the edge element curl-curl formulation over a broad range of examples. The following subsections focus particularly on numerical results demonstrating the accuracy, robustness and flexibility of this formulation.

The examples shown are azimuthally symmetric, solved in a cylindrical coordinate system $(r,\theta,z)$ with an azimuthal dependence of the form $e^{im\theta}$, as given by eq. 16. As this dependence is known a-priori these modes can be solved on a 2D mesh with $\theta$ out of plane. In the following, all figures of mode profiles are thus cross sectional views of the full structure in the ($z,r$) plane.

[TABLE]

For monopole modes ($m=0$) the fields split into modes which can be represented by $A_{\theta}$ alone (transverse electric or TE), or as a combination of $A_{z},A_{r}$, and $\phi$ (transverse magnetic or TM). The TE modes are not susceptible to the challenges discussed previously and are already often solved using nodal basis functions so we focus only on TM modes for $m=0$. To fully prove the suitability of the Lagrangian formulation, particularly in regards to eventual extension to a full 3D finite element implementation, we also demonstrate some examples of dipole ($m=1$) and quadrupole modes ($m=2$). In this case, the problem is fully four-dimensional and all components of the four-potential couple to each other.

For comparison we used COMSOL, a commercially available multi-physics finite element software which includes an edge-element electromagnetic field solver. It is capable of solving axisymmetric in-plane fields on a 2D mesh, allowing for a comparison with our computational implementation in terms of accuracy and problem size. We refer the reader to the COMSOL user manual for exact implementation detailsCOM . While our implementation uses nodal Lagrange elements and COMSOL is using edge elements, in both cases the elements are second order.

Finally, in the convergence plots that follow, we define the error as follows: for the frequency, the error is computed as $\Delta f=\frac{f-f_{\mathrm{theor}}}{f_{\mathrm{theor}}}$ if a theoretical solution exists, or for the ridge waveguide, by the frequency of the problem on a finer mesh than those plotted. For the fields, we calculate the $S_{0}$ (Sobolev Zero) norm of the field error over the entire problem domain: $\Delta E=\left|\frac{E-E_{\mathrm{theor}}}{E_{\mathrm{theor}}}\right|_{S0}$. To keep the plots legible, instead of showing the error of all six fields, we use the averaged error norm, $\Delta F=\frac{1}{6}(|\Delta E_{r}|_{S0}+|\Delta E_{z}|_{S0}+...)$.

V.1 Accuracy

The cylindrical pillbox cavity is a good initial test case as results can be compared to the analytical solution. Figure 1 shows the cross sectional problem geometry and the mode profiles for the TM011 mode with a perfect magnetic boundary condition on the walls. Only the components of $E_{r}$ and the full vector field plot are shown for the sake of brevity. Two distinct solutions for the four-potential are shown, however, corresponding to different values of the gauge fixing coefficient $\alpha$ defined in eq. 11. Changing $\alpha$ numerically perturbs the system, producing a solution with a different residual gauge, $\psi$. Nonetheless, the resonant frequencies and fields calculated from the different solutions for the four-potential correspond to the same mode.

An interesting consequence of calculating different $\mathbf{A}$ for the same mode is that the numerical error is different in each case, as demonstrated by fig. 2. As the mesh is refined, all solutions converge to the same frequency and fields. Plotting this convergence, now for only a few values of $\alpha$, fig. 3 demonstrates similar convergence characteristics for both the nodal and edge elements. The slopes of the linear fits match that predicted from theory for second order elements, converging as O($h^{4}$) where $h$ is the maximum mesh edge length.

A possible downside in solving for the four-potential is that the desired end results are the electromagnetic fields, not the potentials. As the fields are obtained through derivatives of the potential, they are not expected to converge at the same rate as the solution itself. This is also an issue with the curl-curl formulation, as the magnetic field must be calculated from the solution for the electric field or vice-versa. There are methods to resolve or mitigate this issue, for example the superconvergent patch recovery technique often employed to compute stress in structural mechanics problems.Zienkiewicz and Zhu (1992) However, here we take the simplest approach, taking derivatives of the second order basis functions to compute the fields at the mesh nodes, which still produces comparable results with those computed by COMSOL. The convergence of the fields for the mode in fig. 1 are given by fig. 4.

V.2 Flexibility

Moving on to problems where nodal element based solvers using the conventional curl-curl equation fail, fig. 5 plots the solution for a notched pillbox cavity with a perfect electric boundary. There is a singularity in the fields on the corner which conventional nodal field solvers cannot resolve, converging to the incorrect solution even as the mesh is refined. While the solution for the four-potential is continuous, the discontinuity in the fields at the notch is fully captured by $\nabla\phi$ in the Lagrangian formulation, as demonstrated in the figure. In this case, the frequency computed by COMSOL and the four-potential formulation is 120.0 MHz. If instead, we set $\phi=0$ either on the boundary or the entire volume, decoupling the four-potential, we find that instead the frequency computed is 127.1 MHz. The field profile for this (incorrect) mode is shown in fig. 6, now with no singularity at the re-entrant corner.

In fig. 8, the convergence of the frequency is plotted as a function of the number of degrees of freedom solved for. Instead of plotting as a function of mesh size, where we do not expect to obtain a theoretical rate of convergence due to the singularity in any case, we plot as a function of problem size to illustrate another perhaps counter-intuitive result. The absolute accuracy relative to problem size is comparable despite the additional degree of freedom used in the four-potential formulation. This is because edge elements require roughly twice as many degrees of freedom as nodal elements for the same convergence orderMur (1994); Jin (2015). This is in part due to the additional degrees of freedom per mesh element and in part because there are many more edges than nodes in a mesh.

Spherical cavities can also be modeled in 2.5D, presenting another example with re-entrant corners (when approximated as a polygon) but one with a theoretical solution to which we can compare. As we have not yet implemented curvilinear or isoparametric elements, the convergence rate in this case is dominated by the extent to which the curved boundary is approximated by a polygon. The results shown are for the TM331 mode but note that in this case, TM refers to transverse magnetic with respect to $\rho=\sqrt{r^{2}+z^{2}}$, the convention for spherical cavities, and not with respect to $\theta$ so all components of the four-potential must be solved for.

Finally, in addition to field singularities due to singular boundaries and re-entrant corners, the four-potential formulation can model discontinuities at material interfaces without the special treatment typically required to accommodate the jump in the normal fieldJin (2015). Figure 13 shows the field profile for a tapered dielectric lined cavity, for example. $\phi$ and $\vec{A}$ are continuous but $\nabla\phi$ captures the discontinuity in the fields due to the change in $\epsilon_{r}$.

V.3 Robustness

There are two particular aspects to robustness that we consider here: numerical conditioning and the question of spurious modes. A rigorous theoretical analysis is beyond the scope of this paper but we have investigated these issues experimentally. For the eigenmode analysis, we use the definition of condition number for a quadratic eigenvalue given in ref. Tisseur and Meerbergen, 2001. Figures 14 and 15 plot the condition number for the quadratic eigenvalue of the notched pillbox with a perfect magnetic and electric boundary condition, respectively. Plotted as a function of $\alpha$, the need for the gauge fixing coefficient is clear from the singularity in the condition number as $\alpha\rightarrow 0$. There is a similar singularity in conditioning at $\alpha=0$ for the linear system in the driven problem.

The condition number is reasonable for the perfect magnetic boundary condition and there is no strong variation with $\alpha$ or $h_{\mathrm{mesh}}$. The same cannot be said for the impedance or (in the limit of large conductance) perfect electric boundary condition. The addition of the surface integral to impose the perfect electric boundary negatively impacts the condition number. We have found the condition number for a given eigenvalue to scale linearly with the conductance, $Y=\frac{1}{Z}$ and as O($h_{\mathrm{mesh}}^{-2}$) when the impedance boundary is applied. The condition number can be mitigated to some extent by refining the mesh on the boundary while maintaining constant mesh in the interior, but future work will focus on resolving this issue more efficiently by modifying the surface integral or its implementation.

Up to this point, we have focused on specific modes to demonstrate the convergence and stability of the Lagrangian formulation. It is equally important to ensure that in addition to obtaining correct modes, the solved spectrum is free of unphysical modes. Figure  16(a) plots the spectrum for the first several solved monopole modes of the spherical cavity. Comparing with the theoretically expected modes, we note the presence of two unexpected modes. These are not spurious modes in the conventional sense, however, but rather are pure gauge modes. As can be seen from the mode profiles, these are valid solutions for the four potential which result in zero field (to within numerical noise). Unlike spurious modes in the nodal curl-curl formulation, these modes converge as the mesh is refined in a similar manner to the expected modes. Most importantly, the number of these pure gauge modes in a given frequency interval does not increase with problem size. In the driven problem, these modes are not excited by sources, as can be seen in fig. 16(b) (note, only modes with Az on axis are excited so not all resonant modes are present in spectrum).

VI Conclusion

Concluding, we have demonstrated a new finite element formulation to solve time harmonic electromagnetic fields. By encoding the physics of electromagnetism in a different mathematical formulation, the Lagrangian formulation does not suffer from the challenges inherent to the conventional curl-curl equation for $\vec{E}$. In contrast to the curl-curl equation, where $\vec{J}$ is the only driving term, our formulation completely accounts for both $\vec{J}$ and $\rho$. Both Gauss’ law and Ampere’s law are satisfied, not just globally but over individual elements, in contrast to the commonly employed Nédélec edge elements. This is of importance in the analysis of beam driven radiation sources, for example, where the contribution to the fields from the space charge, $\rho$, can be significant even at high frequency.

We show through both theory and experimental results that the nodal elements are the correct basis choice for our 4D formulation. Indeed, our implementation demonstrates that the four-potential formulation easily handles field singularities and discontinuities unlike nodal element curl-curl implementations. We have benchmarked a proof of concept implementation against COMSOL, a state of the art edge element solver, showing that comparable performance can be obtained. Currently, our surface integral for imposing an impedance or perfect electric boundary condition produces accurate results for problem sizes up to around $n_{DOF}=10^{5}$. However, the scaling of the condition number with mesh size and conductance needs to be addressed. Finally, we have demonstrated that this approach is not susceptible to spurious modes though pure gauge modes with zero fields do appear in the eigenmode spectrum.

The Lagrangian formulation provides unique opportunities for the numerical analysis of electromagnetic fields. While here we present some initial results confirming the accuracy, flexibility and robustness of this idea, we believe there is much yet to explore, particularly in the time domain. From a practical point of view, the adoption of the four-potential also offers a straightforward solution for those interested in a nodal field solver. This is not only beneficial in terms of the computational efficiency and simplicity of nodal elements, but given the widespread use of nodal elements in fields from structural mechanics to fluid dynamics, allows for a common framework for multi-physics problems.

Acknowledgements.

This project was funded by U.S. Department of Energy under Contract No. DE-AC02-76SF00515.

Appendix A Variation of the Field Theory Lagrangian

The Lagrangian for the electromagnetic four-potential including the gauge fixing term and free sources is given by eq. 17. For simplicity we assume constant permeability and permitivitty, $\mu$ and $\epsilon$ in a given mesh element thus neglecting derivatives of these parameters, however future work could consider extending this to anisotropic heterogeneous materials even within a mesh element. We also work with $\alpha=\frac{1}{2\xi}$ as the coefficient of the gauge fixing term to keep the notation clean. Finally we are using the metric signature ($+---$) in the following.

[TABLE]

Expanding the four-potential in terms of the components, $\mathbf{A}=(\frac{\phi}{c},\vec{A}$) where $c$ is the speed of light in the medium and, as we are working with time harmonic potentials in the frequency domain, replacing derivatives with respect to time by $i\omega$, we obtain eq. 18.

[TABLE]

From the point of view of numerical stability and so as to work with parameters having the same dimensions, it is better to normalize $\phi$ by $c_{0}$, the speed of light in a vacuum, and work with $k_{0}=\frac{\omega}{c_{0}}$ instead of $\omega$. We thus change to the variable $\tilde{\phi}=\frac{\phi}{c_{0}}$. Similarly for the space charge density, let us define $\tilde{\rho}=\frac{\rho}{c_{0}\epsilon_{0}}$. Finally, let us define a normalized field $\tilde{\vec{E}}=\frac{\vec{E}}{c_{0}}=-ik_{0}\vec{A}-\nabla\tilde{\phi}$. The units of $\nabla\tilde{\phi}$ and $k_{0}\vec{A}$ are now both [V][s][m]-2.

This change of variable results in a common factor of $2\mu_{0}$ over all terms in eq 18 other than the source terms. We divide the entire expression by this factor and subsequently ignore it as in the subsequent analysis we are only interested in setting the variation of this to zero. This gives eq. 19 where $\epsilon_{r}$ and $\mu_{r}$ are the relative permitivitty and permeability. The refractive index, $n=\sqrt{\epsilon_{r}\mu_{r}}=c_{0}/c$ now appears in the gauge fixing term as we normalized by $c_{0}$ rather than $c$ (so as to work with the wavenumber in free space).

[TABLE]

The full action integral comprising the variational formulation, including the surface impedance boundary integral, is given in eq. 20. We are interested in taking the variation over the closed volume $\Omega$ with a perfect magnetic boundary on the surface $\Gamma_{\mathrm{PM}}$ and an impedance boundary on the surface $\Gamma_{\mathrm{Z}}$. The factor $\gamma=\frac{1}{i\omega\epsilon_{0}Z_{0}}$, the reason for which will be made clear at the end of this appendix.

[TABLE]

The variation of eq. 20 with respect to $\tilde{\phi}^{*}$ is taken first:

[TABLE]

[TABLE]

Where we applied the divergence theorem and the identity $(\vec{A}\times\vec{B})\cdot(\vec{C}\times\vec{D})=(\vec{A}\cdot\vec{C})(\vec{B}\cdot\vec{D})-(\vec{A}\cdot\vec{D})(\vec{B}\cdot\vec{C})$ to the impedance integral. Let us define the projection of $\tilde{\vec{E}}$ on the tangent surface as $\tilde{\vec{E}}_{t}^{\prime}=\tilde{\vec{E}}-(\tilde{\vec{E}}\cdot\hat{n})\hat{n}$ (note that this is orthogonal to $\hat{n}\times\tilde{\vec{E}}$). Next we apply the vector identity, $\vec{B}\cdot\nabla\alpha=\nabla\cdot(\alpha\vec{B})-\alpha\nabla\cdot\vec{B}$ to the impedance surface integral:

[TABLE]

We can then apply the divergence theorem on a surface to convert the first term in the impedance boundary to a line integral which vanishes over a closed surface. The corresponding set of equations enforced through the variation with respect to $\phi^{*}$ is then given by eqs. 24 -  26 where we have given the final system of equations in terms of the unnormalized potential and fields. There is also a corresponding set for the complex conjugate terms, obtained through the variation with respect to $\tilde{\phi}$.

[TABLE]

We shall reserve the discussion of these terms until after the variation with respect to $\vec{A}^{*}$, which we now take.

[TABLE]

[TABLE]

[TABLE]

The total set of equations that are satisfied when $S$ is minimized (including eqs. 24 -  26 obtained through the variation with respect to $\tilde{\phi}^{*}$) are as follows. Once again, we have substituted the original expressions for $\phi$, $\vec{E}$ and $\rho$ into these equations.

[TABLE]

Once again, there are another set of equations consisting of the complex conjugate terms, corresponding to the variation with respect to $\tilde{\phi}$ and $\vec{A}$. When the Lorenz gauge is enforced, $\nabla\cdot\vec{A}+i\frac{\omega}{c^{2}}\phi=0$ and the boundary conditions imposed by the natural boundary condition are those corresponding to a perfect magnetic conductor.

[TABLE]

Meanwhile, on the surface with the impedance surface integral added, $\Gamma_{\mathrm{Z}}$, the boundary condition imposed by the variation with $\vec{A}$ is as follows:

[TABLE]

From which we determine the relationship between our original $\gamma$ and the complex impedance on the boundary, $Z$:

[TABLE]