Nyström-Type Technique for Electromagnetic Wave Scattering in Inhomogeneous Material, Plasma and Metamaterial Slabs

We present a simple, stable, and spectrally-accurate quasi-analytical method for studying the reflection, transmission, and radiation of EM waves in the presence of single-layer or multilayer material, plasma, or metamaterial slabs that are inhomogeneous along the normal to the slab interfaces. Our approach formulates the problem as a linear Volterra integral equation of the second kind, discretized using an entire-domain Nyström method with a high-order Gauss-type quadrature. The proposed algorithms ensure mathematically guaranteed spectral convergence, handle electrically thin to electrically thick slabs, and can accommodate discontinuous inhomogeneity profiles. Application examples include one-dimensional band-gap structures, field enhancement in inhomogeneous metamaterial slabs with zero-crossings in permittivity, and dipole-driven plasma antennas.


I. INTRODUCTION
I NHOMOGENEOUS planar layers have diverse appli- cations in integrated optics, photonics, plasma, and microwave engineering [1], [2], [3], [4], [5], [6], [7].Gradedindex slabs, for instance, are used for optical waveguide design and the creation of optical components.Slabs with periodically modulated refractive indices exhibit frequencyselective behaviour, serving as simple one-dimensional EM band-gap structures.Inhomogeneous plasma layers are employed in radio communications, radio astronomy, and stealth applications due to their adaptable reflection, absorption, and transmission characteristics.Other significant applications include geophysical exploration, biomedical imaging, and optical tomography.
A distinct quasi-analytical method, which combines differential and integral equation techniques and exhibits spectral convergence, has recently been used to investigate EM scattering from inhomogeneous isotropic or anisotropic cylindrical [26], [27] or spherical [28] objects.In this study, we extend the application of this method to address scattering by single and multilayered slabs made of materials, plasma, or metamaterials exhibiting inhomogeneity along the normal to the slab interfaces.We consider two types of incident fields: plane waves and spherical waves radiated by elementary dipoles.Inhomogeneity profiles can be smooth or piecewise smooth spatial functions.Our primary goal is a simple, stable, and spectrally accurate algorithm that can provide machine-precision results.We will approach this in two steps.
We start by reformulating Maxwell's equations in the region occupied by the slab as a first-order system of ordinary differential equations (ODEs) with a spatiallyvarying coefficient matrix (Sections II-IV).When addressing the associated initial value problem (IVP), a variety of efficient and reliable numerical schemes are available, typically divided into linear multistep and Runge-Kutta methods.However, selecting the optimal scheme can be a complex task due to the trade-offs between stability, convergence properties, computational complexity, and memory usage.Additionally, these methods are limited by their algebraic order of accuracy and tend to degrade in performance when faced with discontinuous inhomogeneity profiles.
Recognizing these limitations, we depart from the conventional numerical schemes commonly associated with IVPs in the second step.Instead, we transform the system of ODEs into an equivalent linear Volterra integral equation (VIE) of the second kind.The transformation is followed (Section V) by discretization of the resulting VIE using a specialized high-order entire-domain Nyström method, coupled with an appropriate Gauss-type quadrature.
The proposed algorithms ensure mathematically guaranteed spectral convergence, constrained only by the smoothness of the underlying inhomogeneity profiles.When dealing with analytic spatial functions, the convergence is exponential.This stands in contrast to widely used generalpurpose methods that exhibit an algebraic convergence rate of O(h p ); in such methods, the convergence rate is determined by the order p regardless of the regularity of the underlying function.Additionally, our algorithms provide simple closedform expressions for every entry of the Nyström matrix, apply across a wide range of slab thicknesses, and can efficiently handle profile discontinuities.As a result, we attain machine-precision results at low computational and implementation costs.
Throughout this work, we assume a time-harmonic dependence of exp(iωt), which we omit for brevity.

II. SCATTERING OF E-POLARIZED PLANE WAVES
Fig. 1 shows an inhomogeneous slab, region 1, with thickness d and continuous constitutive parameters ε 1 (z) = ε 0 ε r1 (z) and μ 1 (z) = μ 0 μ r1 (z).Sectionally continuous parameters will be addressed in Section VI.The slab extends infinitely in the x and y directions, surrounded by homogeneous regions 0 (z > z 1 = 0) and 2 (z < z 0 = −d) with constitutive parameters (ε 0 , μ 0 ) and (ε 2 , μ 2 ), respectively.Both (ε 1 , μ 1 ) in region 1 and (ε 2 , μ 2 ) in region 2 can be complex to account for potential losses.The primary excitation is an incident E-polarized plane EM wave propagating in the direction of kinc = −k x x − k z ẑ, where In (1), k 0 = ω √ ε 0 μ 0 , and ψ is the incidence angle measured from the z-axis, increasing clockwise.The electric field, with an amplitude of E 0 , is given by Remark 1: In each region, the field components vary as f (z)e ik x x , with the phase factor exp(ik x x) dictated by the incident wave.
In regions 0 and 2, the nonzero component of the total electric field is expressed as where R and T are the refl9ection and transmission coefficients, respectively.In (2b), For a lossy medium in region 2, k 2 2 − k 2 x falls within the lower half-plane of the complex plane.To ensure a physically valid solution, we choose the branch of the square root defining k z2 in such a way that Re(k z2 ) ≥ 0 and Im(k z2 ) ≤ 0. The components of the associated magnetic field are obtained from For the field in region 1, with components F(x, z) = e ik x x F(z) (F = E y1 , H x1 , H z1 ), from Maxwell's equations we obtain an algebraic equation: and an ODE: where the column-vector variable u(z) is defined as and the matrix A 1 (z) is given by

A. SOLUTION OUTLINE
The general solution of ( 5) is: where W 1 (z, z 0 ) is the transition matrix associated with the system matrix A 1 (z).The transition matrix is the solution of the matrix-valued IVP: where I 2 denotes the 2 × 2 identity matrix.
The highly accurate computation of W 1 (z, z 0 ), carried out in Section V, is an integral part of our method.After obtaining W 1 (z 1 , z 0 ), we enforce continuity on ū(z) at z = z 0 and z = z 1 , leading to: where and is the 2 × 2 block matrix: Once we obtain the unknowns R, T, and ū(z 0 ) from ( 10), we can evaluate the field everywhere using (2)-( 4) and (8).

III. SCATTERING OF H-POLARIZED PLANE WAVES
In this section, we investigate the dual problem to the previous one.The incident plane wave is now H-polarized, and its magnetic field is given by In regions 0 and 2, the nonzero component of the total magnetic field is expressed as where R and T are the reflection and transmission coefficients, respectively.The components of the associated electric field can be obtained from with ε = ε 0 or ε 2 .
For the field in region 1, with components F(x, z) = e ik x x F(z) (F = H y1 , E x1 , E z1 ), from Maxwell equations we obtain an algebraic equation: and an ODE: where the column-vector variable υ(z) is defined as The matrix A 2 (z) is the dual of A 1 (z) and is given by The general solution of ( 18) is: where W 2 (z, z 0 ) is the solution of the matrix-valued IVP: The transition matrix W 2 can be computed as in Section V. Enforcing continuity on the vector variable ῡ(z) at z = z 0 and z = z 1 , yields: with h defined in (11).Equation ( 23) is the dual of (10).

IV. ELECTRIC AND MAGNETIC DIPOLE RADIATION
In Sections IV-A and IV-B, we deal with arbitrarily oriented electric and magnetic dipoles p = p x x + p y ŷ + p z ẑ, m = m x x + m y ŷ + m z ẑ, located at r (x , y , z ) over a bilateral [Fig.2(a)] or grounded [Fig.2(b)] slab.Section IV-C focuses on dipoles embedded in either of these slabs, whereas Section IV-D concerns dipoles radiating in the presence of DNG slabs.For simplicity, we assume that region 2 (z < z 0 ) in the structure of Fig. 2(a) has the properties of vacuum (ε 0 , μ 0 ).However, the analysis readily extends to include the general case where a homogeneous medium (ε 2 , μ 2 ) occupies that region.

A. DIPOLES ON TOP OF AN INHOMOGENEOUS SLAB
We start with the problem of electric and magnetic dipoles radiating on top of the inhomogeneous slab of Fig. 2(a).

1) FIELD REPRESENTATION INSIDE THE SLAB
In the 2D Fourier-transform domain defined by the pair:

+∞
−∞ e i(k x x+k y y) f (x, y, z) dx dy,(24b) Maxwell's equations for region 1 split into two algebraic equations: together with four first-order scalar differential equations having the components of the column-vector variable ū1 (z) = Ẽx1 (z), Ẽy1 (z), Hx1 (z), Hy1 (z) (25) as the unknowns.Here, for simplicity, we use fq1 (z) as shorthand for fq1 (k x , k y ; z), where f can be either E or H, and q can be either x or y.Let After suitable manipulation, we obtain the following ODE: where with The general solution of ( 27) is: where W 3 (z, z 0 ), the transition matrix of the linear differential system (27), satisfies the IVP: with I 4 denoting the 4 × 4 identity matrix.The computation of W 3 (z, z 0 ) is detailed in Section V.The term w1 (z 0 ) encountered in (30) is elaborated in Section IV-A4.Once w1 (z) is determined, we can evaluate ū1 (z) using (26): 2) FIELD REPRESENTATION IN REGION 2 In the homogeneous source-free region 2, the field can be described using scalar electric (π e2 ) and magnetic (π m2 ) Hertz potentials [29]: Both π e2 and π m2 satisfy homogeneous Helmholtz equations: with Fourier-transform counterparts: Using, for simplicity of notation, the shorthand πq2 (z) for πq2 (k x , k y ; z), we obtain: Here, a 2 = a 2 (k x , k y ) and b 2 = b 2 (k x , k y ) are coefficients to be determined, and Selecting the branch of the double-valued function γ as in (38), ensures the validity of the radiation condition at z → −∞.

3) FIELD REPRESENTATION IN REGION 0
The field in region 0 can be written as a sum of two components: The first component, denoted by the subscript d, represents the primary field that would be generated by the dipole source radiating in an unbounded space with parameters (ε 0 , μ 0 ).We introduce the vector variable wd (z), defined similarly to w 1 (z) but with the index 1 replaced by d in both (25) and (26).It can be shown [29] that where and cd is defined as follows: • For the electric dipole case: Here, s = sign(z − z ).
• For the magnetic dipole case: The second component, denoted by the subscript s, represents the secondary field, also known as the scattered field.It satisfies the source-free Maxwell equations and can be expressed in terms of two scalar Hertz potentials, π s e0 and π s m0 , using formulas similar to (33) and (34).The Fourier transforms π s e0 and π s m0 are given by where a 0 = a 0 (k x , k y ) and b 0 = b 0 (k x , k y ) are unknown coefficients to be determined.Following a similar analysis as in Section IV-A2, we arrive at where ws (z) is defined in a similar manner to w 1 (z), with the index 1 replaced by s in both ( 25) and ( 26).

4) FINALIZING THE SOLUTION
To determine the yet-unknown coefficients {a 0 , b 0 , a 2 , b 2 } and subsequently find the Fourier transforms of the fields, we impose continuity on both the electric and magnetic fields' x and y components-equivalent to enforcing continuity on w(z)-at z = z 0 and z = z 1 .Using (30), we obtain a 4 × 4 system of linear algebraic equations: Once these coefficients are determined, we can evaluate ws (z) and w2 (z) using ( 47) and (39), respectively.Subsequently, we can obtain w1 (z) by using (30) with the condition w1 (z 0 ) = w2 (z 0 ).The field can then be computed everywhere by applying the inverse Fourier transform (24a).

5) FAR FIELD
In the radiation zone (r λ 0 ), the generated field consists of two parts: a space wave component, representing the 3D radiation field in regions 0 and 2, and a surface wave component, characterizing the 2D waveguide modal field across all three regions.The evaluation of these field components can be carried out analytically via contour integration techniques detailed in [30] and [31] (also refer to [32], [33]).Further details regarding their application to our specific problem can be found in Appendix A.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. DIPOLES ON TOP OF A GROUNDED SLAB
Refer to the configuration of Fig. 2(b).Now, region 2 is perfectly electrically conducting, and as a result, both E x and E y vanish at z = z 0 .Setting the continuity of w(z) at z = z 1 leads to the following equation analogous to (48): The expression of the far field in region 0 remains unchanged.

C. HORIZONTAL DIPOLES EMBEDDED IN A SLAB
Here, we consider electric and magnetic dipoles Region 0 is now source-free, and its field ( Ē0 , H0 ) can be expressed in terms of two scalar potentials, π e0 and π m0 , with Fourier transforms: Here a 0 and b 0 are coefficients to be determined.The corresponding vector variable w0 (z)-defined as in (26) but with the index 1 replaced by 0 in both ( 25) and ( 26)-is given by For the vector variable w1 (z) associated with the field in region 1 (z 0 ≤ z ≤ z 1 ), we can write The field representation in region 2 remains unchanged.By applying the boundary conditions at z = z 0 , z = z 1 , and z = z , we obtain the linear algebraic system: where For dipoles embedded in the grounded slab of Fig. 2(b), we obtain the following linear algebraic system analogous to (60): The expression of the far field in region 2 (case of bilateral slab) remains unchanged.For region 0, we use (49) with F θ (θ, ϕ) and F ϕ (θ, ϕ) equated to F s θ (θ, ϕ) and F s ϕ (θ, ϕ) as defined in (51).

D. DNG SLABS
Our analysis extends to arbitrary isotropic inhomogeneous layers.Initially, we focus on the structure shown in Fig. 2(a) (dipoles over a bilateral slab).We define hDPS = [a 0 , b 0 , a 2 , b 2 ] as the column vector containing the unknown coefficients encountered in (48) for a DPS slab with ε 1 (z) > 0 and μ 1 (z) > 0. We designate the associated transition matrix as W DPS 3 (z 1 , z 0 ).In the case of the DNG counterpart with −ε 1 (z) and −μ 1 (z), we use hDNG and WDNG

3
. By substituting (42) into (48), we obtain: where, for q ≡ DPS and for q ≡ DNG, q is the 4 × 4 block matrix The property In the far scattered field, where Consequently, ( DNG ) * = DPS , and (64) can be restated as: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
From ( 63) and (66), using (43), we deduce: As a result, the relationship between the radiation patterns FDNG (θ, ϕ) and FDPS (θ, ϕ) of the far scattered field, whose components are given by ( 51) and (52), is: For the structure in Fig. 2(b) (dipoles over a grounded slab), we substitute C(γ ) in (65) with and use the definitions in (54) for a 2 and b 2 .The final resultsequations ( 67) and (68)-remain unchanged.When a dipole is placed within any of the slabs in Fig. 2, all the previously discussed results and equations remain unchanged.However, in (67) and (68), the phase term −k 0 (z − z 1 ) cos θ should be omitted.
Under the specified conditions, transitioning from a DPS slab to its corresponding DNG slab does not alter the magnitude of the radiation pattern of the far scattered field.

V. NUMERICAL SCHEME FOR COMPUTING THE TRANSITION MATRICES
In this section, we present an overview of the numerical scheme used to compute the transition matrices W 1 , W 2 , and W 3 needed in equations ( 8), ( 12), ( 21), ( 23), ( 30), (48), ( 55), (60), and (62).Instead of focusing on a specific case, we provide a unified approach for solving the IVP where I represents the identity matrix of the same order as A.
In (69a), A can be A 1 , A 2 , or A 3 , and correspondingly, W represents W 1 , W 2 , or W 3 .
To solve the IVP in (69), we transform it into a secondkind Volterra integral equation (VIE): by integrating (69a) from z 0 to z and using (69b).To simplify notation, we introduce the change of variable: We then define: With these definitions, we express (70) as: To discretize the VIE (73), we employ the Legendre-Gauss-Radau (LGR) quadrature with nodes at the zeros is the k-th degree Legendre polynomial.The corresponding weights are: Using the expansion: with coefficients [27]: we derive the discrete counterpart of (73) as: Here [27], Subsequently, applying (76) at t = t 1 , t 2 , . . ., t N yields a system of linear algebraic equations: where Once we find {F(t i )} N i=1 from (78), we can evaluate F(t) at any t in (−1, 1) either using (76) or by interpolation [27]: The property W(z 1 , z 0 ) = F(1) = F(t N ) allows for the direct evaluation of W(z 1 , z 0 ) from (78) without the need for interpolation, justifying our preference for the LGR quadrature method in this paper.This property is also shared by other alternatives, including Gauss-Legendre-Lobatto and Clenshaw-Curtis rules.

A. ERROR ESTIMATE
Let F(t) be the exact solution of (73) and F N (t) its approximation obtained from either (76) or (80).A convergence analysis (see [35], [36]) provides the following error estimate: If F ∈ C m (−1, 1), then, for a sufficiently large N, F− c is a constant independent of N), indicating exponential convergence.These properties demonstrate that the proposed technique exhibits the typical convergence behaviour of a spectral method.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The algorithm's spectral convergence rate depends on, and benefits from, the smoothness of the underlying functions, distinguishing it from general-purpose methods that have fixed algebraic convergence rates (O(h p )); in such methods, p remains constant regardless of the regularity of the function.
Discontinuous (non-smooth) profiles can be addressed in the way outlined in Section VI.

B. SINGULAR PROFILES
For later reference in Section VII-D, we will discuss simplezero and double-zero inhomogeneity profiles, resulting in corresponding simple-pole and double-pole system matrices Here, z p represents the pole, and m is its multiplicity.The function With z = z(t) as in (71), we define: For ˜ (τ )F(τ ), we use an expansion of the form (74). Following a procedure similar to the one that led from (69) to (78), we derive a system of linear algebraic equations analogous to (78): for i = 1, . . ., N, where with . When t p is in the interval (−1, t), the integral in (84) should be interpreted as the Cauchy principal value for m = 1 and in the Hadamard finite-part sense for m = 2.
Once we find {F(t i )} N i=1 from (82), we can evaluate F(t) at any point, for instance, via interpolation based on (80).We can evaluate C (m) k (t) using the recurrence relation:

1) EVALUATION OF C
This leads to the following expression: The term C (0) k (t) that arises in (85) for m = 1 is given by: with C k (t) defined in (77).The initial values are: 1 t; t p = 1 + t + t p C (1) (2) 0 t; t p ∂t p ,

A. HANDLING DISCONTINUOUS PROFILES
The latter observation offers a simple and efficient method for handling single-layer slabs with discontinuous inhomogeneous profiles: Modelling the slab as a stack of layers by introducing artificial interfaces at the locations of the discontinuities and then applying (89).The algorithm achieves spectral convergence, with the rate depending on the smoothness of the parameters within each sublayer.

B. THICK SLABS
While the model in Fig. 3 primarily focuses on multilayer slabs, it also offers advantages when applied to single, electrically thick slabs.In the case of a single slab, the required quadrature order N increases with the electric size of the slab.As k 0 d becomes larger, the size of the Nyström matrix (N 2 ) and, consequently, the computational demands can become quite high and potentially prohibitive.To address this challenge, a thick slab can be artificially divided into a stack of L moderate-thickness layers, and the equation (89) can be used.This approach optimally maintains N at a manageable level, significantly enhancing the stability and computational efficiency of the algorithm.For more details, refer to Section VII.

VII. NUMERICAL RESULTS AND COMPARISONS
In this section, we present numerical examples to validate the algorithms and demonstrate the high accuracy and flexibility of the proposed method.For scattering problems involving E (H) polarized plane-waves, we assume that (ε 2 , μ 2 ) = (ε 0 , μ 0 ), and use an incident field amplitude of E 0 = 1V/m (H 0 = 1A/m), unless stated otherwise.
The computer programs were implemented in Mathematica 10.4, which employs machine precision arithmetic as the default setting.It uses double-precision floating-point numbers with approximately 16 decimal digits of precision.However, Mathematica also supports variable precision arithmetic, allowing users to specify their desired precision for numerical calculations.Precision can be set either globally, using the global variable $PreRead, or on a per-expression basis using the SetPrecision function.In our work, we consistently used a precision setting of 17 digits for all examples, except for Table 1 and Fig. 8.

A. CONVERGENCE ANALYSIS
In Fig. 4, we demonstrate the convergence of the algorithm by plotting the base-10 logarithm of the relative errors of |R| and |T| against the quadrature order N.These results pertain to the E-case with k 0 d = π and ψ = π/4.We consider six distinct inhomogeneity profiles, labelled as Case I through Case VI, each of which has exact analytical solutions available (see Appendix B).The relative error for X = |R| and X = |T| is defined as = |(X − X exact )/X exact |.The figure illustrates exponential convergence, consistently achieving an error on the order of 10 −16 , comparable to the round-off error, with the appropriate selection of N.
For a deeper insight into the convergence behavior of the method, Table 1 presents the relative error of |R| for various quadrature orders N with parameter values matching those in Fig. 4.These computations use highprecision arithmetic.Notably, Table 1 reveals a phenomenon of super-convergence: doubling the value of N yields more than double the number of correct digits in the computed reflection coefficient.This observation, combined with the results from Fig. 4, suggests that the algorithm's accuracy is only limited by the employed machine precision.These findings emphasize the exceptional precision and reliability of the method.
In Fig. 4, the slab had a thickness falling within the resonance region.However, the results shown in Fig. 5, with k 0 d at 20π, 200π, 500π , emphasize that the algorithm can handle significantly thicker slabs.By employing the model shown in Fig. 3 and appropriately selecting the number of artificial layers L and the quadrature order N, we consistently achieve an error on the order of 10 −16 in all cases, underlying the reliability of the method for very thick slabs.
In all previous examples, the profiles were continuous spatial functions.To explore the algorithm's behavior in the presence of discontinuous parameters, we now consider a slab with μ r1 = 1 and a relative permittivity given by: In Fig. 7, we show the base-10 logarithm of the relative error for | F(θ, ϕ)|, representing the radiation pattern of an electric dipole in the configuration of Fig. 2(a).We used the model from Fig. 3  In Fig. 8, we compare the convergence of our algorithm with that of the widely used SA method (staircase approximation) for the E-case with parameters ψ = π/4, k 0 d = 20π , ε r1 = 3 + cos [k 0 (z + d)], and μ r1 = 1.In Fig. 8(a), the relative error is plotted against N 2 L for our algorithm and against L for the SA method.Notably, SA requires nearly L = 10 6 layers to achieve a relative error of 10 −13 , whereas our algorithm, with L = 20 and N = 70, attains an error of 10 −100 .Moving to Fig. 8(b), we compare the time required for both methods to reach a set accuracy level, further highlighting the superiority of our approach.

B. APPLICATION EXAMPLE: PHOTONIC CRYSTAL
As mentioned in the Introduction, inhomogeneous slabs with periodically modulated refractive indices can function as 1D EM band-gap structures.In Fig. 9, we explore this possibility with d = 1m, ε r1 = 2 + 15 sin 4 (30π z) and μ r1 = 2.In Fig. 9(a), we show the first three band gaps in the transmittance |T| 2 for ψ = 0, indicating equal transmittance for both E and H polarizations. Focusing on the lower portion of the first band gap, in Figs.9(b) and 9(c), we depict transmittance variations with the angle of incidence in the E and H case, respectively.Notably, within the range 22 < k 0 d < 32, transmittance remains consistently below 10 −8 for both polarizations, regardless of the angle of incidence.

C. SCATTERING ON A PLASMA SLAB
We now consider an inhomogeneous cold collisional plasma layer with relative permittivity Here, ν en is the frequency of electron-neutral collisions, and ω p is the angular plasma frequency: ω p = e 2 n e (z)/ε 0 m, where n e , e, and m stand for the plasma density, electron charge, and electron mass, respectively.For generality, we assume that both n e and ν en vary spatially with distributions: shown in Fig. 10.Here, n 0 is a constant, and the parameter α can be set to either α = 0 (in which case ν en (z) = ν 0 remains constant) or α = 1.

D. FIELD ENHANCEMENT IN EPSILON-CROSSING-ZERO FILMS
Field enhancement occurs in structures with near-zero parameters [37].To illustrate, consider a thin (k 0 d << 1) homogeneous slab with permittivity ε 0 ε r , positioned in air and illuminated by an H-polarized wave.The boundary conditions E z (z − 1 ) = E z (z + 1 )/ε r and E z (z + 0 ) = E z (z − 0 )/ε r signify an amplification of |E z | within the slab, inversely proportional to ε r .As ε r approaches zero, this phenomenon intensifies, leading to significant energy confinement.
To investigate the field enhancement effect in inhomogeneous slabs with permittivity profiles that exhibit zero crossings, Fig. 16 presents the distribution of |E z (z)| within a plasma slab surrounded by air in the H-case, with parameters   permittivity, defined as: corresponds to (90) with ν 0 = 0 in (91).In Fig. 16(a), we select the circular frequency as ω = ω p (−d/2) = 5.64146 × 10 10 rad/s.With this choice, ε r1 (z) takes the form ε r1 (z) = (d + 2z) 2 (d 2 − 4dz − 4z 2 )/d 4 (inset), exhibiting a double zero at z = −d/2.As a result, substantial field confinement occurs at the slab's mid-plane.The results in Fig. 16(a) were obtained using the algorithm outlined in Section V-B (case m = 2).
In Fig. 16(b), we have chosen ω = 5.5 × 10 10 rad/s.For this specific choice, ε r1 (z) has zeros at z/d = −0.4208244and z/d = −0.5791756(see inset), leading to twin field confinement at these two planes.The results in Fig. 16(b) were obtained using the algorithm in Section V-B (case m = 1).

F. SOME FEATURES OF A DIPOLE DRIVEN PLASMA ANTENNA
At the frequency of f = 8 GHz, in Fig. 20 This alignment is easily explained by examining the plots in Figs.14(b) and 15(b), which reveal that, at f = 8 GHz, the transmittance is extremely small for any angle of incidence.This indicates that the interface z = z 0 = −d behaves like a perfectly conducting plane.Based on the results in Fig. 20, we can conclude that, at the selected operating frequency, the bilateral plasma slab of Fig. 2(a) serves as a practical implementation of the grounded plasma slab of Fig. 2(b), eliminating the need for any conductive elements.

VIII. CONCLUSION
We have introduced a simple yet highly accurate numerical scheme for studying scattering of plane and spherical electromagnetic waves on inhomogeneous materials, plasma, and metamaterial slabs.While our primary focus has been on scattering and radiation, the algorithms can also address the related propagation problems, as well as complex frequency eigenvalue and eigenvalue lasing problems [40].Furthermore, these algorithms can be extended to slabs containing complex inhomogeneous media that can be anisotropic, chiral, or bi-anisotropic.Such an extension would enable investigations into intricate EM problems across a broad spectrum of innovative materials and configurations.

APPENDIX A FAR SCATTERED FIELD
Here, we provide an overview of the application of the method of contour integration to our specific problem.We focus on a vertical electric dipole (p = ẑp) located at r = ẑz above the slab in Fig. 2(a), and limit our discussion to the analytical evaluation of the Hertz potential π s e0 as defined in Section IV-A3.The approach for evaluating the remaining Hertz potentials and the field within the slab is similar.
A closer examination of (78) reveals that the transition matrix W 3 follows the structure noted in [11]: The elements w ij are analytic functions of λ 2 (or equivalently, of γ 2 ).Once we determine a 0 from (48) using (93), we can derive π s e0 from (46) as follows: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where, To evaluate the inverse Fourier transform π s e0 of π s e0 , we introduce a change of variables: ) and follow the procedure outlined in [31].When z 1 = 0, as in Fig. 2, we obtain: The original inverse transform path C 0 takes the form shown in Fig. 21, accounting for the presence of the two branch points ±k 0 of the multivalued function γ (λ) and the surface wave poles of the integrand.These poles correspond to the roots λ = ±β 1 , . . ., ±β n of D[γ (λ)] = 0, representing the dispersion equation for the supported TM y modes.(Note: The pole configuration illustrated in Fig. 21 is specific to the lossless case.When the slab exhibits loss, these poles shift away from the Re(λ) axis.)The branch cuts of γ extend to infinity as illustrated in the figure.Additionally, the Hankel function introduces a branch point at λ = 0 with the corresponding branch cut extending to infinity along the positive Im(λ) axis.

A.SPACE-WAVE COMPONENT
The first term in (97) which contains the integral over C 2 represents the space-wave component of the potential.A steepest descent evaluation of this integral, as outlined in [31] and valid for large distances from the dipole, results in the closed-form expression: The correctness of (98) has been checked via an independent derivation using the stationary phase integration technique.

B.SURFACE-WAVE COMPONENT
The second term in (97) represents the surface-wave component of the potential.Here, for large distances from the dipole, it is appropriate to use the large-argument asymptotic form: In obtaining the Res λ→β j , we employ the chain rule: The resulting expression for the surface-wave component is: . The term D (γ j ) can be easily obtained from (95) provided the derivative w ij (γ ) can be computed.

1)COMPUTATION OF {β j } n j=1
As previously mentioned, the elements w ij with 1 ≤ i, j ≤ 4 of the transition matrix W 3 are analytic functions of λ 2 (or γ 2 ).Consequently, D(γ ), as defined in (95), is an analytic (i.e., pole free) function of γ .This property enables us to  compute the roots {γ j } n j=1 of D(γ ) = 0 using well-established algorithms, such as the one described in [41].Once the set of roots {γ j } has been obtained, we compute β j as β j = γ 2 j + k 2 0 .As an illustrative example, Table 2 displays the normalized propagation constants β/k 0 for the TM y modes supported by the bilateral slab of Fig. 2

APPENDIX B EXACT ANALYTICAL SOLUTION FOR SPECIAL CASES
We will focus on the scattering problem in the E-case, noting that the H-case can be addressed through duality.Let Ē1 = ŷE y1 (x, z) = ŷe ik x x Ẽy1 (z).
Let V + (z) and V − (z) be two linearly independent solutions of (102).Then where c 1 and c 2 are constant coefficients.In terms of Ẽy1 , we can obtain Hx1 from By enforcing continuity of both E y (z) and H x (z) at the interfaces z = z 0 and z = z 1 , using (2), ( 3), (103), and (104), we obtain a 4 × 4 linear system of algebraic equations from which R, T, c 1 , and c 2 can be determined.At this stage, the solution procedure is purely formal because V ± (z) cannot be expressed in terms of known special functions for general inhomogeneity profiles.In Table 3, we provide a summary of special profiles for which exact  Case II: where D ν (z) is the parabolic cylinder function, and: g ± (z) = ±1 + i 2α 2 3/4 (α 1 + 2k 0 α 2 z), Case III: Case V: Here MathieuC[a, q, z]/MathieuS[a, q, z] denote the even/odd Mathieu function with characteristic value a and parameter q.The Mathieu functions are solutions to the equation y + [a − 2q cos(2z)]y = 0.
Case VI: where Y 1 denotes the Bessel function of the second kind, and with q(z) = az k x 2 − bk 0 2 + i bk 0 2 − k x 2 .For additional inhomogeneity profiles with exact solutions to the wave equation, see [29,Sec. 4.3.4].

FIGURE 2 .
FIGURE 2. (a) Bilateral inhomogeneous slab excited by an arbitrarily oriented elementary electric or magnetic dipole.(b) Grounded inhomogeneous slab.

s
(z) denotes the jth scalar component of ws (z).•In region 2:

TABLE 1 .
Convergence of relative error for |R| in the E-case with k0d = π and ψ = π/4.Constitutive parameters εr1(z) and μr1(z) in Case I-Case III are the same as in Fig. 4.

2 54
d 2 /(z + d) 2 , − d 2 < z < 0. This profile has a second-kind (jump) discontinuity at z = −d/2.Following the framework outlined in Section VI, we use the model of Fig. 3 with L = 2, introducing an artificial interface at z = ζ 1 = −d/2, and use (89).In Fig. 6, we present the relative errors of |R| and |T| as N increases for the E-case, with k 0 d = 5π and ψ = π/4, demonstrating exponential convergence similar to our previous examples.

FIGURE 6 .FIGURE 7 .
FIGURE 6. log 10 of the relative error for |R| and |T | versus N in the E-case with a discontinuous permittivity profile (inset) at k0d = 5π and ψ = π/4.

FIGURE 21 .
FIGURE 21.Path of integration in the complex λ− plane.
ae bz /b), where J λ (•) denotes the Bessel function of the first kind.Case IV: V ± (z) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Upon differentiating with respect to γ , we obtain: k (t) = I. j ; γ , i = 1, . . ., N, (101)