A mixed finite-element formulation for the modal analysis of electromagnetic waveguides featuring improved low-frequency resolution of transmission line modes

Abstract. This paper presents an improved finite-element formulation for axially uniform electromagnetic waveguides. It allows for both dielectric and conduction losses and covers the entire range from optics down to the static limit. Propagation coefficients of small magnitude, particularly those of transmission line modes in the low-frequency regime, are computed much more accurately than with previous approaches.


Introduction
Since electromagnetic waveguides constitute generalized transmission line systems, they are of fundamental importance for the design of electromagnetic devices.Because of complex geometries and inhomogeneous materials, analytical solutions for the modal field patterns and propagation coefficients are often difficult to obtain or even unavailable.Numerical field simulation methods provide a powerful remedy.In particular the finite-element (FE) method stands out for its high convergence rates and great flexibility in modeling geometry and materials.
Early approaches suffered from the occurrence of spurious modes, but with the advent of H(curl) conforming basis functions these problems were overcome.A great variety of FE formulations have been proposed: Some are in terms of the electric (Lee et al., 1991) or magnetic field intensity (Valor and Zapata, 1995), while others employ a magnetic vector potential and a electric scalar potential (Bardi and Biro, 1991;Polstyanko and Lee, 1995).They all provide accurate results for high frequencies, but very few can handle the low-frequency (LF) (Vardapetyan and Demkowicz, 2003) or even the static case (Lee et al., 2003;Farle et al., 2004).Moreover, they incorporate electric conductivity σ by means of an equivalent imaginary part ε of the electric permittivity (Lee, 1994), which breaks down when the angular frequency ω tends to zero, due to ε = σ ω . (1) All methods above lead to generalized algebraic eigenvalue problems for the square of the propagation coefficient.As will be detailed in Sect.4.2, this strongly amplifies numerical round-off error in propagation coefficients of small magnitude, e.g.waveguide modes close to cut-off and, more importantly, transmission line (TEM or quasi-TEM) modes in the LF regime.
To overcome these limitations, we propose in Sect. 2 the prototype of a mixed-field formulation in terms of the electric field intensity E and the magnetic flux density B. Its distinguishing feature is that the sought eigenvalue is the propagation coefficient itself rather than its square.The suggested approach incorporates conduction losses quite naturally and minimizes round-off errors in propagation coefficients of small magnitude.In Sect. 3 we establish stability in the LF regime, by imposing suitable constraints.
The ability of a single formulation to predict the propagation characteristics of inhomogeneous transmission lines accurately from statics up to microwave frequencies is of great practical importance for mixed-signal analysis, particularly in the context of model-order reduction (Polstyanko et al., 1997;Bertazzi et al., 2002;Schultschik et al., 2008).
Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V.

R. Baltes et al.: Mixed finite-element formulation for electromagnetic waveguides 2 Formulation
We consider an axially uniform waveguide pointing in the ẑ direction.Its boundary , with unit outward normal vector n, is assumed to consist of perfect electric walls E and perfect magnetic walls H . Let be connected and E nonempty.We denote the wavenumber, speed of light, and characteristic impedance, respectively, of free space by k 0 , c 0 , and η 0 .The relative magnetic permeability µ r , the relative electric permittivity ε r , and the electric conductivity σ are assumed to be scalar-valued.Thanks to uniformity along the waveguide axis, the modal fields are given by waves propagating in ±ẑ direction.For the negative ẑ direction, we have the decomposition (Pozar, 2005, p. 93) where subscripts t and z denote transversal and longitudinal components, respectively, and γ is the propagation coefficient.Thus the differential operators ∂ z and ∇ simplify to By applying Eqs. ( 2)-( 5) to Faraday's law and Ampere's law, we arrive at a homogeneous boundary value problem, given by the set of partial differential equations subject to the boundary conditions (13)

Finite-element representation
Let H 1 , H curl , H div , and H 0 denote the Sobolev spaces with respect to the operator in superscript (Boffi et al., 2013, pp. 4).We define the subspaces and denote the corresponding finite-element spaces by The transversal and axial fields E t , E z , B t , B z are discretized by basis functions of lowest order (Zhu and Cangellaris, 2006, pp. 19), for their respective function spaces.Let N N denote the number of free nodes, N E the number of free edges, and N F the number of faces in the mesh.Then Note that Eqs. ( 19) and ( 20) have been scaled by c 0 , for numerical stability.The FE representations of the differential operators grad, curl, and div are given by the incidence matrices G ∈ R N E ×N N , R ∈ R N F ×N E , and D ∈ R N F ×N E .Also, in two dimensions, we have (Zhu and Cangellaris, 2006, pp. 26) Since Eqs. ( 8) and ( 9) are in terms of fields represented on the primal mesh, they are discretized in strong form.In contrast, Eqs. ( 10) and ( 11) are based on the fields H = µ −1 B, D = εE, and J = σ E. These reside on the dual complex and are thus considered in weak form.The resulting generalized algebraic eigenvalue problem reads At this point it can be seen that the eigenvalue is the propagation coefficient γ itself rather than its square.The vector and the matrices in Eq. ( 22) are given by In absence of magnetic (dielectric) losses, we have µ r ∈ R + (ε r ∈ R + ) everywhere, and the corresponding matrices T Bt and T Bz (T Et and T Ez ) are positive definite.In contrast, T σ t and T σ z are positive semi-definite, because σ > 0 holds in the conductive region only.

Low-frequency behavior
For k 0 = 0, the eigenvalue problem of Eq. ( 22) reduces to The matrix A + γ B turns out to be singular, for any value of γ .We will analyze this effect and propose a suitable regularization.

Analysis of instability
We decompose the waveguide cross-section into the lossless region ll , with σ = 0, the lossy subdomain l , with σ > 0, and their common interface 0 : In view of Eq. ( 23), the eigenvalue problem for the static case, Eq. (30), is equivalent to −Rx Et = 0, (33) Consider the vectors with x F ∈ R N F arbitrary and any Plugging Eqs. ( 36) and (37) into Eqs.( 32)-( 35) shows that these are satisfied for any arbitrary value of γ .The vectors form the eigenspace corresponding to γ .To regularize the eigenvalue problem and eliminate the non-physical eigenvectors, we first analyze the vectors u γ .Let E γ and B γ denote the fields represented by Eqs. ( 36) and ( 37), respectively.In view of Eq. ( 38), E γ is nonzero in the lossless region only and, by Eq. ( 4), it is a gradient field.Similarly, it can be shown that the magnetic field strength H corresponding to B γ is a gradient, represented on the dual mesh.Plugging E γ and B γ into Eqs.( 6) and ( 7) and taking the divergence yields It can be seen that, for k 0 = 0, these fields are no longer forced to be source-free.

Stabilization
To stabilize the method, we impose the conditions similarly to (Polstyanko et al., 1997;Hiptmair et al., 2008), by introducing Lagrange multipliers p F and p N .The discrete representations of Eqs. ( 42) and ( 43) read where the matrix ĨN selects the nodes belonging to the lossless regions ll .
Care must be taken if the structure is bounded by electric walls only ( H = ∅).In this case, the transpose of the curl matrix R T , i.e. the gradient operator on the dual mesh, has a nontrivial nullspace of dimension one.It consists of constant fields and is represented by a coefficient vector x c with

Rectangular Waveguide
We consider a rectangular waveguide (RWG) of dimensions 22.86 mm×11.43mm with perfectly conducting walls and homogenous material properties (ǫ r = 4, µ r = 1) in both the lossless case (σ = 0) and in presence of Ohmic losses (σ = 5 S/m).Our main goal is to demonstrate the correct function of the present approach and compare it to the A-V potential formulation of (Farle et al., 2004), a state-of-the-art method that solves for γ 2 and models conductivity via complex per-255 mittivity (1).According to (Pozar, 2005, p. 108), the analytical solution for the propagation coefficient γ mn is given by where a and b denote the width and height of the wave-  where J i is the Jacobian matrix of face i.In view of Eq. ( 23), the corresponding kernel vector of A reads It is easy to verify that the fields of u c satisfy Eqs. ( 32)-( 35) and ( 45).Moreover, Eq. ( 44) implies Thus there exists a zero eigenvalue corresponding to a constant magnetic field strength H pointing in ẑ direction.To exclude the eigenpair (γ = 0, x = u c ), we impose the orthogonality constraint by introducing a Lagrange multiplier p c .Thus, the stabilized eigenvalue problem for the case H = ∅ becomes with I 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 −T Et 0 0 0 0 0 0 0 T Ez 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 All eigenvectors of ( 50) that exhibit non-vanishing Lagrange multipliers, i.e. non-zero p N , p F or p c , are known to span a subspace without physical meaning.This means the relevant subvectors do not solve the original eigenvalue problem (22).Rather than to purge such unwanted eigenvectors from the results, the authors prefer to shift the associated eigenvalues to infinity, so that they do not pollute the spectrum.This is easily accomplished, by setting the columns of B that belong to Lagrange multipliers to zero.The modified matrix B reads and the final form of the eigenvalue problem is given by

Rectangular Waveguide
We consider a rectangular waveguide (RWG) of dimensions 22.86 mm×11.43mm with perfectly conducting walls and homogenous material properties (ǫ r = 4, µ r = 1) in both the 250 lossless case (σ = 0) and in presence of Ohmic losses (σ = 5 S/m).Our main goal is to demonstrate the correct function of the present approach and compare it to the A-V potential formulation of (Farle et al., 2004), a state-of-the-art method that solves for γ 2 and models conductivity via complex per-255 mittivity (1).According to (Pozar, 2005, p. 108), the analytical solution for the propagation coefficient γ mn is given by where a and b denote the width and height of the wave- All eigenvectors of Eq. ( 50) that exhibit non-vanishing Lagrange multipliers, i.e. non-zero p N , p F or p c , are known to span a subspace without physical meaning.This means the relevant subvectors do not solve the original eigenvalue problem Eq. ( 22).Rather than to purge such unwanted eigenvectors from the results, the authors prefer to shift the associated eigenvalues to infinity, so that they do not pollute the spectrum.This is easily accomplished, by setting the columns of B that belong to Lagrange multipliers to zero.The modified matrix B reads and the final form of the eigenvalue problem is given by 4 Numerical examples

Rectangular waveguide
We consider a rectangular waveguide (RWG) of dimensions 22.86 mm×11.43mm with perfectly conducting walls and homogenous material properties ( r = 4, µ r = 1) in both the   [10 −10 , 10 10 ] Hz, and Fig. 3b gives the error in propagation coefficient as a function of the number of unknowns.
In the lossless case, the solutions of both the E-B and A-V methods are in very good agreement with the theory, over the whole frequency range; see Fig. 1.Fig. 3a demonstrates that the convergence rates of both formulations are the same, for f = 0 Hz and f = 5 GHz, respectively.In the lossy case (σ=5 S/m), it can be seen from Fig. 2 that the E B approach works over the entire frequency range whereas the A-V scheme fails to converge below 10 −8 Hz, due to breakdown of (1).Above this threshold, the results of both methods agree very well with analytical results.The numerical noise visible in Fig. 2b for phase coefficients of very small magnitude, Im γ < 10 −12 rad/m, is insignificant because, as seen in Fig. 2a, the corresponding attenuation coefficients are more than 14 orders of magnitude larger, Re γ > 120.Fig. 3b indicates that, within their respective range of validity, both numerical methods exhibit the same rate of convergence, independently of frequency.In case of the E-B scheme, this also holds in the static case, f = 0 Hz.

Shielded Microstrip Line
The purpose of this example is to demonstrate the advantage of formulating the eigenvalue problem in terms of γ rather than γ 2 , with respect to round-off error in propagation coefficients of small magnitude.Again, we compare the E-B and A-V methods, by means of the lossless shielded microstrip line shown on the inset of Fig. 4. In contrast to the RWG of Section 4.1, the present structure supports a quasi-TEM come as a surprise, because the latter formulation was designed to be low-frequency stable (Farle et al., 2004).
To clarify the situation, we present in Table 1 a comparison for not only the quasi-TEM wave but also the first two box modes.It can be seen that only the quasi-TEM mode is 300 affected; the results for box modes are correct, even at 0 Hz.This behavior results from the fact that the static limit of |γ| is zero for the quasi-TEM mode but non-zero for all others.At sufficiently low frequency values, |γ qTEM | becomes so much smaller than all other |γ| values that the eigenvalue solver is 305 unable to resolve γ qTEM to sufficient accuracy, due to numerical noise.Note that the gap in eigenvalue because the E-B formulation solves for γ and the A-V 310 scheme for γ 2 .Fig. 5 presents the eigenvalue ratio as a function of frequency.For both formulations, the eigenvalue solver produces significant round-off error for |λ 1 /λ 2 | < 10 −12 . . . 10 −9 .However, the frequency which this happens at is 10 4 Hz for the A-V method, whereas the proposed E-B 315 scheme produces accurate results down to 10 −3 Hz, thanks to improved eigenvalue ratio in (55).lossless case (σ = 0) and in presence of Ohmic losses (σ = 5 S/m).Our main goal is to demonstrate the correct function of the present approach and compare it to the A−V potential formulation of (Farle et al., 2004), a state-of-the-art method that solves for γ 2 and models conductivity via complex permittivity, Eq. ( 1).According to (Pozar, 2005, p. 108), the analytical solution for the propagation coefficient γ mn is given by where a and b denote the width and height of the waveguide, respectively.Figures 1 and 2 present the dispersion curves of the TE 10 and TE 20 modes in the frequency range [10 −10 , 10 10 ] Hz, and Fig. 3b gives the error in propagation coefficient as a function of the number of unknowns.
In the lossless case, the solutions of both the E − B and A − V methods are in very good agreement with the theory, over the whole frequency range; see Fig. 1. Figure 3a demonstrates that the convergence rates of both formulations are the same, for f = 0 Hz and f = 5 GHz, respectively.
In the lossy case (σ =5 S/m), it can be seen from Fig. 2 that the E − B approach works over the entire frequency range whereas the A − V scheme fails to converge below 10 −8 Hz, due to breakdown of Eq. ( 1).Above this threshold, the results of both methods agree very well with analytical results.The numerical noise visible in Fig. 2b for phase coefficients of very small magnitude, Im γ < 10 −12 rad/m, is insignificant because, as seen in Fig. 2a, the corresponding attenuation coefficients are more than 14 orders of magnitude larger, Re γ > 120. Figure 3b indicates that, within their respective range of validity, both numerical methods exhibit the same rate of convergence, independently of frequency.In case of the E −B scheme, this also holds in the static case, f = 0 Hz.

Shielded microstrip line
The purpose of this example is to demonstrate the advantage of formulating the eigenvalue problem in terms of γ   [10 −10 , 10 10 ] Hz, and Fig. 3b gives the error in propagation coefficient as a function of the number of unknowns.
In the lossless case, the solutions of both the E-B and A-V methods are in very good agreement with the theory, over the whole frequency range; see Fig. 1.Fig. 3a demonstrates that the convergence rates of both formulations are the same, for f = 0 Hz and f = 5 GHz, respectively.In the lossy case (σ=5 S/m), it can be seen from Fig. 2 that the E B approach works over the entire frequency range whereas the A-V scheme fails to converge below 10 −8 Hz, due to breakdown of (1).Above this threshold, the results of both methods agree very well with analytical results.The numerical noise visible in Fig. 2b for phase coefficients of very small magnitude, Im γ < 10 −12 rad/m, is insignificant because, as seen in Fig. 2a, the corresponding attenuation coefficients are more than 14 orders of magnitude larger, Re γ > 120.Fig. 3b indicates that, within their respective range of validity, both numerical methods exhibit the same 280 rate of convergence, independently of frequency.In case of the E-B scheme, this also holds in the static case, f = 0 Hz.

Shielded Microstrip Line
The purpose of this example is to demonstrate the advantage of formulating the eigenvalue problem in terms of γ rather than γ 2 , with respect to round-off error in propagation coefficients of small magnitude.Again, we compare the E-B and A-V methods, by means of the lossless shielded microstrip line shown on the inset of Fig. 4. In contrast to the RWG of Section 4.1, the present structure supports a quasi-TEM come as a surprise, because the latter formulation was designed to be low-frequency stable (Farle et al., 2004).
To clarify the situation, we present in Table 1 a comparison for not only the quasi-TEM wave but also the first two because the E-B formulation solves for γ and the A-V 310 scheme for γ 2 .Fig. 5 presents the eigenvalue ratio as a function of frequency.For both formulations, the eigenvalue solver produces significant round-off error for |λ 1 /λ 2 | < 10 −12 . . . 10 −9 .However, the frequency which this happens at is 10 4 Hz for the A-V method, whereas the proposed E-B 315 scheme produces accurate results down to 10 −3 Hz, thanks to improved eigenvalue ratio in (55).rather than γ 2 , with respect to round-off error in propagation coefficients of small magnitude.Again, we compare the E −B and A−V methods, by means of the lossless shielded microstrip line shown on the inset of Fig. 4. In contrast to the RWG of Sect.4.1, the present structure supports a quasi-TEM mode, the phase coefficient of which is known to depend linearly on frequency in the LF regime.Figure 4 gives results for both numerical methods.While the E −B data exhibit the expected behavior, the phase coefficients produced by the A − V method stagnate for frequencies below 10 4 Hz.This may come as a surprise, because the latter formulation was designed to be low-frequency stable (Farle et al., 2004).
To clarify the situation, we present in Table 1 a comparison for not only the quasi-TEM wave but also the first two box modes.It can be seen that only the quasi-TEM mode is affected; the results for box modes are correct, even at 0 Hz.This behavior results from the fact that the static limit of |γ | is zero for the quasi-TEM mode but non-zero for all others.At sufficiently low frequency values, |γ qTEM | becomes so much smaller than all other |γ | values that the eigenvalue solver is unable to resolve γ qTEM to sufficient accuracy, due to numerical noise.Note that the gap in eigenvalue because the E − B formulation solves for γ and the A − V scheme for γ 2 .Figure 5 presents the eigenvalue ratio as a function of frequency.For both formulations, the eigenvalue solver produces significant round-off error for |λ 1 /λ 2 | < 10 −12 . . . 10 −9 .However, the frequency which this happens at is 10 4 Hz for the A − V method, whereas the proposed E − B scheme produces accurate results down to 10 −3 Hz, thanks to improved eigenvalue ratio in Eq. ( 55).

Critical discussion and conclusions
We have presented a mixed finite-element formulation for electromagnetic waveguides.Stability in the LF regime is achieved by imposing source-free conditions for the electric and magnetic flux densities, with the help of Lagrange multipliers.
Since the sought eigenvalue is the propagation coefficient itself rather than its square, round-off errors in propagation coefficients of small magnitude are significantly reduced.This is of great practical importance for the mixed-signal analysis of transversally inhomogeneous transmission line structures, such as microstrips or coplanar waveguides: For the example of Sect.4.2, the LF range is extended by 7 orders of magnitude, compared to γ 2 approaches.Another feature of the proposed method is that Ohmic losses enter the formulation directly rather than in the form of equivalent complex permittivity.
The advantages of the new method do not come for free, though.Table 2 shows that the number of variables is three times as large as in the A − V case.Nevertheless, the nonzeros in Â + Ĉ just increase by 50 %, and those in B actually decrease.Iteration counts of the eigs solver of MATLAB-R2013a are 50 % higher for the E − B approach.As a result, solution times for the lossy RWG are six times as long.For the microstrip, that ratio is twelve because, in the lossless case, the system matrices of the A − V method become realvalued, whereas those of the E −B scheme remain complex.
260 guide, respectively.Figs. 1 and 2 present the dispersion curves of the TE 10 and TE 20 modes in the frequency range

Fig. 4 :
Fig. 4: Shielded microstrip line: LF behavior of propagation coefficient of quasi-TEM mode.Inset shows cross-section of structure.Dimensions are in mm.

Figure 3 .
Figure 3. Convergence analysis for the TE 10 mode.

Fig. 4 :
Fig. 4: Shielded microstrip line: LF behavior of propagation coefficient of quasi-TEM mode.Inset shows cross-section of structure.Dimensions are in mm.

Figure 4 .
Figure 4. Shielded microstrip line: LF behavior of propagation coefficient of quasi-TEM mode.Inset shows cross-section of structure.Dimensions are in mm.

Table 1 :
Comparison of first three modes of microstrip line

Table 1 :
Comparison of first three modes of microstrip line This behavior results from the fact that the static limit of |γ| is zero for the quasi-TEM mode but non-zero for all others.At sufficiently low frequency values, |γ qTEM | becomes so much smaller than all other |γ| values that the eigenvalue solver is 305 unable to resolve γ qTEM to sufficient accuracy, due to numerical noise.Note that the gap in eigenvalue |λ 1

Table 1 .
Comparison of first three modes of microstrip line.