Planar Multilayer Structure Analysis : an Educational Approach

This paper discusses how symbolic computation combined with a circuit model can be used for analyzing planar multilayer structures, in a manner suitable for educational approach. Working in the Fourier domain, expressions for the transversal spectral Green’s functions are evaluated in compact, closed form using the symbolic computation capability of the Mathematica package. Printed antennas were analyzed through the method of moments. Further validation was achieved with the IE3D and HFSS packages.


I. INTRODUCTION
Sources positioned on the interfaces of planar multilayer structures can be used as a practical model in a wide variety of applications such as geophysical investigations, remote sensing, optoelectronics, microwaves, military surveillance and antenna theory [1], [2].For accurate results, the analysis of these structures requires numerical methods such as the Method of Moments (MoM), Finite Elements or Finite Differences.It is presently widely accepted that MoM-based algorithms are suitable for rigorous numerical analysis of printed structures of small to medium sizes (in terms of wavelength) stacked up in layers [1].However, for the application of this method, either in the spectral or in the spatial domain, the corresponding Green's functions need be derived.Due to their complex calculations, an essential question is posed: which is the most proper approach for determining the entire set of Green's functions, particularly for educational purposes?
In the spatial domain, Green's functions for multilayer media are traditionally represented by Sommerfeld's integrals.Due to the oscillatory nature of these integrals, their numerical computation is inefficient and time-consuming [1].Consequently, this approach is not recommended for a first graduate course on applied electromagnetism.On the other hand, spectral Green's functions can be derived in closed form if appropriate techniques are utilized [3]- [5].However, the calculations of spectral fields in multilayer structures are usually tedious and error-prone if done by hand, what is especially true in the analysis of structures with anisotropic and bi-isotropic materials [6], [7].
A planar multilayer structure composed of N + 2 isotropic, linear and homogenous layers stacked up in the z direction is shown in Fig. 1.The layers are assumed to be unbounded along the x and y directions.The lower layer, having complex permittivity ε g and complex permeability µ g , is denoted the ground layer and occupies the negative-z region.The next N layers are characterized by thickness ℓ n , complex permittivity ε n and complex permeability µ n , where 1 ≤ n ≤ N. The planar interface z = d N separates the N-th layer from free space (the upper layer).Metallic patches are printed at arbitrary positions on each one of the N + 1 interfaces of the structure.The surface current densities on these elements are given by , where boldface letters represent vectors, and x and y are the unit vectors along the x-and the y-directions, respectively.These current densities are considered the virtual sources of the electromagnetic field within each layer.Since the layers are free of sources (the sources are located on the interfaces of the multilayer structure, as depicted in Fig. 1), the electric field in the spatial domain of a monochromatic wave can be written, according to spectral techniques, as a superposition of plane waves, as follows [10] where the function E (k x , k y , z) is the spectral electric field, and k x and k y are the spectral variables.
For the n-th confined layer, located between the interfaces z = d n-1 and z = d n , the spectral fields E n (k x , k y , z) and H n (k x , k y , z) can be written as the superposition of two plane waves traveling in the ± z directions, that is where e nτ (k x , k y ) and h nτ (k x , k y ) are the amplitudes of the spectral fields, . Since the propagation function is given by z k i nz e τ for the n-th confined layer then, for a plane wave traveling in the +z direction, the square root in the propagation constant k nz1 needs to be adequately evaluated; in this case, Re{k nz1 } ≤ 0 and Im{k nz1 } ≥ 0. Similarly, for the plane wave traveling in the -z direction, the square root in k nz2 is calculated considering that Re{k nz2 } ≥ 0 and Im{k nz2 } ≤ 0. With these restrictions on the signs of the real and imaginary parts of k nzτ , Sommerfeld's radiation condition [14] is satisfied.From Maxwell's equations and using the symbolic computation capability of the Mathematica package, the amplitudes of the spectral fields in the x and y directions (e n xτ , e n yτ , h n xτ , and h n yτ , with τ = 1 or 2) can be expressed in terms of the field amplitudes in the z direction, e n zτ and h n zτ , according to [5].They are listed below and the Mathematica command lines utilized for their determination are given in the Appendix.
In order to establish the full-wave equivalent circuit for the n-th layer, firstly the x-and y-component of the spectral fields, given by equations ( 2)-(3), are specified at the z = d n interface, producing the following set of four equations in the unknown amplitudes e n zτ and h n zτ where with . Again, symbolic computation was performed in the Mathematica package.As an example, the calculation of the first row of the matrix n Z ~ is given in the Appendix.
Consequently, for the n-th confined layer, a two-port full-wave equivalent circuit, as shown in Fig. 2, can be established.
An equivalent procedure is employed to evaluate the fields at the free space and the ground layer interfaces.The boundary conditions are then applied for the planar interfaces, resulting in the circuit elements shown in Table I.Then, the multilayer structure depicted in Fig. 1 can be modeled by the circuit illustrated in Fig. 3. where, ) ( . Since the propa gation functions for free space and the ground were chosen to be , respectively, Sommerfeld's radiation condition is satisfied if the square roots in the propagation constants k 0z and k gz are calculated as follows: Re{k 0z } ≥ 0, Im{ k 0z } ≤ 0, Re{k gz } ≥ 0 and Im{k gz } ≤ 0.
Before finishing this section, it is important to point out that the theory developed here can also be applied to planar multilayer structures where the metallic patches are positioned inside the confined layers.
For these cases, the metallic patch can simply be considered to be at the interface between two adjacent layers having the same electrical characteristics.Hence, the developed theory can be utilized without any change.

III. APPLICATIONS
As a first application, a printed Yagi-Uda antenna fed by a delta-gap generator was investigated.In this case, the multilayer structure, shown in Fig. 4, is composed of the free space-substrate-free space layers.The radiator consists of two narrow flat dipoles printed on the dielectric substrate of thickness d 1 , electric permittivity ε 1 and magnetic permeability µ 0 .The length and width of the dipoles are denoted where η 0 is the intrinsic impedance of free space and ε r1 is the relative electric permittivity of the substrate.
Notice that the transversal Green's functions in (19) to ( 22) are equivalent to the ones presented in [16], where they were evaluated analytically, i.e., without employing a full-wave equivalent circuit to represent the planar multilayer structure.Once the Green's functions are derived, the next step is to determine the Yagi antenna's electric characteristics, such as its input impedance and radiation pattern.This involves solving integral equations constrained by the pertinent boundary conditions.In the present case, the following equations are obtained by enforcing the tangential electric field component being zero on the perfect conducting surface of the antenna elements, Using the MoM procedure -the most widely utilized numerical technique for rigorous analysis of printed geometries on multilayer planar media -to solve the system of equations ( 25)-( 26), the electric current densities J xv (x, y) are expanded in a set of basis-functions as follows, , on the driven element (27) , on the parasitic element (28) where J x1m and J x2n are the basis-functions, and I x1m and I x2n are coefficients to be determined.
After applying Galerkin method [11] and Parseval's theorem [12], the following linear system is established where J xv (k x , k y ) is the Fourier transform of J xv (x, y), S 1 is the driven element surface, p = 1, 2, … , M, and Piecewise-linear sub-domain rooftop basis-functions (taking into account the edge condition) were used for expanding the surface current densities J xv (x, y), so that As a delta-gap generator is used to feed the driven element, then is the Dirac's delta function.Consequently, the excitation equation ( 30) can be rewritten as with M being an odd number.
Using this approach, a Mathematica-based CAD software capable of performing the analysis of printed antennas was implemented.For greater efficiency in the numerical calculation of the elements Z pm , …, Z qn of the impedance matrix [Z] (equation ( 29)), mathematical procedures -like parity analysis of the Green's functions, change of the coordinate system (rectangular to polar) and asymptotic extraction technique [13] -were applied.To illustrate their effect, after changing the coordinate system -(ρ , φ) are the polar coordinates -and undertaking the parity analysis, equation (30) assumes the simpler form [14] For applying the asymptotic extraction technique, first the asymptotic term is calculated, resulting Consequently, Z pm can be evaluated from The first double integral in (42) depends on the operating frequency, whereas the second one does not.
These are the well-known Sommerfeld's integrals, whose integrands exhibit singularities in the form of branch points and poles, so that their numerical calculation requires careful attention.The poles, often complex, correspond to surface and leaky waves that can be excited in the dielectric layer.The number of poles depends on the thickness of the layer, its electrical permittivity and the wave number.An efficient way to calculate these integrals is the use of a deformed path of integration, with adaptive routines in the vicinity of the poles and branch points.Details on the numerical integration can be found in [15] - [17].In this work, a triangular path was chosen to calculate the elements Z pm , …, Z qn of the impedance matrix.
Finally, as the antenna is fed directly by a delta-gap generator, it is important to point out that the dielectric thickness poses no limitation on the application of the present theory.
Since M is an odd number, the antenna input impedance can be promptly calculated using the relation where E 0 was selected as 1 V.
The radiation pattern was calculated using the following expressions, derived from the stationary phase method [14]:   The agreement between our results for the input impedance and radiation pattern with those simulated in HFSS and IE3D certifies the accuracy of the technique presented in this work.
As a second application, the three-layer structure (perfect ground plane -substrate -free space) shown in Fig. 8 was analyzed.The radiator in question consists of a flat dipole printed atop the dielectric substrate of thickness d 1 , electric permittivity ε 1 and magnetic permeability µ 0 , such that 2w 1 denotes the dipole width and L 1 its length (L 1 >> 2w 1 ).The transversal Green's functions at the dipole interface z = d 1 are easily computed from the circuit model depicted in Fig. 9, as follows [14] 1 1 Fig. 9. Full-wave equivalent circuit for the structure depicted in Fig. 8.
Graphics for the input impedance and radiation pattern of a flat dipole (L 1 = 53.134mm and 2w 1 = 3.0 mm,

Fig. 4 .
Fig.4.Geometry of a Yagi-Uda antenna excited by a delta-gap generator.
denotes the tangential electric field component created by the antenna feeder.Since L v >> 2w v , (v ∈ {1, 2}), the unknown consists of just the x-components of the surface current densities J xv (x, y) excited on the dipoles.

Fig. 7 .
Fig.7presents the results for the yz-plane radiation pattern calculated at 2.0 GHz.

Fig. 8 .
Fig. 8. Geometry of a flat printed dipole excited by a delta-gap generator.

Fig. 10 .Fig. 11 .
Fig. 10.Dipole input impedance.Again, the results for the dipole input impedance and radiation pattern compared with simulations performed in HFSS and IE3D validate the accuracy of our CAD.