Modeling of Electromagnetic Fields in Parallel-Plane Structures : A Unified Contour-Integral Approach

A unified reciprocity-based modeling approach for analyzing electromagnetic fields in dispersive parallelplane structures of arbitrary shape is described. It is shown that the use of the reciprocity theorem of the time-convolution type leads to a global contour-integral interaction quantity from which novel both timeand frequency-domain numerical schemes can be arrived at. Applications of the numerical method concerning the time-domain radiated interference and susceptibility of parallel-plane structures are discussed and illustrated on numerical examples.


Introduction
The analysis of high-speed multilayered printed circuit boards (PCBs) can be in principle carried out with the aid of standard full-wave numerical techniques (e.g.[1]).The multi-scale nature of such complex systems along with the ever increasing edge rates of transmitted signals, however, lead to computational efforts that are still beyond the capability of the available hardware within an acceptable computational time.A way out of this limitation is splitting the problem configuration to elementary components that admit the (physics-based) equivalent network representation [2], [3].An integral part of such equivalent networks is the parallel-plane impedance matrix that represents the signal transmission characteristics of a (parallel-plane) cavity in a multilayered PCB.Introducing a unified, reciprocity-based computational methodology for electromagnetic (EM) analysis of such a general parallel-plane structure is hence exactly the main objective of this paper.
Previous works on the subject are largely focused on applications of the classic frequency-domain (FD) contour integral method (CIM) [4,Chapter 3] to analyzing transmission parameters and radiated EM emissions of general power-ground structures [5] and its further developments addressing non-uniform field distributions along the port periphery [6], the inclusion of substrate inhomogeneities [7] and the efficient speed-up techniques [8], for instance.More recently, the contour-integral approach has also been formulated directly in the time domain (TD) [9], [10], thereby introducing a novel TD integral-equation numerical method.This computational technique is formulated using the EM reciprocity theorem of the time-convolution type [11,Sec. 28.2] and is known as the time-domain contour integral method (TD-CIM).It has been recently shown that the TD-CIM is an expedient tool for studying TD EM transmission/reception reciprocity properties of parallel-plane structures [12] and provides an efficient means for evaluating their space-time mutual EM coupling [13] and pulsed EM plane-wave responses [14].In particular, it has been demonstrated that these calculations are extremely fast for rectangular parallelplane structures (e.g.[14,Sec. 4.2]) for which the ray-type TD analytical description is available in closed form [15].
Firstly, the problem configuration that is analyzed is described in the following Sec. 2. A generic, complex-FD reciprocity-based relation from which the CIM can be derived both in real-FD and TD is then introduced in Sec. 3.This general formulation represents the starting point for all numerical schemes presented throughout this paper.Subsequently, the relation between the general formulation and the classic FD-CIM is shown in Sec. 4. In this section, the outcomes of two FD numerical schemes are examined with respect to a FD closed-form analytical description concerning a rectangular parallel-plane structure.Subsequently, Sec. 5 introduces a numerical solution of the TD reciprocity relation.The computational scheme is validated using an alternative (three-dimensional) computational method.This section is also supplemented with the Appendix, where some general features of the typical (Boltzmann-type) relaxation models are discussed.Furthermore, applications of the TD-CIM to the evaluation of pulsed EM radiated emissions and susceptibility of a parallel-plane structure are briefly described.Finally, the introduced methodology is validated by a numerical experiment employing the fundamental property of selfreciprocity.This experiment combines the TD-CIM with an independent numerical tool and is conducted directly in TD.

Problem Description
To localize a spatial point in the problem configuration we employ the orthogonal right-handed reference frame defined via its origin O and the standard basis {i 1 , i 2 , i 3 }.The position vector is then The time coordinate is t.The partial differentiation with respect to x 1 , x 2 and t is denoted by ∂ 1 , ∂ 2 and ∂ t , respectively.The continuous time-convolution is denoted by * t and the time-integration operator is ∂ −1 t .The Heaviside unit step function is H(t).The Dirac delta distribution is denoted by δ(t).
We shall analyze a parallel-plane structure shown in Fig. 1.The structure consists of a dielectric slab placed in between two perfectly electrically conducting (PEC) planes.The thickness of the slab d is assumed to be so small with respect to the spatial support of the excitation pulse, which implies that only x 3 -independent EM field components {E 3 , H 1 , H 2 } = {E 3 , H 1 , H 2 }(x, t) are dominant.The EM field components are excited t = 0. Prior to this instant these fields are zero throughout the problem configuration.The vanishing thickness of the slab allows to define a pulsed voltage V (x, t) = −dE 3 (x, t).The spatial solution domain is then bounded by a (time-invariant) domain Ω extending over the PEC plate.The corresponding bounding contour is denoted by ∂Ω whose outward unit normal vector is ν = ν (x).
The EM properties of the analyzed structure are described by the dielectric relaxation function κ(t) whose generic form is written as where χ = χ(t) is the (causal) relaxation function (see the Appendix), i.e. χ(t) = 0 for all t < 0, and is the (realvalued and positive) electric permittivity.The structure is assumed to be instantaneously-reacting in its magnetic behavior, which implies the (impulsive) form of its magneticpermeability relaxation function µ 0 δ(t).The corresponding EM wave speed is c = ( µ 0 ) 1/2 with (.) 1/2 > 0.
The structure is placed in the linear, homogeneous and isotropic embedding D 0 whose EM properties are described by real-valued and positive constants { 0 , µ 0 }.The EM wave speed in the embedding is c 0 = ( 0 µ 0 ) 1/2 > 0.
The one-sided Laplace transformation of a causal and bounded space-time quantity f (x, t) is defined as with {s ∈ C; Re(s) > 0} such that f does exist.The corresponding inversion is denoted by L −1 .A limiting case is arrived at via s = iω + δ as δ ↓ 0 with {ω ∈ R; ω ≥ 0} being the angular frequency.This special case is of interest within the scope of the (real-)FD analysis.Examples of the latter will be given in Sec. 4.

Problem Formulation
Application of the complex-FD reciprocity theorem of the time-convolution type [11,Sec. 28.4] to the actual field state and the testing field state (denoted by superscript 'B') and to domain Ω leads to the following interaction quantity where χ Ω (x) is the characteristic function χ Ω (x) = {1, 1/2, 0} for x ∈ {Ω, ∂Ω, Ω } (Ω denotes the complement of Ω in R 2 ), Ĵ3 is the electric-current volume density and ∂ Ĵ = Ĥ2 i 1 − Ĥ1 i 2 is the electric-current surface density.In (3), the test wave fields are linearly related to their source via for κ = {1, 2}, where Ĝ∞ (r, s) is the (bounded) fundamental solution of the two-dimensional modified Helmholtz equation in R 2 and r (x|x T ) is the Eucledian distance between the points specified by position vectors x and x T .Note that when specifying a testing-field quantity, x| x S stands for the field and source position vectors, respectively.For a general form of the TD counterpart of Ĝ∞ (r, s) we refer the reader to the Appendix.
To get a boundary-contour relation, the testing surface electric current density is next applied to the periphery of the circuit, Ĵ , where δ(x − x S ) being the Dirac delta distribution operative along x S ∈ ∂Ω.In this way we get the following reciprocity relation for κ = {1, 2}.
Reciprocity relation (6) with equations ( 7) and ( 8) serve as the basis for the numerical methods whose description follows.

Real-FD Numerical Solution
The classic FD-CIM is a well-established numerical technique for the efficient analysis of parallel-plane structures of arbitrary shape [4].The following subsections illustrate that the standard FD-CIM formulation is a special case of the starting complex-FD reciprocity relation (6).
In order to solve the resulting Fredholm integral equation of the second kind numerically, the circuit's rim is first approximated by a set of line segments ∂Ω ∪ N m=1 Ω [m] .Subsequently, upon employing the piecewise-constant expansion, the equation is cast into its matrix form [4, Eq. (3.2)]

Point-matching Solution
In the first step, the unknown voltage distribution along the approximated circuit's periphery is expanded in terms of the rectangular functions Π, i.e.
where v[m] are (yet unknown) the voltage expansion coefficients of vector V .
Upon enforcing the equality in (9) at isolated points located at the centers of the dividing segments denoted by x [m;c] , for m = {1, • • • , N }, we end up with (cf.[4, Eq. (3.26)]) for all S m and for all S = m with γ = 0.57721... being Euler's constant.The latter expression has been found using the small-argument expansion of the Hankel function (see [16, Eq. (9.1.8)]).

Pulse-matching Solution
Instead of applying the point-matching procedure we can start over with the real-FD version of the reciprocitybased relation (6) and associate the corresponding testingsource density with the rectangular function, i.e. we let . This choice in combination with the piecewise-constant expansion (11)  algebraic equations (10) with in which x T (λ for all S m, where x [n] localizes the position of the n-th nodal point.Similarly to the previous section, the diagonal terms are handled analytically.After a few steps of algebra we get for all S = m.Finally it is noted that equations ( 12)-( 13) can be understood as a special case of ( 16)- (17) to which the 1-point Gaussian quadrature (see [16, Eq. (25.4.30)]) is applied.

A Numerical Example
In this section we shall analyze a rectangular parallelplane structure of dimensions L = 0.10 (m) and W = 0.15 (m) (see Fig. 2).
The thickness of the structure is d = 1.50 (mm).The EM properties of the dielectric slab are described by its electric permittivity = 4.50 0 and magnetic permeability µ = µ 0 .The circuit is assumed to show low losses such that the corresponding (complex-valued) wavenumber k can be approximated according to [4, Sec.2.2.1]where the dielectric loss is accounted for via tan(δ) = 0.0045, the skin depth of the conductor is found from δ s = 2/ω µ 0 σ with conductivity σ = 5.80 • 10 7 (S/m).The structure is activated using the excitation vertical port that has its center at x S = (75.0,112.5) (mm) and radius 1.50 (mm).All the calculations that follow are performed in the frequency range {50 ≤ f = ω/2π ≤ 2000} (MHz) at 200 uniformly-spaced frequency points.The circuit's boundary is discretized such that max n (| Ω [n] |) < 0.12 c/ max( f ).
For validation purposes, the input impedance is also evaluated using the eigenfunction-expansion approach via the double-summation formula (see [4,Sec. 2.4.1]),namely with e m = 1 for m = 0 and e m = √ 2 for m 0 and k m = mπ/L, k n = nπ/W with F S mn = F P mn , where where we take W S 1 = W S 2 = 1.0 (mm).In the actual calculations, the summations in (25) are truncated to N = M = 10 3 .The obtained results are summarized in Fig. 3.In Figs.3a and  3b we have shown the point-matching and pulse-matching solutions, respectively, together with the FD response calculated according to the analytical solution (24).As can be observed, the both solutions agree with the reference well.In order to clearly assess the accuracy of the numerical solutions, the absolute error of the calculated input impedance with respect to the referential solution (24) has been plotted in Fig. 3c.Apparently, the calculated error curves attain their peak values at the cavity resonance frequencies.Comparing the point-matching and pulse-matching approaches, the latter solution leads, except for the very high-frequency part of the frequency range, to more accurate results.The difference is most evident at low frequencies.Moreover, it has been observed that doubling the number of the integration from K = 6 to K = 12 does not imply a significant precision increase.Finally, as the pulse-matching solution requires computation of integrals, one should carefully consider whether the improvement with respect to the point-matching solution is worth the effort for practical purposes.

TD Numerical Solution
The first step in solving the TD counterpart of (6) numerically consists of the discretization of the bounding contour ∂Ω ∪ N m=1 Ω [m] and of the time axis T = {t k ∈ R; t k = k t, t > 0, k = 1, ..., NT }.The piecewise-linear space-time expansion of the voltage distribution in terms of the triangular functions T, i.e.
where v [m]  [k] are (yet unknown) the voltage expansion coefficients, along with the piecewise linear spatial expansion of the testing current density, i.e.
allow to cast (the TD counterpart of) ( 6) into the following system of algebraic equations (cf.[9, Eq. ( 14)]) that can be solved in a step-by-step manner for all and 3D-array whose elements are given as [9, Eq. ( 16)] whose elements read [9, Eq. ( 17)] for all S = {1, • • • , N } and all t ∈ T , where |∂S| denotes the arc length of ∂S ⊂ ∂Ω.The time-dependent functions in (30) and (31) are where I(t) is the source signature of the electric current injected into the circuit periphery and I(t) = 0 for t < 0. For a special case of an instantaneously-reacting dielectric slab (see equation ( 42)) we have Γ(s) = s, which simply yields Combining equations ( 28) with ( 29)-(34), the spacetime voltage distribution along the discretized circuit periphery and time axis is readily evaluated.In contrast to the real-FD solution that requires the matrix inverse at a large number of frequency points in combination with the inverse FFT, the marching on-in-time scheme (28) furnishes the stable TD voltage response in a given time window T via a single matrix inversion of the (time-independent) matrix For more details on the typical relaxation models and their behavior we refer the reader to the Appendix.Finally note that a dedicated numerical technique for handling L −1 has been proposed and successfully validated in [10].

Pulsed EM Emissions and Susceptibility
With the space-time voltage distribution at our disposal, the pulsed EM emissions radiated from a thin parallel-plane structure follow from [12, Eq. ( 9)] where E T;∞ is the amplitude radiation characteristic of the electric-type [11,Sec. 26.12], ξ is the unit vector of observation and τ = i 3 × ν is the unit vector tangential to ∂Ω.
Remarkably, the fundamental property of selfreciprocity allows to go yet another step further and carry out the radiated susceptibility analysis of the structure.Namely, if the emitting circuit is activated by the electric-current pulse I (t) applied at x S we may write [12,Eq. (22)] where V (x S , t) is the pulsed voltage response of the structure to a pulsed incident plane wave defined by its pulse shape e i (t) and by its direction of propagation and polarization specified by unit vectors β and α, respectively.
In a related work [14], the combination of equations (37)-(39) with the TD ray-type field expansion (see [15]) has proved to be computationally very efficient.Moreover, it has been recently shown that the voltage response given by (38) can be interpreted as the Thévenin voltage in the equivalent Kirchhoff-network representation [17].

A Numerical Example
In this section we shall analyze the irregularly-shaped parallel-plane structure shown in Fig. 4. The thickness of the structure is d = 1.50 (mm).The EM properties of the dielectric slab are described by its dielectric relaxation function (43) with = 4.50 0 and σ = 0.02 (S/m).
The structure is at x S = (25, 37.5) (mm) excited by the electric-current pulse described as where we take I m = 1.0 (A) and ct w = 0.10 (m) (see Fig. 5a).The radius of the excitation port is 1.0 (mm).The resulting pulsed voltage response has been evaluated at two field In order to validate the computational scheme described in Sec. 5, the TD voltage responses have also been calculated using the finite integration technique (FIT) as implemented in CST Microwave Studio .The resulting pulse shapes are shown in Fig. 6.As can be seen, the pulse shapes agree well.The observable discrepancies can be largely attributed to different boundary conditions along ∂Ω.While the TD-CIM model prescribes the perfect magnetic wall along the circuit's rim, the referential (three-dimensional) FIT model has the 'open' boundary condition on the faces of its square bounding box.Thanks to the spatial dimensionality reduction from 3 to 1 that is implicit in any CIM formulation, the reduction of computational costs with respect to standard direct-discretization computational techniques is typically (at least) 3 orders of magnitude.
In the final example, the evaluated space-time voltage distribution along ∂Ω × T is used in ( 37) to calculate the TD radiated amplitude at a chosen direction specified by ξ = (0, √ 2/2, √ 2/2).For validation purposes, the FIT model has been now irradiated by a uniform EM plane wave whose pulse shape e i (t) is related to the excitation electric current according to (39) (see Fig. 5) and whose direction of propagation and polarization are described by unit vectors β = −ξ = (0, − √ 2/2, − √ 2/2) and α = (0, − √ 2/2, √ 2/2), respectively.Figure 7 shows the FIT-based plane-wave voltage response V (x S , t) together with the TD-CIM-based, co-polarized radiation amplitude in the backward direction α • E T;∞ (−β, t).A good agreement between the calculated signals thus validates the relations discussed in Sec.5.1 as well as the numerical solution given in Sec. 5 once more.

Conclusions
A unified reciprocity-based contour-integral approach for analyzing transfer and EM emission/reception characteristics of a dispersive parallel-plane structure has been presented.It has been shown that the complex-FD reciprocitybased relation of the time-convolution type is a convenient starting point for proposing new CIM computational schemes in both real-FD and TD.Examples of such numerical schemes have been given and validated using a real-FD closed-form solution and TD numerical full-wave simulations.Subsequently, a few applications of the TD-CIM have been discussed.In particular, it has been shown that the TD radiated emissions can be readily calculated either from the voltage space-time distribution along the circuit's periphery or from the pulsed plane-wave response in the relevant receiving situation.An evidence of the property of self-reciprocity has been furnished via a numerical experiment employing both the described TD-CIM and the (three-dimensional) FIT.
excitation coefficients and U and H are [N × N] 2D-arrays.Clearly, the corresponding [N × N] impedance matrix directly follows from (10) as Z = U −1 • H, where U −1 is the matrix inverse for U.

Fig. 3 .
Fig. 3. Input impedance of the rectangular circuit evaluated using FD-CIM and the eigenfunction expansion method (EIG-E).(a) Point-matching solution with EIG-E.(b) Pulse-matching solution with EIG-E.(c) Absolute error of the FD-CIM solutions with respect to the referential EIG-E formula.