Introducing the Scattering Element Method

In this paper we introduce the scattering element method (SEM) as a new generic denomination for a special type of electromagnetic field simulation. The SEM is characterized by the fact that every spatial element of the simulation domain is defined by a scattering matrix. Classically, this type of simulation method is known as transmission line matrix (TLM) method, where the unit cell is modeled with transmission lines. In this paper, we consider the two-dimensional case and present an alternative approach for defining/modelling the two-dimensional unit cell. This approach samples plane waves directly in the spatial domain. This wave sampling concept leads to a new unit cell, which is referred as wave sampling matrix (WSM). It turns out, that a SEM grid with this type of unit cell has improved dispersion properties compared to the classical cell. A grid with the WSM can be about 1.5 times coarser to obtain the same 1% phase velocity error. We show how the WSM is embedded in the classical grid. This demonstrates that the WSM cannot be derived with the classical transmission line approach and thus justifies the term SEM as a new generic denomination. Finally, the performance of WSM and “classical” frequency-domain TLM cells are compared in a numerical example determining the cutoff frequencies of a dielectric-loaded rectangular waveguide.


I. INTRODUCTION
The purpose of this paper is to introduce the scattering element method (SEM) as a general generic denomination for electromagnetic field simulations, in which scattering matrices are used to discretize the simulation domain.Every spatial element in those simulations is defined by a scattering matrix.In the scientific community, this type of simulation method is known as transmission line matrix (TLM) method and is defined in the time domain (TD) [1], [2], [3], [4], [5], [6], [7], [8], [9], [10] as well in the frequency domain (FD) [11], [12], [13], [14], [15].In this paper, we consider the two-dimensional case.The classical model of the two-dimensional unit cell of the (TD/FD-) TLM method, consists of a network of four transmission lines.The matrix coefficients of the TLM are the reflection and transmission coefficients of this network.
The reason for introducing the SEM is that we have found a new valid solution for the scattering matrix of the twodimensional unit cell.Our approach is based on the sampling of waves in the spatial domain.The scattering matrix from this approach is referred as wave sampling matrix (WSM) and is a second solution in the frequency domain in addition to the classical FDTLM.The WSM cannot be derived from the classical transmission line approach, but, as we show in this paper, the WSM is embedded in a lattice of FDTLM cells.The main differences between lattices made of FDTLM unit cells and WSM unit cells are the dispersion properties.It turns out that a 1% error of phase velocity can be achieved with a coarser lattice of WSM cells compared to FDTLM cells.To the authors' knowledge, this wave sampling approach and the associated WSM and their properties have not yet been documented in the international literature.Therefore, this paper is not (re)defining the FDTLM, as defined in e.g.[11], [12], [13], [14], [15].With the WSM we have found a new unit cell in the frequency domain that cannot be derived using the classical transmission line model.However, all cells (TD-/FDTLM and WSM) have in common that they are described by scattering FIGURE 1. Overview of the scattering element method (SEM): The SEM is defined in the time and frequency domain and includes two approaches for describing the scattering matrix of the (2D) unit cell.The classical approach models the unit cell with transmission lines.This results in the classical transmission line matrix (TLM) with the symbol S TLM .A transformation of the TLM into the frequency domain results in the FDTLM with the symbol SFDTLM .The tilde indicates that a velocity error has to be corrected before getting the FDTLM with symbol S FDTLM .The second approach of the SEM samples waves in the spatial domain.The result is the new wave sampling matrix (WSM).The symbol of this matrix is S. The FDTLM also emerges from the wave sampling approach.
matrices, and therefore it seems appropriate to unite TLM and WSM under the more general designation Scattering Element Method (SEM).
Fig. 1 shows an overview of the SEM.The SEM is defined in the time domain (TDSEM) as well as in the frequency domain (FDSEM), while the classical TLM method computes the electromagnetic fields in the time domain [6], [7], [8].The scattering matrix of the TLM unit cell is referred as S TLM .This matrix can be transformed into the frequency domain and is referred to as SFDTLM .The tilde indicates that a grid of this kind of unit cell computes electromagnetic waves with an incorrect phase velocity (like the classical TLM S TLM [4, p. 507]).This velocity error can be corrected, which leads to the unit cell S FDTLM .This matrix can be interpreted as the two-dimensional version of [11, p. 19 ff.], [12].Employing the generalized wave sampling approach, an alternative frequency-domain unit cell with its scattering matrix S is introduced in the following, which we call the wave scattering matrix (WSM).We will also show that the velocity-corrected FDTLM unit cell S FDTLM similarly emerges from the wave sampling approach.
This paper is structured as follows.In Section II we discuss all the different types of unit cells of the SEM from Fig. 1 and show in detail the derivation of the WSM.In Section III we present the dispersion properties of FDSEM grids with the different types of unit cells and show numerical results.Section IV shows how the WSM is embedded in a FDSEM grid with unit cells of type S FDTLM .In Section V we show how the wave sampling approach can be generalized and derive the unit cell S FDTLM .In Section VI we show how we use the FDSEM algorithm (with both cells FDTLM and WSM) to calculate the cutoff frequency of a dielectric-slab-loaded waveguide.Section VII concludes this paper.

II. THE UNIT CELLS OF THE SEM
This section elaborates the description of the different types of unit cells of the two-dimensional SEM in Fig. 1.Each of these cells solves the two-dimensional wave equation (in the time or frequency domain) in a small spatial element of size in each direction (x and y).The material inside this spatial element is supposed to be linear, lossless, isotropic and homogeneous.The permittivity and permeability are given by ε = ε 0 ε r and μ = μ 0 μ r .The phase velocity in this material is ( In the SEM, the solution of the wave equation for this small spatial element is represented by a 4 × 4 scattering matrix S x which is generally given by 1 with the incoming (a    In the following subsections, we present these different types of unit cells.First, we repeat the classical TLM with its matrix S TLM .This matrix can be transformed into the frequency domain, leading first to SFDTLM (with velocity error) and in succession to S FDTLM (without velocity error).At the end of this section, we show the derivation of the WSM S.

A. THE CLASSICAL TLM
The classical TLM method is defined in the time domain and computes the impulse response of the simulation area.The model of the unit cell of the TLM is depicted in Fig. 3(a).It consists of two crossing transmission lines.Each line is modeled with lumped elements: C (capacitance per unit length) and L (inductance per unit length).The inductance per unit length is separated into two equal parts L /2 to ensure symmetry.Each single transmission line represents the properties of the material to be simulated (C ≡ ε and L ≡ μ).The phase velocity of the single line is [4, p. 505] As can be seen in Fig. 3(a), the total capacitance per unit length in the center of the cell (at the crossing point) is 2C .Therefore, this unit cell represents a material with twice the dielectric constant of the material to be simulated [4, p. 505]: The maximum phase velocity of this cell (and the TDTLM grid) is ( The scattering matrix of this cell is given by [4, p. 502] The coefficients can be calculated with standard transmission line theory [4, p. 501] and the equivalent transmission line model, which is depicted in Fig. 3(b).The reflection coefficients are and the transmission coefficients τ can be calculated with assuming a lossless material.Many of these cells form the TLM grid.The pulse/wave propagates in this grid.The calculation is done in the time domain.In every time step, the pulse/wave transmits through exactly one TLM cell.An excellent description of this algorithm is given in [3].The velocity error from (5) will be corrected by adjusting the time step by a factor of √ 2 as described in [4, p. 550, 584f.].

B. THE CLASSICAL TLM IN THE FREQUENCY DOMAIN
The transformation of the TLM unit cell into the frequency domain is straightforward.Assuming that a wave enters the unit cell shown in Fig. 3(b) at port 1 (left-hand side), this wave travels with the phase velocity v l through the first transmission line of length /2.At the intersection of the four lines, part of the wave is reflected back in the direction of port 1 with the amplitude of = −1/2 from (7).The other part of the wave passes in equal parts in the direction of ports 2, 3 and 4, each with the amplitude of τ = 1/2 from (8).Each of these reflected and transmitted waves travels again at the phase velocity v l through a transmission line of length /2 before leaving the cell.This results in a relative phase shift between the input signal and the output signals of where f is the frequency of the wave and k is the wave number.The equivalent of the TDTLM scattering matrix from ( 6) is in the frequency domain.The tilde in (9) denotes that the velocity error is not yet compensated.Like the TLM cell from (6), this cell also describes a material with twice the dielectric constant and thus the maximum phase velocity of the cell is also v g from (5).
To correct this error, the phase velocity of the single transmission line must be increased so that This results in a new relative phase shift between the input signal and the output signals of e The scattering matrix of the TLM in the frequency domain with compensated velocity error is therefore The only difference between SFDTLM and S FDTLM is the phase velocity, which is embedded in the phase term.A detailed comparison in terms of the propagation properties of these two cells is given in Section III.

C. THE 2D WAVE SAMPLING MATRIX (WSM)
The derivation of the wave sampling matrix starts with the symmetrical spatial element shown in Fig. 2 and the scattering matrix of a four-port network from (2).
Taking advantage of the symmetry properties of this model, the scattering matrix from (2) can be simplified: There are only three S-parameters R (reflection), T (through) and D (diagonal) that have to be determined.This calculation can be split into two parts.The magnitudes |R|, |T | and |D| can be derived from the condition that S is unitary for a lossless medium and the phases will be calculated with the wave sampling approach.Both will be described in the following subsections.

1) MAGNITUDES OF R, T AND D
The matrix S for a lossless medium is unitary and satisfies the condition [17, p. 183]: with the Hermitian transpose S H and the identity matrix I. (13) can be reduced to by using S from (12).With the general identity of some two complex numbers N and M, N * M + NM * = 2|N||M| cos(ϕ N − ϕ M ) the two ( 15) and ( 16) can be formulated as The ( 14), ( 17) and ( 18) form a nonlinear underdetermined system of equations with six unknowns.This system of equations can be solved with an additional assumption: where A is an arbitrary amplitude.This assumption means, that the spatial element, which is described by S, acts like a Huygens source.The incoming power from one side/port of the element is distributed equally in all directions/ports (cf.[4, p. 498ff.]).
A can be calculated by using (19) in ( 14): From both ( 17) and ( 18) follows (with ( 19)) a fixed phase relation between R and T : Under these conditions, the phase ϕ D does not influence the condition that S is unitary and therefore remains an independent parameter to be determined.With the results from ( 20) and ( 21) the expressions for R, T and D from ( 12) can be simplified to 2) PHASES OF R, T AND D The starting point for the derivation of the phases of R, T and D is illustrated in Fig. 4. We consider a one dimensional wave equation in the x-direction where v = 1 / √ εμ (see also (1)) and f is the frequency.The x-direction has been chosen arbitrarily.With this choice, the phases of R, T and D of the first and the third column of the matrix S can be calculated (y-direction: second and fourth column).In the following discussion, we are initially interested only in the first column of S.
The solution of the wave ( 23) is the wave function ψ. ψ is again a superposition of a forward ψ + and a backward ψ − travelling plane wave in x-direction [18, p. 38], In the following discussion, the absolute amplitudes of these waves are irrelevant.In the scattering formalism, only the relative amplitudes of R, T and D are needed and these amplitudes are 1 /2 (given in ( 22)).
The idea of the wave sampling approach is, that the four ports of the scattering matrix sample the wave function ψ at different "points" in space (port 1 at x 1 , port 2 and 4 at x 2 and port 3 at x 3 ), as depicted in Fig. 4. On every x-coordinate, we get: with x ∈ {x 1 , x 2 , x 3 }.The square brackets indicate that the continuous wave function is now sampled at discrete points.
We then assign the complex wave amplitudes (a and b) of each port to the forward and backward traveling waves at each sample point.For port 1 at x 1 , we get Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
and for port 3 at a 1 and b 3 are complex wave amplitudes pointing in the positive x-direction (cf.Fig. 4).Therefore, they are associated with ψ + .b 1 and a 3 are complex wave amplitudes pointing in the negative x-direction and therefore are associated with ψ − .From ( 26) and ( 27) we can derive the following fractions, which correspond to the phase factors of R and T and have to be equal where = x 3 − x 1 was used.(29) shows the expected lengthdependent phase shift caused by wave propagation through the spatial element of length .However, the phase shift in (28) depends on the absolute coordinate x 1 , which is arbitrary.( 22) solves this problem.It states that both phase shifts of R and T have to be equal, and that defines For port 2 and 4 at x 2 , the assignment of the travelling waves to the outgoing wave amplitudes (b 2 and b 4 ) is ambiguous: it is not clear how the outgoing wave b 2 in y-direction is related to the forward and backward traveling waves ψ + [x 2 ] and ψ − [x 2 ] in x-direction.We get the following two relations for port 2 (port 4 accordingly), which correspond to D and has to be equal to e jϕ D : with /2 = x 2 − x 1 .Analog to (28) and (29) the formula in (31) shows the expected phase shift, while the phase shift in (32) depends again on x 1 .By applying the definition from (30) in (32) this indeterminacy can be corrected: b2 Thus, the assignment of the outgoing wave amplitudes at ports 2 and 4 to the travelling plane waves is also no longer ambiguous.
We can summarize the wave sampling matrix (WSM) as follows, where ( 22), ( 29) and (33) were used.The WSM from (34) is one of two possible solutions for describing a wave propagation with the scattering formalism in the two-dimensional case.The other solution is the FDSEM cell of type S FDTLM (as it was derived above from a transmission-line model).This cell can also be derived by a generalization of the wave sampling approach as described in Section V.

D. MATERIALS IN THE FDSEM
Lossy Materials and material transitions can be simulated in the FDSEM using different wave numbers (k) and an additional transition matrix (T) [4, p. 252] Considering the material transition in Fig. 5 the wave numbers of these two materials are: The parameters ε and μ can also be complex to simulate lossy materials.Depending on which unit cell is used (S or S FDTLM ), the wave numbers of these cells must be adjusted according to the materials in the simulation area, e.g. S 1 with k 1 and S 2 with k 2 according to Fig. 5.An additional matrix must be placed at the material interface, where is the ordinary reflection coefficient between two materials with different wave impedances ( . T is a 2 × 2 matrix and has no spatial extension in the simulation area.A detailed illustration of the material transition and its modelling in the (2D) FDSEM is shown in Fig. 5. Nonlinear and anisotropic materials have not been tested with the new type of unit cell (WSM) but are considered future work.However, such materials are documented for the classical (3D) TLM in [19], [20].

III. DISPERSION RELATIONS
Dispersion relations describe the difference between the phase velocity v (from (1)) of a plane wave (with wave length λ) in the medium to be simulated and the phase velocity v g in the numerical grid.In general, this difference depends on the ratio of /λ and the angle of wave propagation with respect to the axis of the grid.Results of a numerical field simulation can be interpreted as good when the ratio holds at any angle.This condition gets fulfilled better and better the finer the grid of the numerical simulation or the smaller the value of /λ is chosen.
In the following sections, we show the dispersion relations of an FDSEM grid with all three types of unit cells: SFDTLM from (9), S FDTLM from (11) and S from (34).For all grids, we present the extreme cases of the dispersion.These extreme cases occur when a wave propagates either diagonally or axially in the numerical grid [4, p. 507].The phase velocities of these waves are referred to as v g,diag and v g,axial .
Furthermore, we show numerical experiments for each of the three grid cell types.The setup is depicted in Fig. 6.The numerical grid has a dimension of 400 × 400 FDSEM cells (excluding the absorber triangles).The material within this region is vacuum.A point source is placed in the center of the grid.The simulations are performed in the frequency domain.The frequency of the source sweeps, so that /λ has a definition range of 0.05 ≤ /λ ≤ 0.35 ( is fixed at 5 mm) and thus covers the range of practical relevance.The phase of the circular wave fronts is detected at two concentric circular contours (blue lines in Fig. 6).To obtain an absorbing boundary, we placed triangles with an absorbing material (ε absorb = ε 0 (1 − 0.25), μ absorb = μ 0 (1 − 0.25)) at the edges of the simulation area.The material parameters were determined empirically and are selected so that from (36) is zero.The triangles have a height of 60 • and the baseline has a length of 20 • .This is a simplified version of a perfectly matched layer [21], [22].Instead of modifying the material parameter towards the simulation edges, we use absorbing triangles like in an anechoic chamber.The outer ports of the simulation area also have a reflection coefficient equal to zero.

A. DISPERSION RELATIONS OF SFDTLM
The dispersion properties of the classic TLM cell (here SFDTLM ) are already well documented [4], [23].
Pulses/waves, which propagate diagonally to the coordinate axes/transmission lines, are not affected by dispersion.The amplitudes of the transmitted wave components add up in such a way that the original shape is preserved (cf.[4, p. 503 Fig. 4]).The propagation velocity for all values of /λ is always [4, p. 505] with v l from (3).
Pulses/waves which propagate axially (0 • /90 • ) in the grid are affected by dispersion.When describing an axial wave propagation (e.g., x-axis), the SEM grid reduces to a transmission line aligned in x-direction, which has open stubs on both sides.The stubs have a length of /2 and repeat at intervals of along the transmission line [4, p. 504 Fig. 6].The topology is similar to that of a strip line stub filter.A transformation of the open circuit to an arbitrary reactance takes place across the stubs.This impedance transformation depends on both the wavelength λ l = v l/ f inside the transmission line and the length of the transmission line, and thus on the size of the cell.If the length of the stub is λ l/4, the open circuit transforms to a short circuit and no wave propagation takes place.Thus, the SEM grid has a cutoff [4, p. 505].If the length of the stub is less than λ l/4, the transformation results in a capacitive reactance.This now adds to the capacitance per unit length of 2C of the cell and thus influences the propagation velocity of the wave, resulting in the dispersion in the SEM grid [23].The dispersion relation for the axial wave propagation in the SEM grid is given with: Now we apply (3), so that λ l = λ and get [4, p. 507], [9] Fig. 7 shows the results of the numerical simulation.In this simulation, we use the setup of Fig. 6 with FDSEM cells SFDTLM from (9).It can be seen that in the case of diagonal wave propagation (45 • ) the ratio v g/v = 1 / √ 2 holds for every value of /λ.This property is also given in (38).The dispersion curve of the axial wave propagation follows exactly the curve which is given in (40).The cutoff is at a value of /λ = 0.25.In the cutoff region, the simulated dispersion curve for 0 • and 90 • differs from the analytical (green) curve.In this region, the grid cannot simulate a correct wave propagation.A detailed discussion of the wave behavior beyond the cutoff frequency is not in the scope of this paper.The values for /λ > 0.25 are plotted to have a uniform illustration of the dispersion curves of all cells as in Figs. 8 and 9.The smaller the value of /λ becomes, the better the dispersion property of the grid becomes.However, the value of v g/v = 1 / √ 2 is not exceeded (see velocity error in Section II-B).

B. DISPERSION RELATIONS OF S FDTLM
The general dispersion relations of this unit cell are similar to those of SFDTLM .The difference is that with (10) the phase velocity of the transmission lines has been increased by a factor of √ 2. This leads to a diagonal phase velocity in the grid of where (10) was inserted in (38).From (10) it also follows that This result can be inserted in (39) to get the axial phase velocity of the grid [16, p. 76]: The results of the numerical simulation are shown in Fig. 8.We use the setup shown in Fig. 6 with the FDSEM unit cell S FDTLM from (11).It can be seen that the diagonal wave propagation is not affected by dispersion.The phase velocity is constant v g,diag = v for all values of /λ.Axial wave propagation is affected by dispersion.The curve from the numerical simulation follows very well the analytical curve from (43).The cutoff is shifted to /λ = 0.25 √ 2 ≈ 0.35.For values /λ < 0.1, the velocity error is less than 1%.

C. DISPERSION OF THE WSM S
To get an idea of the dispersion relations of a FDSEM grid with unit cells of type S (WSM) from (34), let us first discuss the simulation results.Based on these results, the dispersion relations can be modeled by adapting the formulas of the above sections.
The simulated dispersion curves are depicted in Fig. 9.It is noticeable that the dispersion properties here behave inversely to the previous ones.This means that v g,axial is not affected by dispersion, while v g,diag shows a curve that depends on /λ.Furthermore, the dispersion curves are less curved in the range shown in Fig. 9 and therefore, the 1%-error limit is fulfilled now at values of /λ < 0.15.
A model that describes the dispersion behavior of the FD-SEM grid with these cells is depicted in Fig. 10.It shows in blue an FDSEM cell of type WSM S embedded in a small FDSEM grid with cells of type S FDTLM (orange dashed lines).This small orange grid is rotated by an angle of 45 • with respect to the blue WSM cell.The WSM cell (blue) has a dimension of × .Each of the FDTLM cells (orange) has a dimension of / √ 2 × / √ 2, which follows from Pythagoras.In Section IV, we will give a detailed description of how the WSM can be calculated with the model from Fig. 10.
The dispersion curve for the diagonal wave propagation of an FDSEM grid with cells of type S, can be evaluated by using (42) in (39).Additionally, the new size of the cell has to be considered → / √ 2. This results in [16, p. 77] The dispersion properties of the axial wave propagation of the WSM (blue) correspond to those of the diagonal wave propagation of S FDTLM (orange).These are given in (41).This means that applies here.The cutoff of the FDSEM grid with the WSM S unit cell is at a value of /λ = 0.5, which means that the wave must be sampled through the grid at least twice per wavelength.Thus, the cutoff is exactly at the Nyquist limit.

IV. HOW THE WSM IS EMBEDDED IN A S FDTLM GRID
In this section, we give a detailed mathematical description of the model shown in Fig. 10.We show how the WSM (blue in Fig. 10) can be calculated from the small grid of FDSEM cells of type S FDTLM (orange in Fig. 10).The four ports of the WSM are numbered with Roman numerals, while the eight outer ports of the grid are numbered with Arabic numerals.The first step is to calculate the overall scattering matrix of the four connected grid cells.This can be done with a matrixbased approach [16, p. 33 ff.], [24] or manually with Mason's method [25].The result is a 8 × 8 scattering matrix (49) The symmetry properties of the model shown in Fig. 10 are encoded in the gray scaled colors of (46).For example, the two transmissions from port 1 to port 3 and 4 are equal (D 0 ).The same coefficient is obtained for the transmission from port 1 to port 7 and 8. Now we can use the matrix S 0 to calculate the matrix coefficients of the WSM S. We show the procedure in detail on the reflection coefficient S I,I from port I, which is (according to the numbering of the ports in Fig. 10) the ratio of The only ports that contribute to S I,I are port 1 and 2. The incoming signal a I of port I follows from a 1 and a 2 with: The reference planes of a 1 and a 2 have to be shifted towards the reference plane of port I.This is considered with the prefactor e − jk /4 , which is the phase shift through a single transmission line (see Fig. 10).The outgoing signal b I of port I follows from b 1 and b 2 with: In this case, the reference planes of b 1 and b 2 have to be shifted backwards to the reference plane of port I (positive exponent).Now we can substitute b 1 and b 2 and get (first with standard S-parameter notation): Now we can use (52) to replace (a 1 + a 2 ) in (55) to get After some modifications we get from (56) with R 0 and R 1 from (47) and (48): which is equal to R from the WSM in (34).
The same procedure can be applied to For S II,I only the ports 1 and 2 (for a I ) and the ports 3 and 4 (for b II ) are to be considered.a I is given in (52) and b II follows with: It follows with the S-parameter notation from (46) that Now we can derive where ( 52) and (49) used.S II,I is equal to D from the WSM in (34).
The last S-parameter, that has to be calculated, is For S III,I only the ports 1 and 2 (for a I ) and the ports 5 and 6 (for b III ) are to be considered.a I is given in (52) and b III follows with, We can rewrite (63) with the notation from (46) and get The final result follows with where ( 52) and (50) were used.S III,I is equal to T from the WSM in (34).With the three (57), ( 61) and ( 65) we verify the model, which describes how the WSM is embedded in the FDSEM grid with cells of type S FDTLM .

V. GENERALIZATION OF THE WAVE SAMPLING APPROACH
In this section, we show how to generalize the wave sampling approach from Section II-C.The goal is to show that the FDSEM unit cell of type S FDTLM also emerges from this approach.Just like in Section II-C, a scattering matrix shall also describe a spatial element of size × in a linear, lossless, isotropic and homogeneous material.Instead of using a one-dimensional wave equation such as (23), we consider here the two-dimensional case: The solution can be found in terms of separation of variables and is generally given by [17, p. 20 (three-dimensional case)] ψ x and ψ y are the wave functions in x-and y-direction.These wave functions are a superposition of a forward (+) and a backward (−) travelling plane wave in the respective direction (cf.( 24)): From (68), we get four individual wave functions which are associated to the wave vectors These four wave functions are shown in Fig. 11 along with the spatial element.As the material is assumed isotropic, the propagation coefficients of ψ x and ψ y have the same magnitude.Introducing ψ x and ψ y written similar to ( 24) into (69), it follows that In general, the wave vectors from (70) can be rotated by an angle ϕ against the axes of the SEM coordinate system, because the coordinate system of the general wave solution of (66) shown in (67) does not need to be aligned with the SEM grid axes.This situation is depicted in Fig. 12(a).Now we apply the wave sampling approach to this situation: the four waves (see Fig. 11) propagate at an arbitrary angle ϕ against the axes of the SEM coordinate system (see Fig. 12(a)).The ports of the spatial element (which will be described by a scattering matrix) sample the wave function from (68) or (69).From this, we get the following equation, in which the two terms correspond with the incoming (a) and outgoing (b) complex wave amplitudes of port 1 and 3 of the scattering matrix as follows: x (with some modulation in y-direction).In the same way, we can find an analog equation for port 2 and 4 of the scattering matrix: It should be noted that the plane waves ψ x and ψ y are only dependent on one of the coordinate directions of (66), respectively, i.e. ψ x [x] and ψ y [y].For example, at port 1 this would read ψ x [x 1 ] and ψ y [y 2 ] (see Fig. 11).Equations ( 71) and (72) demonstrate the fundamental concept of the wave sampling approach.They show the direct connection between the complex wave amplitudes a i , b j and the wave functions.Both (71) and ( 72) are evaluated according to the grid positions of the a i , b j defined in Fig. 11.
It turns out that in terms of the wave sampling approach, only two angles ϕ = 0 • (with matrix S) and ϕ = 45 • (with matrix S FDTLM ) are valid solutions for an isotropic material.These solutions are illustrated in Figs.12(b) and 13 for S and, similarly, Figs.12(c) and 11 for S FDTLM .In the following, a more detailed interpretation is given.

A. DISCUSSION OF S
In this approach we sample the Eigen-waves, which propagate parallel to the axes of the coordinate system of the SEM grid, separately.For each direction, we only consider the one-dimensional wave (23)  the y-direction can be moved (virtually) to the middle of the cell, which is illustrated in Fig. 13.The corresponding matrix coefficient is D = 1 2 e − jk /2 from (33).The thus resulting arrangement of sampling points increases the sampling rate, which results in improved dispersion properties, as described in Section III.
For a given simulation domain, this type of matrix solves the two-dimensional wave (66).A direct sampling of the two-dimensional wave function (67) is not necessary for the definition of S. The two directions are combined by the matrix formalism itself.

B. DISCUSSION OF S FDTLM
In this approach the two-dimensional wave function (67) is used to define the matrix S FDTLM .Due to the square geometry of the unit cell, this approach samples the Eigen-waves, which propagate diagonally to the SEM-coordinate system.Both directions are not independent of each other (see.Fig. 11).A virtual shift of sampling points is not possible.
A careful examination of the calculation from (69)-(71) shows, that k 1x = −k 2x must apply.The solution is an angle of ϕ = 45 • (see Fig. 12).The scattering matrix of the unit cell of type S FDTLM can be found by substituting the wave vectors from Fig. 12(c) into the wave functions (69) and calculating the desired fractions of a scattering matrix with (71) and (72).The amplitudes are given with (20).Some details for these calculations can be found in [16, p. 135 ff.].
It is interesting to note that this diagonal wave sampling approach corresponds to the TLM model in Fig. 3, where the transmission lines are oriented in parallel to the coordinate axes.The reason is that the wave at the center of the cell of Fig. 3 is sampled via a section of transmission line (e.g. in x-direction at ports 2 and 4).In contrast, the WSM approach in Fig. 13 samples the wave directly at the center of the cell.

VI. NUMERICAL EXAMPLE
In this section we demonstrate how we use the 2D-FDSEM to calculate the cutoff frequencies of an asymmetrical dielectricslab-loaded rectangular waveguide [26, p. 411   The model is depicted in Fig. 14.
There are two different types of modes that can propagate in this waveguide: longitudinal-section electric (LSE) and longitudinal-section magnetic (LSM) [26, p. 411 ff.].The characteristic equations of the modes are: with and k 0 = 2π f √ ε 0 μ 0 .The mode index n has a definition range of: N 0 for LSE and N for LSM-modes.At the cutoff frequencies the propagation constant k z is zero.To obtain a theoretical reference, we solve the (transcendental) characteristic ( 73) and (74) numerically in a frequency range of 1 GHz ≤ f ≤ 5 GHz with a frequency resolution of f = 10 kHz.The results ( f theoretic ) are listed in Table 1 .
A schematic of the FDSEM model of the dielectric-slabloaded waveguide is illustrated in Fig. 15.The cells S 1/2 have a dimension of 1 mm × 1 mm.We simulate this model with both types of unit cells: WSM and FDTLM.The material and the material interface are modelled as described in Section II-D.The waveguide modes (LSE and LSM) can be selected by the boundary conditions of the waveguide: 1 = 3 = −1, 2 = 4 = 1 for LSE and 1 = 3 = 1, 2 = 4 = −1 for LSM.Due to the two-dimensional approximation of the waveguide, the field variations in the z-direction cannot be simulated and for this k z = 0.The model therefore represents a two-dimensional resonator.The cutoff frequencies of the waveguide are the resonance frequencies.When the frequency sweep passes over a resonance frequency, a phase jump of ±π is observable in the simulated field distribution.These phase jumps can be used to detect the cutoff frequencies in the following way [16, p. 85f.]: 1) Calculate the phase difference at position p in the grid between adjacent frequency steps f n and f n+1 : 2) Add these phase difference over all positions P: and divide by π to get the number of pixels with a phase jump.The maximum peak height equals the number of cells (3200 in our example).The curve of H peak ( f n ) can be interpreted as a mode spectrum.Both effects (low peak height and false resonances) can be explained by the excitation of the resonator.We excite the resonator with a point source at a random coordinate (x, y) with L d < x < L x and 0 < y < L y .If other coordinates are chosen, the individual peaks vary in height.However, the dominant resonances are retained, and false resonances are shifted to another frequency.
The mode spectrum in Fig. 16 shows the results of a coarse frequency scan with f = 10 MHz and the WSM unit cell.After this scan we use a finer frequency resolution ( f = 10 kHz) and simulate the waveguide in a smaller bandwidth around the detected resonances.This procedure is repeated with the FDTLM unit cell.The results for both types of unit cells are listed in Table 1.It can be seen that the relative error with the unit cell of type WSM is lower than the rel.error with the unit cell of type FDSEM.

FIGURE 2 .
FIGURE 2. A two-dimensional SEM unit cell of size × in a linear, lossless, isotropic and homogeneous medium with ε r and μ r .In the SEM, this cell will be described with a scattering matrix S x , where x stands either for TLM or FDTLM.a and b are incoming and outgoing complex wave amplitudes (cf.[16, p. 52 Fig. 3.2]).

FIGURE 3 .
FIGURE 3. Transmission line models of the TLM unit cell: (a) Model with lumped elements (L and C ) (cf.[4, Fig. 5 (b) p. 504]), (b) model with four transmission lines.Each line has a length of /2, a characteristic impedance of Z 0 and a phase velocity of v l .These two models describe a spatial element in of size in x-and y-direction.

FIGURE 5 .
FIGURE 5. Material interface: Physical model and realization in the two-dimensional FDSEM.

FIGURE 6 .
FIGURE 6. Setup for determining the dispersion properties of the FDSEM grid with unit cells of type: SFDTLM , S FDTLM and S. The numerical grid is 400 × 400 cells large (excluding the absorber triangles).Every cell has a dimension of = 5 mm.A point source is placed in the center of the grid.The background material is vacuum and the absorber triangles have a permittivity of ε absorb = ε 0 (1 − 0.25) and a permeability of μ absorb = μ 0 (1 − 0.25).The blue lines indicate where the phase of the circular wave is detected.

FIGURE 10 .√ 2 × / √ 2 .
FIGURE 10.A FDSEM unit cell of type WSM S (in blue) embedded in a grid of four FDSEM unit cells of type S FDTLM (in orange) (cf.[16, p. 77 Fig. 3.13]): The WSM with the four ports I to IV has a size of × .Every cell of the orange grid has a size of / √ 2 × / √ 2. The outer ports of this grid are numbered from 1 to 8.

FIGURE 12 .
FIGURE 12. Wave vectors (cf.[16, Fig. 3.3]): (a) The wave vectors k 1 and k 2 are rotated by an angle of ϕ.This is not a solution in terms of the wave sampling approach.(b) The wave vectors k 1 and k 2 are rotated by an angle of ϕ = 0 • .This describes the FDSEM unit cell of type S. (c) The wave vectors k 1 and k 2 are rotated by an angle of ϕ = 45 • .This describes the FDSEM unit cell of type S FDTLM .

y 3 .
corresp.w. a 1 and b 3 + ψ − x • ψ + y + ψ − y corresp.w. b 1 and a (71) a 1 and b 3 are complex wave amplitudes pointing in +xdirection.The corresponding wave function is ∝ ψ + x (with some modulation in y-direction).a 3 and b 1 are complex wave amplitudes pointing in −x-direction.The corresponding wave function is ∝ ψ −

y corresp. w. a 2 and b 4 + 4 .
ψ + x + ψ − x • ψ − y corresp.w.b 2 and a Both directions (x and y) are handled independently of each other.That means, that the sampling points b 2 and b 4 for the x-direction and b 1 and b 3 for

FIGURE 13 .
FIGURE 13.Detailed model for the derivation of S: (a) Sampling a one-dimensional wave function in x-direction.The sampling points of port 2 and 4 are (virtually) shifted in the middle of the cell, so they are horizontally aligned with port 1 and port 3 at y 2 .(b) Sampling a one-dimensional wave function in y-direction.Here the sampling points of port 1 and 3 are vertically aligned with port 2 and 4 at x 2 .

FIGURE 14 .
FIGURE 14. Cross-section of an asymmetrical dielectric-slab-loaded rectangular waveguide with the dimensions: L x = 80 mm, L y = 40 mm and the dielectric slab with: L d = 20 mm and ε r = 4.82.

Fig. 16
shows the mode spectrum of the example considered here.The peaks of the resonances are clearly visible.The horizontal black line indicates the minimum peak height above which resonances are detected.As can be seen in Fig. 16, there are weak resonances with a lower value of H peak , e.g. at 3.64 GHz (LSM) or 4.29 GHz (LSE) and there are also two false resonances at ≈ 2.2 GHz and ≈ 4.6 GHz visible.

FIGURE 16 .
FIGURE 16.Mode spectrum H peak of the dielectric-slab-loaded waveguide calculated with unit cells of type WSM.The horizontal black line indicates the threshold for the height of the resonance peaks.
S 11 to S 22 are located in the upper-left corner of (46).By using the notation from (46) in (54) we get b ff.].The lossless waveguide has a cross-section of 80 mm × 40 mm