Transformation Optics Combined With Line-Integrals for Fast and Efficient Mode Matching Analysis of Waveguide Devices

Mode matching method is here applied to waveguide fields calculated in a transformation optics framework. The proposed approach allows to obtain modes in terms of an expansion onto analytical basis functions for a waveguide of general cross-section provided that the contour can be expressed as a single-valued function in polar coordinates. Mode matching surface integrals are then easily reconducted to line-integrals on the contour thanks to the analytic expansion. The effectiveness of the proposed technique is validated against full-wave results from a commercial simulator and measurements available in literature. The proposed method shows a substantial speed-up, especially if the contour contains arcs of circumference, with a comparable accuracy. An in-depth study on the numerical convergence of mode matching and on the accuracy of the line integral formulation is also provided.


I. INTRODUCTION
Mode matching (MM) is a well-known full-wave accurate technique to solve scattering problems in devices constituted by waveguides of different cross-sections [1], [2]. The basic brick of a MM solution is the analysis of the step, that is, the junction between two waveguides of different cross-section. The scattering by the step is solved by applying the method of moments, enforcing tangential field continuity. The unknown fields at the step are expanded in terms of the two sets of modes of the two waveguides. Modes can be either analytically known or computed numerically. The so-obtained scattering matrix of the step can be integrated easily in CAD environments for fast design and optimization of complex structures [3], [4].
MM requires the computation of the so-called projection or coupling matrix [5], [6]. The entries of such a matrix can be expressed both by means of surface or line-integrals (SI and LI, respectively) [7], [8]. The computation of the projection matrix is straightforward in case the modes of both waveguides are known analytically, but this is true only for very few cross-section geometries like rectangular, circular, elliptical [9] or triangular [10].
For more generic waveguides, of practical interest for the realization of devices, mode fields must be computed numerically. Consequently, also the integrals for the coupling matrix must be computed numerically. For example, in differential phase shifters, thick irises realized with circular waveguide sections having one cut or two parallel cuts are a common solution to obtain birefringence [11], [12], [13] and numerical evaluation of SIs is relatively time consuming.
In the past, the contour-integral mode-matching technique [14] (CIMM) has been used for iris shapes with only moderate variation of the cross-section from a circle. Thanks to the specific field expansion in terms of cylindrical wave functions [14], double integrals are reduced to contour integrals, but no LI formulation was considered.
The LI approach has been applied in the past when modes were known analytically [15], obtained via boundary-integral resonant-mode expansion (BI-RME) [16] or by finite element method (FEM) [6].
MM is here applied to modes calculated in a transformation optics (TO) framework. TO allows to express each mode in a hollow waveguide of generic cross-section in terms of a linear combination of analytic basis functions in a circular waveguide (CW) filled with an anisotropic medium, once a spatial mapping is defined [17], [18]. While in [17], [18] the focus was limited to computing cut-off wavenumbers and field maps, present work expands [17], [18] by resorting to TO computed modes for the implementation of a MM solution for the waveguide step. To this aim, as a further element of novelty, the LI formulation of MM developed in [7], [8] is exploited.
TO modes are computed numerically as an expansion onto a set of full-domain analytically-known basis functions. This can be at a premium over completely numerical approaches (FEM) [19], which are still suitable but heavier computationally. Furthermore, TO is particularly efficient in case the cross-section is at least partially composed by arcs of concentric circumferences, since for these sectors the CW modes applies. This simplification is not possible in the approach presented in [14]. This geometrical feature is indeed fulfilled by typical waveguide cross-sections used for phase shifters, as the two analyzed in the following.
Accelerating the filling of the projection matrix is important since MM converges slowly [5], [20], [21] and hence a large number of basis functions, i.e., mode fields, and consequently integrals, are necessary to reach convergence.
It is also shown in the following that LI can lead to unacceptable large errors if the modes appearing in the integrals have eigenvalues of similar value. In this case it is advisable to drop LIs and resort to SIs as outlined in [6]. A threshold to determine when SIs should be used will be discussed in the following.
The article is organized as follows. MM and its LI formulation are reviewed and described in Section II, thus providing a thorough overview of well-known MM theory. In Section III the proposed method is applied to modes evaluated exploiting a TO approach. Results are then presented in Section IV for both a single iris in a CW and a CW polarizer. Conclusions are drawn in Section V.

II. MODE MATCHING AND ITS LINE-INTEGRAL FORMULATION
Let us consider a hollow waveguide with a uniform crosssection along z-axis. The filling medium is vacuum, characterized by a permeability μ 0 , and a permittivity ε 0 . The generic transverse electric and magnetic fields (e jωt time-dependence implied) can be written in terms of modes as [22]: In (1) and (2) e i (h i ) is the transverse electric (magnetic) field of mode i, γ i the pertinent propagation factor and a i (b i ) the complex coefficient of the mode propagating in positive (negative) z-direction. The propagation constant γ i is given by: where k 0 = 2π f /c is the wave-number, f the frequency, k c,i is the cut-off wavenumber of mode i, f c,i = c k c,i /(2π ) the respective cut-off frequency, and c the speed of light in vacuum. Transverse electric and magnetic fields can be expressed for transverse electric (TE) modes as [7]: and for transverse magnetic (TM) modes as: In (4) and (5), ϕ i and ψ i are the frequency-independent Hertz-type potentials, TE Y i = γ i /(jωμ 0 ) and TM Y i = (jωε 0 )/γ i the modal admittances, C i and D i two arbitrary normalization constants and i z the unit vector in z-direction. Hertz-type potential ϕ i (ψ i ) is solution of the 2D Helmholtz equation obtained imposing homogeneous Neumann-type (Dirichlet-type) boundary conditions along the waveguide contour.
To model the step, a further index is introduced top right of the symbol: either (w) for the wider waveguide or (s) for the smaller waveguide. Let us place the origin of a Cartesian reference at the step, with (w) extending towards z < 0 and (s) extending towards z ≥ 0 (Fig. 1).
The cross-sections and contours of the two waveguides are (w,s) and L (w,s) , respectively. In the following the problem will be restricted to the case in which the cross-section of waveguide (s) is contained in the cross-section of waveguide (w), that is, (s) ⊂ (w) . This is at no loss of generality since if (s) ⊂ (w) , the step can be seen as two cascading steps, one from (s) to a smaller waveguide of zero length and section (s) ∩ (w) followed by a waveguide of section (w) .
Normalization constants C i and D i in (4) and (5) are chosen so that, for both family of modes, it is: MM returns the generalized scattering matrix (GSM) of the single waveguide step by enforcing the continuity of the transverse fields at the discontinuity. At z = 0 (1) and (2), specialized in each waveguide, become: where summations are truncated to the first N (w,s) modes, so as to comprise all modes in waveguide (w, s) whose cutoff frequency is below a fixed value f SPM [15]. The choice of f SPM will be matter of a parametric analysis in the following.
The MM standard implementation is [1]: 1) enforce continuity of the transverse fields on (s) and vanishing of the transverse electric field in waveguide (w) on (w) \ (s) . This leads to the electric-field continuity equation (EFCE) and the magnetic-field continuity equation (MFCE); 2) project the continuity equations onto a suitable set of weights. Following [8], the electric field modes of the larger waveguide is considered for projecting the EFCE and the magnetic field modes of the smaller waveguide for projecting the MFCE. The following system of algebraic equations is obtained: . . ] are column vectors of dimension N (w) ; a (s) and b (s) , defined in an analogous way, are column vectors of dimension N (s) ; Y (w) (size N (w) × N (w) ) and Z (s) (size N (s) × N (s) ) are diagonal matrices containing the wave admittances Y (w) i and the wave impedances Z (s) i of the modes for the two waveguides; P = [P i j ] is the projection or coupling matrix. Entries of the projection matrix are Since transverse electric fields in (4) and (5) are frequency independent, entries of the projection matrix (10) are also frequency independent. Thanks to (6a) and the properties of the projection [23] it is: and the equality holds in case an infinite number of modes for waveguide (w) is considered. Then, the voltage-wave GSMŜ is: being: where I (w) and I (s) are the identity matrices of dimensions N (w) × N (w) and N (s) × N (s) , respectively. From (12) the power-wave GSM can be obtained as: This will be simply referred to as GSM in the following. Once the GSM for a single step is computed, the GSM for more complex devices involving several step discontinuities can be easily obtained considering the delay matrices for straight uniform waveguides and cascading the GSMs [1].
The most computationally expensive part is the computation of the projection matrix entries (10). It turns out, that only three different types of integrals must be calculated [8], namely: TE e (w) , TE e (s) , TM e (w) , TM e (s) , and TM e (w) , TE e (s) ; since TE e (w) , TM e (s) vanishes thanks to the properties of Hertz-type potentials [8]. It should be noted that coupling integrals involve only the transverse electric field both for EFCE and MFCE since Y (s) The integrals for the LI formulation of MM are [7] TE e (w) i , TE e (s) where k (w,s) c are the cut-off wavenumber of waveguide (w) and (s), respectively, i ν is the normal unit vector in outward direction, and i l = i z × i ν the tangent unit vector to the contour (Fig. 2). Notice that the LIs (15) are in general extended to the entire contour L (s) of waveguide (s). However, if the intersection L (s) ∩ L (w) between the contours of the two waveguides is not empty, the integrals effectively extend to the non-shared contour L (s) \ (L (s) ∩ L (w) ) [7], thus reducing the computational effort.
Relations (15a) and (15b) are only valid when the eigenvalues of the two waveguides are different. When they are equal they cannot be applied, but even if they differ by a small amount, the numerical error propagation leads to results with unacceptable errors [7], [16]. To avoid this, in the present work a threshold th is defined and, if |k (s) − k (w) | < th, then LIs are avoided and full SIs used [6]. These SIs are numerically computed over a triangular mesh of the waveguide cross-section via a zero-order approximation of the integrand. Specific values of the threshold are given in the examples of Section IV.

III. MODE FIELDS FROM TRANSFORMATION OPTICS AND LINE-INTEGRALS IN POLAR COORDINATES
The following two subsections will briefly recall how modes are obtained in a TO framework and then cast the LI MM formulation in polar coordinates, which are the natural coordinate system for TO application.

A. MODE FIELDS FROM TRANSFORMATION OPTICS
In the TO framework, the waveguide contour is expressed in polar coordinates in terms of a continuous single-valued functionρ(φ) of the angular coordinate φ (Fig. 2). The generic cross-section can then be mapped onto a CW of unit radius and anisotropic filling in a new reference (ρ,φ) whereρ = ρ/ρ(φ) andφ = φ. Hertz-type potentials are given by [18]: being c e mn , c o mn , d e mn and d o mn the expansion coefficients associated to basis function mn, and χ mn (χ mn ) the n-th zero of the Bessel function of the first kind and order m, J m (·) (or of its derivative); M, N truncate the summations. Functions r m (ρ), l ι mn (φ), and τ ι with ι = {e, o} are special functions derived in [18] and summarized in Appendix A. Index e (o) is used to indicate functions that are even (odd) with respect to x-axis.
If for a certain interval of the angular coordinate the waveguide contour is an arc of circumference centered in the origin, (16) and (17) reduce to linear combinations of CW modes and the well-known scattering problem for CW is resumed.
In case the waveguide contour exhibits an axis of symmetry, a fictitious magnetic or electric wall can be placed along it, thus halving the domain. Modes are then subdivided as even or odd. In the case of a symmetry along x-axis, coefficients c e mn and d e mn are zeros for odd modes while coefficients c o mn and d o mn are zeros for even modes. If the y-axis is also an axis of symmetry, then this further symmetry can be exploited, reducing the computational domain to one quarter of the full domain.

B. LINE-INTEGRAL FORMULATION IN POLAR COORDINATES
The expressions in polar coordinates for unit vectors i ν and i l in (15) are [18]:  being v a strictly positive function of the angular coordinate.
Integrals (15) have the form: where (20) can be recast in angular coordinates: the radial coordinate being evaluated on the waveguide contour. Specialized expressions for g(ρ(φ), φ)v are given in Appendix A. If, for instance, the waveguide cross-section is symmetric with respect to the x-axis (see Fig. 3), then it isρ(φ) =ρ(−φ) and (19c) implies v(φ) = v(−φ). Consequently (21) vanishes when g(ρ(φ), φ) is an odd function of the angular coordinate and, if it is an even function, it reduces to: Analogous considerations hold for further axes of symmetry. Generally speaking, if N AS is the number of axes of symmetry, the integration interval reduces to [0, π/N AS ] [15]. Summarizing, only modes with symmetries matching the geometrical symmetries exhibit a non-vanishing coupling integral.

IV. RESULTS
In this section the proposed approach is applied to the characterization of two devices: 1) a thick iris in CW. The iris cross-section ( Fig. 3(a)) presents a single cut, hence in the following referred to as 'circle with one cut' (C1C); 2) a full polarizer comprising the cascade of several irises in CW. Each iris cross-section ( Fig. 3(b)) presents two cuts, hence in the following referred to as 'circle with two cuts' (C2C). The numerical results attained with the proposed method are compared with a FEM solution (ANSYS HFSS [24]). For the second device measurements are also available [25]. In both cases, when treating the single step, C1C and C2C are the smaller waveguide (s) while the larger waveguide (w) is the CW, whose radius equals the radius R of the circular sectors of C1C and C2C. The modes of C1C and C2C are evaluated with TO while those of CW are analytical.
The expressions for the contours of C1C and C2C in polar coordinates are similar and fillets are taken into account [26]. Geometry C1C is characterized by one axis of symmetry (x, hence N AS = 1) while C2C by two orthogonal axes of symmetry (x and y, hence N AS = 2). The definition of the contour, restricted to its minimum domain, is, in both cases:  The parameters defining the cross-section are: R the outer radius of the waveguide,x the distance from the center of the circumference and the cut(s), r b the fillet radius. The other parameters in (24) are depicted in Fig. 3.

A. CIRCLE WITH A CUT
Geometrical parameters for the CW-C1C case are: R = 10 mm for both CW and C1C,x = 5 √ 3 mm ≈ 8.66 mm, and r b = 1 mm; the iris thickness is 5 mm and the distance from the ports to the iris on both sides is 40 mm (Fig. 4). The dimensions of C1C are the same as in [18]. The threshold to switch from LIs to SIs is set to th = 60 m −1 .
The device is analyzed in the 8 − 12 GHz bandwidth, central frequency being f 0 = 10 GHz. For the CW, the cut-off frequency of the fundamental TE 11 and of the first higher mode are 8.78 GHz, and 11.47 GHz, respectively. The fundamental mode of the CW is indicated as TE c 11 or TE s 11 , according to its polarization, respectively inductive or capacitive. The iris causes in the C1C section a splitting for the impinging TE 11 according to its polarization and the cut-off frequency too splits assuming two different values: 8.67 GHz, 9.11 GHz. The first higher cut-off frequency does not split and is 11.69 GHz.
Firstly, numerical convergence of MM is studied. Ratio f SPM / f 0 for the mode selection criterion varies from 4 to 12 with a step of 1 and summations in (16) and (17) are truncated to M = N = 5, 10, 15, 20, 25. We define a normalized truncation error as where Q is the quantity whose convergence is analyzed and Q ref is the reference, that is, the value assumed by   The corresponding numbers of modes considered in CW and C1C are then N (w) = 32, 142, 313 and N (s) = 32, 136, 303. Vertical grey lines marks the cut-off frequencies: solid ones for CW and dashed for C1C. Curves for f SPM / f 0 = 8, 12 are almost superimposed. Fig. 7 shows the magnitude and phase of the reflection and transmission coefficients for both polarizations of the fundamental CW mode considering f SPM / f 0 = 8 and M = N = 20.
With these values the normalized truncation errors in Fig. 5(a) and (b) are comparable and around 0.5 %. For frequencies greater than 8.78 GHz, transmission coefficient magnitude gets close to 0 dB, and reflection coefficient magnitude decreases steadily with frequency. The good agreement between results from HFSS and MM is apparent. Table 1 shows the computing times for the construction of the coupling matrix. Absolute times considering SIs and the combination of line and surface integrals (LSIs) point out the effective speed-up. Last row shows the ratio between the LSI times and the full SI times. Absolute times are relative to an Intel Core i7 processor working at 4 GHz with 16 GB RAM memory.
In the LSI case, the elapsed time to compute only the required SIs is estimated multiplying the total elapsed time for SIs times the fraction of integrals that are actually required. Indeed, due to the series expansion for TO modes, it is not possible to compute only some entries of the coupling matrix without calculating the SIs for CW modes and all the basis functions of the TO expansion.
From Table 1 it is apparent that the computation of LSIs is around an order of magnitude faster than the computation of the corresponding SIs. Hence the usefulness of applying the LI formulation to compute the entries of the projection matrix is proved. Furthermore, the computation of the LSIs is still competitive in case the ratio f SPM / f 0 is increased.

B. POLARIZER
To better appreciate the advantage of the proposed MM approach over a full-wave FEM simulation, the device presented in [25] is analyzed. The device exploits eight C2C irises, symmetric with respect to the center of the device. Geometrical details are in [25] and are reported in Fig. 8 for clarification. Distance between ports and first irises is not indicated in [25] and is considered equal to 20 mm in present analysis. In our approach, for a more realistic modelling [18], a fillet with a radius 1/10 of the radius of the CW is also considered in the C2C cross-sections.

TABLE 2. Elapsed Times of HFSS and MM for the Analysis of the Polarizer in [25] and Speed-Up of MM With Respect to HFSS in the Last Row
The frequency bandwidth is 27 − 37 GHz and center frequency is f 0 = 32 GHz. Summations (16) and (17) are truncated at M = N = 25. This value is greater than that of the first example since C2C cross section differs more from a circular one than C1C, hence more bases are necessary to achieve a similar accuracy [17]. Notwithstanding the increase in the number of basis functions, computation time does not increase since thanks to symmetries only a quarter of the cross-section is analyzed in this case. Threshold for SIs in case of TE-TE and TM-TM couplings can be estimated by scaling the threshold from the first example. Being the radii of CWs 10 mm for first example and 4.69 mm for the second one, the threshold is th = (10/4.69) 60 m −1 ≈ 130 m −1 . Fig. 9 shows the results of the analysis for f SPM / f 0 = 7, 8, 9, compared with both ANSYS HFSS and measurements reported in [25]. The results for f SPM = 9 match measurements, except for some minor deviations in the magnitude of the reflection coefficients below −40 dB. These differences in the working bandwidth are due to manufacture tolerances [25] and to measurement noise. The magnitude of the reflection coefficient for TE c 11 inductive polarization converges faster than that for TE s 11 capacitive polarization. In all cases MM matches measurements better than the full-wave simulations by HFSS, as stated also in [11].
Computing times for HFSS and MM, and the speed-up of MM in comparison with HFSS are shown in Table 2. Settings for HFSS are: phase convergence at 32 GHz less than 0.5 • and discrete frequency sweep with 401 points. Notice that MM is about two orders of magnitude faster than HFSS even in the case f SPM / f 0 = 9.

V. CONCLUSION
An implementation of MM based on TO computed modes is presented. The proposed technique allows for a LI based formulation of MM leading to substantial computing time reduction, about one order of magnitude with comparison to SI based MM and about two orders of magnitude with respect to a FEM approach.
Present method can be applied to waveguides of general cross-section provided that the waveguide contour is expressed as a single-valued function of the angular coordinate.
Numerical convergence of the proposed technique is very good and its accuracy is verified against both full-wave results from a commercial software and measurements available in literature. The excellent agreement proves the effectiveness of the method.

A. EXPLICIT EXPRESSION OF TO FUNCTIONS
Expressions for r m (ρ) and l ι mn (φ) in (17) and their derivatives are [18]: where, here and in the following, prime sign indicates the derivative with respect to the argument,ρ = ρ/ρ(φ), and κ = 20 as in [18]. From now on the φ-dependence ofρ(φ) will be omitted. Expressions for , u ι , and τ ι in (26) are: As in the main text, waveguide (w) is CW and waveguide (s) is analyzed with the TO approach. CW modes are indexed with p and q. To implement MM, these are indexed with a single index i. Mapping from the single to the double index is indicated with p i , q i . TO modes in waveguide (s) are indexed with index j and in this case no mapping to a double index is necessary. Since analyzed geometries present at least one axis of symmetry, index ι is added to ψ (w,s) (ϕ (w,s) ) to indicate its ι-symmetry. Then the expressions for terms in (15) evaluated along the waveguide contour ρ =ρ(φ) are  are the specialized expressions for g(ρ(φ), φ)v in (21). In (28) and (29) χ pq andχ pq are the q-th zeros of the Bessel function of order p and of its derivative, respectively, and F pq , G pq are the normalization constants for CW modes where the double prime in (30b) indicates the second derivative with respect to the argument. Notice that the correction term r m (ρ) and its derivative r m (ρ), defined in (26a) and (26b) and introduced in [18], do not appear in the expressions of functions (28) and (29). Moreover, ifρ(φ) = R the terms (28a), (29a), and (29b) vanish, proving that LIs are to be computed only over the non-shared contour L (s) \ (L (s) ∩ L (w) ).