Transverse (2D) electric fields in waveguides

9. Transverse (2D) electric fields in waveguides#

\(\newcommand{\veps}{\varepsilon} \newcommand{\og}{\omega} \newcommand{\ii}{\imath} \newcommand{\Jcl}{J_{\pi/2}} \newcommand{\d}{\partial} \newcommand{\rot}{\mathop{\mathrm{rot}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\curl}{\mathop{\mathrm{curl}}} \newcommand{\dive}{\mathop{\mathrm{div}}} \)

Waveguides are important in modern technology: they confine and guide electromagnetic energy along a prescribed path, enabling applications such as light propagation in fiber-optic communications. Accurately computing the fields inside waveguides is essential for the design of antennas, microwave circuits, photonic devices, particle accelerators, etc.

To model waveguides, let us start by recalling the Maxwell system (from a prior notebook). The Maxwell system, for the time-harmonic electric (\(E\)) and magnetic (\(H\)) fields, reads as follows:

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

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 operating frequency \(\og\), the electric permittivity \(\veps\), the magnetic permeability \(\mu\), and the electric conductivity \(\sigma\).

Consider an infinite waveguide with translational symmetry in one direction and electromagnetic propagation through it. A typical example we have in mind is 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.

9.1. The decoupling: TE & TM systems#

The three-dimensional (3D) curl can be written in terms of the two-dimensional (2D) curl when there is translational symmetry. To show this, 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

\[\begin{split} \rot \phi = \Jcl (\grad \phi) = \begin{bmatrix} \d_y \phi \\ -\d_x \phi \end{bmatrix} \end{split}\]

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, then \(\d_z E\) vanishes and the application of the 3D curl reduces to applications of the 2D curl and rot as follows:

\[\begin{split} \curl E = \begin{bmatrix} \rot E_z \\ \curl E_{xy} \end{bmatrix} \end{split}\]

where \(E_{xy} = E_x e_x + E_y e_y\). This 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{split} \begin{aligned} -\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{aligned} \end{split}\]

and the three-equation TE system

\[\begin{split} \begin{aligned} -\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{aligned} \end{split}\]

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

\[\tag{1} \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.

9.2. An infinite cylindrical waveguide#

Consider an infinite cylindrical tube-shaped waveguide (akin to an optical fiber) placed along the \(z\) axis in 3D. The configuration has translational symmetry in the \(z\)-direction. The waveguide 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. Here is a rendering of the cross section geometry.

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);

Question for discussion: Perfect electric boundary condition on the outer wall of the 3D waveguide, with unit outward normal vector \(n\), is

\[ E \times n = 0. \]

What boundary condition does this imply for the transverse electric field \(E_{xy}\) on the boundary of the 2D cross section?

We now proceed 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{aligned} \rot \curl E_{xy} - (\og^2\veps\mu - \ii \og\mu\sigma) E_{xy} & = - \ii \og\mu J_{xy}^a, && \text { in } \Omega, \end{aligned} \]

    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{aligned} E_{xy} \cdot t & = 0, && \text { on } \d\Omega, \end{aligned} \]

    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,

    \[\begin{split} s = \left\{ \begin{aligned} & 0.1 && \text{ if } r < 1.5 \\ & 300 && \text{ otherwise}, \end{aligned} \right. \end{split}\]

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\).

9.3. Weak form of the TE system#

Multiply equation (1), which in our situation reads as \(\rot \curl E_{xy} - (\og^2\veps\mu - \ii \og\mu\sigma) E_{xy} = f\), by a test function \(v\) and integrate by parts. The previously mentioned boundary condition on \(E_{xy}\) is incorporated as an essential boundary condition into the space.

Then 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

\[\tag{2} (\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.

This 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.

Questions for discussion:

  • How do you integrate \(\rot\) by parts?

  • How will you prove that there is a unique solution to the above weak formulation (2)?

  • Does the finite element approximation converge?

9.4. Finite element solution#

Here is a simple code for computing the Galerkin approximation using Nedelec finite elements.

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 while most modern optics texts applies the time-harmonic ansatz with \(e^{-\ii \og t}\), in quantum physics, one often finds the ansatz used is \(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 a wave character in space.

Question for discussion:

  • 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?

  • Explore the relationship between the time-dependent wave equation and equations like \(\rot \curl u - k^2 u = f\) we have been solving.

9.5. Summary#

We have seen

  • a waveguide with translational symmetry along its longitudinal direction,

  • the symmetry decouples the 6-equation Maxwell system into smaller TE and TM systems,

  • the 2D Maxwell system with \(\rot\) and \(\curl\),

  • a finite element method for the TE system using Nedelec elements,

  • how waves get absorbed as they enter a lossy medium with \(\sigma>0\),

  • time-harmonic solutions that are wavy in space and in time.

Contents