Transverse fields in waveguides

$\newcommand{\veps}{\varepsilon} \newcommand{\og}{\omega} \newcommand{\ii}{\imath} \newcommand{\Jcl}{J_{\pi/2}} \renewcommand{\d}{\partial} \DeclareMathOperator{\rot}{rot} \DeclareMathOperator{\grad}{grad} \DeclareMathOperator{\curl}{curl} \DeclareMathOperator{\dive}{div} $

Recall the Maxwell system for the time-harmonic electric ($E$) and magnetic ($H$) fields,

\begin{align} -\ii \og \mu H + \curl E & = 0, \\ \ii \og \veps E + \curl H & = - \sigma E - J^a. \end{align}

It represents six equations for the six unknown electric and magnetic field components. The system is driven by $J^a,$ the externally applied current, and we are given the frequency $\og,$ the electric permittivity $\veps$, the magnetic permeability $\mu$, and the electric conductivity $\sigma$.

Consider electromagnetic propagation through a medium with translational symmetry in one direction. A typical example is an infinite waveguide or an optical fiber whose longitudinal axis represents the direction of its translational symmetry.

In the presence of translational symmetry in a direction, the six-equation Maxwell system decouples into two systems, each of three equations, known as the TE (transverse electric) and the TM (transverse magnetic) systems. In this activity, we derive the TE system for a lossy waveguide, its weak formulation and finite element approximation, and numerically solve it.

The decoupling

The three-dimensional (3D) curl can be written in terms of the two-dimensional (2D) curl when there is translational symmetry. To do so, let $e_x, e_y$ and $e_z$ denote the coordinate unit vectors in $x, y$ and $z$ directions. Recall that the 2D curl is obtained by composing divergence with a rotation -- we defined it in a previous notebook using the operator $\Jcl$ that rotates vectors clockwise by 90 degrees as follows: $$ \curl u = \dive (\Jcl u) = \d_x u_y - \d_y u_x. $$ We will also need the 2D rot operator, which represents rotated 2D gradient of a scalar-valued function $\phi$, defined by
$$ \rot \phi = \Jcl (\grad \phi) = \begin{bmatrix} \d_y \phi \\ -\d_x \phi \end{bmatrix} $$

If a smooth 3D vector field $E = E_x e_x + E_y e_y + E_z e_z$ has translational symmetry in $z$-direction (i.e., if $\d_z E$ vanishes), then the application of the 3D curl reduces to applications of the 2D curl and rot as follows: $$ \curl E = \begin{bmatrix} \rot E_z \\ \curl E_{xy} \end{bmatrix} $$ where $E_{xy} = E_x e_x + E_y e_y$.

The above identity is all we need to witness the decoupling. Replacing every $\curl$ in the 3D Maxwell system using this identity, simple calculations show the following.

If the domain, $J^a, \veps, \sigma, \mu$ all have translational symmetry in the $z$-direction, then the six-equation 3D Maxwell system decouples into the three-equation TM system \begin{align} -\ii \og \mu H_{xy} + \rot E_z & = 0, \\ \ii \og (\veps + \frac{\sigma}{\ii\og}) E_z + \curl H_{xy} & = - J^a_z, \end{align} and the three-equation TE system \begin{align} -\ii \og \mu H_z + \curl E_{xy} & = 0, \\ \ii \og (\veps + \frac{\sigma}{\ii\og}) E_{xy} + \rot H_z & = - J^a_{xy}. \end{align}

The first system is called the Transverse Magnetic (TM) system since it does not contain the longitudinal component of the magnetic field $H_z$, only the transverse magnetic components in $H_{xy}$. Similarly, the second system is called the Transverse Electric (TE) system since it does not contain the longitudinal component of the electric field $E_z$.

Eliminating $H_z$ from the last system, we obtain a system solely for the components of the electric field transverse to the direction of symmetry, namely $E_{xy}$ satisfies $$ \rot\mu^{-1} \curl E_{xy} - (\og^2\veps - \ii \og\sigma) E_{xy} = - \ii \og J_{xy}^a. $$

In the remainder we focus on computing the TE field $E_{xy}$ for a specific example. Henceforth we drop all the $xy$ subscripts.

An infinite cylindrical waveguide

Consider an infinite cylindrical tube-shaped waveguide (such as an optical fiber) placed along the $z$ axis in 3D. The configuration has translational symmetry in the $z$-direction. It has a circular cross section $\Omega$ in the $xy$-plane, enclosing two layers of lossy dielectric materials with positive conductivity, and the entire cross section is enclosed by a perfect conductor.

import ngsolve as ng
from ngsolve import curl, dx, x, y, CF
from ngsolve.webgui import Draw
from netgen.occ import Circle, Glue, OCCGeometry, X, Y

r = 2; r0 = 1.5
c = Circle((0,0), r).Face()
c.edges[0].name='out'
inner = Circle((0,0), r0).Face()
outer = c - inner
outer.faces.name = 'polymer'
inner.faces.name = 'core'
domain = Glue([inner, outer])
g = OCCGeometry(domain, dim=2)
mesh = ng.Mesh(g.GenerateMesh(maxh=0.1))
Draw(domain);

This class activity is to solve for the TE field in such a waveguide, given the following parameters:

• Both layers of materials in the waveguide have the same dielectric properties except for differing conductivity. Both $\mu$ and $\veps$ are constant functions. We may therefore multiply through the previous equation for the TE field by $\mu$.

\begin{align} \rot \curl E_{xy} - (\og^2\veps\mu - \ii \og\mu\sigma) E_{xy} & = - \ii \og\mu J_{xy}^a, && \text { in } \Omega, \end{align}

where $\Omega$ is a disk having a nondimensional radius of 2 units. The material properties after nondimensionalization are given below.

• Since the tube waveguide is enclosed by a perfect conductor, we may use the perfect electric boundary condition on $\d\Omega$. How will you express the perfect electric boundary condition using the transverse field?

• The boundary condition is

\begin{align} E_{xy} \cdot t & = 0, && \text { on } \d\Omega, \end{align}

where $t$ is a unit tangent to the boundary. This is an essential boundary condition, i.e., it is imposed in the finite element space. We incorporate it using the name we gave previously to the outer boundary:

X = ng.HCurl(mesh, order=4, type1=True, complex=True, dirichlet='out')

• After nondimensionalization, we are given that the wavenumber $k^2 = \og^2 \veps \mu$ is given by

$$ k = 15. $$

k = 15

• The tube has a low-loss central region and a absorbing outer polymer layer. This is modeled by setting a piecewise constant function $$s = \og\sigma \mu.$$

• In this example, $$

s = \left{ \begin{aligned} & 0.1 && \text{ if } r < 1.5 \\ & 300 && \text{ otherwise}, \end{aligned} \right. $$ where $r = \sqrt{x^2 + y^2}$ is the radial distance in the transverse plane. Here is a plot of $s$:

s = mesh.MaterialCF({'core': 0.1, 'polymer': 300})
Draw(s, mesh);

• The source is a time-harmonic current pulse (with only transverse components), set centered in the waveguide so that

$$ f = -\ii \og J^a = 10 \ii \exp(-100( x^2 + y^2)) e_x. $$

Here is a plot of the pulse:

f = 10j*ng.exp(-100*( x*x + y*y))
pls = ng.CF((f, 0))
Draw(f.imag, mesh, 'pulse');

Let us compute a finite element approximation to the resulting TE field $E$.

Weak form of the TE system

Multiply the equation $\rot \curl E_{xy} - (\og^2\veps\mu - \ii \og\mu\sigma) E_{xy} = f$ by a test function $v$ and integrate by parts. How will you integrate $\rot$ by parts?

We obtain the following weak formulation for the TE field $E$: $\newcommand{\Ho}{\mathring{H}} \newcommand{\om}{\Omega}$ Find $E_{xy} \in \Ho(\curl, \om)$ satisfying

$$ (\curl E_{xy}, \curl v) - ((k^2 - \ii s) E_{xy}, v) = (f, v) $$ for all $v \in \Ho(\curl, \om).$ Here $(\cdot, \cdot)$ denotes the complex $L^2$ inner product.

Exercise I: Prove that there is a unique solution to the above weak formulation by applying the complex Lax-Milgram theorem whenever $s>0$. Apply the Cea lemma and develop a finite element convergence theory for this example.

Finite element solution

The above weak form directly leads to a finite element method using the $H(\curl)$-conforming Nedelec elements, in much the same way as we have seen for other boundary values problems. Without much ado, we proceed with the simple code.

u,v = X.TnT()
a = ng.BilinearForm(X)
a += (curl(u) * curl(v) - (k**2 - 1j*s) * u*v) * dx
f = ng.LinearForm(X)
f += pls * v * dx
E = ng.GridFunction(X, name='TE field')

with ng.TaskManager():
    a.Assemble()
    f.Assemble()
    
E.vec.data = a.mat.Inverse(X.FreeDofs()) * f.vec

The computed $E$ is a complex vector field. The next plot shows the real part of the computed TE field (zoom in to see the directional arrows).

Draw(E.real, mesh, vectors={'grid_size': 100});

We observe from this solution plot that the TE field appears to decay after it enters the outer absorbing region of higher conductivity. This is the effect of higher conductivity $\sigma$ of the outer layer, which makes it lossy.

Recall that the physical electric field is (not complex, but) given by the real part of the complex time-harmonic field multiplied by $\exp(-\ii\og t)$, i.e., the physical (real) TE field is $$ \mathcal{E}(x, y, t) = \mathrm{Re}( E_{xy}(x, y) e^{-\ii \og t} ). $$

You can get an idea of this harmonic time variation by asking ngsolve for a time-harmonic animation of the computed field. (Note that not everyone applies the time-harmonic ansatz with $e^{-\ii \og t}$, often it is used with $e^{+\ii \og t}$.)

Draw(ng.Conj(E[0]), mesh,  animate_complex=True);

Observe that the solution clearly displays a wave character in space. Of course, we expected wave character in time, due to the time-harmonic assumption, but the solution also shows wave character in space. Why does our TE system give solutions that are wavy in space, even when the source we provided was a pulse with no spatial wave character? The next exercise may help in understanding this.

Exercise II: Suppose a 2D vector field $u = u_x e_x + u_y e_y$ solves $$\rot \curl u - k^2 u = f$$ for some $f$ with $\dive f =0$ and some constant wavenumber $k$. Show that each component of the vector field $U(x, y, t) = u(x, y) e^{-\ii k t}$ solves the wave equation: $\d_{tt} U - \Delta U = f e^{-\ii k t}$.

Exercise III: Consider a configuration of data and 3D domain that is rotationally symmetric about the $z$-axis and examine Maxwell solutions that satisfy $\d_\theta E = \d_\theta H = 0$, where $\theta$ is the azimuthal angle in the cylindrical coordinate system. (The system need not have translational symmetry in the $z$-direction.) For this case, prove an analogue of the decoupling shown in this notebook, formulating the analogues of TE and TM systems in cylindrical coordinates.

---

$\ll$ Table Of Contents
$\ll$ Jay Gopalakrishnan