Electromagnetic Scattering by . Stratified Inhomogeneous Anisotropic Media

An analytical formulation is presented for the computation of scattering and transmission by general anisotropic stratified material. This method employs a first-order state-vector differential equation representation of Maxwell’s equations whose solution is given in terms of a 4 x 4 transition matrix relating the tangential field components at the input and output planes of the anisotropic region. The complete diffraction problem is solved by combining impedance boundary conditions at these interfaces with the transition matrix relationship. A numerical algorithm is described which solves the state-vector equation using finite differences. The validation of the resultant computer program is discussed along with example calculations.


T HE CONSIDERATION of electromagnetic interaction
with anisotropic materials has been a topic of considerable interest in recent years. For the most part, this interest has been manifested in a diverse realm of particular research areas such as ionospheric and magnetospheric propagation, crystal and semiconductor physics, opto-electronics and composite material scattering, as well as many others. In addition, the employment of anisotropic materials in electromagnetic engineering applications is becoming increasingly widespread.
Some notable examples are in microstrip and active device substrates, novel antenna radome structures and radar absorbing materials.
In many instances of practical importance the anisotropic material structure can be approximated as having an inhomogeneity in only one dimension. Such stratified media appear in the description of optical coatings, earth strata, ferrite absorbing coatings, ionosphere models, etc. Early published efforts by AbelCs [I] and Jones [2] in considering stratified media were resticted to the isotropic case. An excellent review of propagation through isotropic stratified media is given by Born and Wolf [3]. By extending some of the original concepts employed in the works of Abelb and Jones, the case of stratified anisotropic material was first considered by Teitler and Henvis.[4] and later by Berreman [ 5 ] . In these two efforts, a 4 x 4 matrix method is introduced and applied to the analytical treatment of propagation through a single layer of various special case anisotropic materials bounded on both Manuscript received March 21, 1986.  Our goal in this effort was to develop a general purpose formulation and efficient numerical algorithm for scattering and transmission calculations involving stratified anisotropic media, [ 101. By combining this algorithm with a numerical optimization procedure it will become possible, in a continuing effort, to synthesize multiple layered structures to achieve desirable polarization filtering characteristics for a postulated range of incident field wave vectors (specified by incident angle and temporal frequency). Of special interest is the case of scattering by coated metallic surfaces.
The methodology described herein is an extension of the technique developed by Ruck et al. [ 1 11 for use with isotropic material. Maxwell's equations are cast into a 4 X 4 matrix formulation which is equivalent, but not identical, to that in [ 5 ] . Solution procedures based upon both eigenfunction expansions and a new finite-difference approach are considered. A computer program is described and the various means of validating its performance are discussed along with accompanying sample computations.

GENERAL FORMULATION
The geometry being considered is shown in Fig. 1, wherein a one-dimensionally inhomogeneous planar anisotropic region is bounded by homogeneous isotropic material half-spaces for z < 0 and z > d . A monochromatic plane wave is obliquely incident from the left half-space (region 1) with arrival angle 0, in the xz plane of incidence. Using an exp ( j u t ) time convention, the complex phasor fields may be written in the separable product form where k, = kl sin is the x-component of the incident field wave vector R1. Maxwell's source-free equations incorporate 3 X 3 complex permittivity and permeability tensors, E ( Z ) and p(z), that contain the effects of dielectric and magnetic losses, as well as electrical conductivity, and using i, s, and t superscripts to denote respective incident, scattered and total fields, we define 2 x 2 complex scattering . and transmission matrices, S and T, by . . . .

E'",(O)
Note that the total field for z 2 d is also the transmitted field into region 3. By partitioning the 4 X 4 transition matrix into four 2 x 2 quadrant submatrice-, After substituting (1) into (2), followed by cancellation of the common exponential factor, the resulting six scalar component equations can be further reduced to four independent equations through the algebraic elimination of the zcomponents of and fi. By defining a state vector in terms of the transverse field components of i? and fi in (l), the resulting four coupled linear first-order ordinary differential equations can be expressed in matrix notation by The complex elements of the 4 X 4 r-matrix are given in the Appendix.
The overall diffraction problem, including specular scattering and refractive transmission, can be solved by first undertaking the computation of the complex 4 X 4 transition matrix A. This matrix provides a linear relationship between the state vectors at the left and right planar boundaries of the inhomogeneous media.
there results from (5) and (6) Ei We introduce wave impedance matrices for the homogeneous half-spaces in regions 1 and 3 which relate the transverse field vectors, where, for rn = 1 or 3, with characteristic impedance, and In providing this relation between the fields at the interfaces, the A-matrix incorporates all of the effects of the intervening material composition in region 2. The computation of the A-matrix is obtained via the numerical solution of (4), which may be implemented through any of a variety of algorithms. Two of the possible techniques will be considered in the next section for the case of layered media.
Using the hypothetical A-matrix we can solve the diffraction problem in a direct manner. By introducing tangential Snell's law, as given by ( I~c ) , is still valid for the case of anisotropic media in region 2.
The S and T matrices may now be found through simple matrix manipulations, after substituting (11) into (lo), where we have defined new impedance matrices, For the important case of a perfect conductor in region 3 , with Z3 = 0, the T-matrix becomes zero and the scattering matrix simplifies to We will now consider the major computational effort in solving the diffraction problem, which is in obtaining the transition matrix A.

TRANSITION MATRIX COMPUTATION
Consider the special case of the anisotropic media in region 2 being composed of N'homogeneous layers, as illustrated in Fig. 2, where z1 = 0 and Z.hr+ = d. The incident field halfspace (region 1) will be denoted as the zeroth layer while the left half-space (region 3) is enumerated as the N + 1 layer.
Our initial goal will be to describe procedures for obtaining local transition matrices relating the state vectors at adjacent "n" and ''n -1" interfaces, Since r ( z ) = rn is constant within the nth layer, the resultant state vector solutions of (4) will have the form of nonuniform plane waves +n(z)= v n exp ( 7 ,~) ( 1 8) where 7, is a constant vector. Acceptable complex z-directed wavenumbers 7, are found by first substituting (1 8) into (4) to produce the eigenvalue system with I being the 4 X 4 identity matrix. The eigenvalue wavenumbers may be found as the roots of the quartic characteristic equation obtained by setting the determinant of the bracketed matrix in (19) to zero. An alternate derivation of this dispersion relation is given by Damaskos and Ushlenghi [ 131. In addition to direct numerical root solving routines, such as the Newton-Raphson method, there exist several more efficient procedures for evaluating eigenvalues (e.g. QRfactorization) based upon matrix transformations [14]. The eigenvectors v; are found as solutions of (19) for each of the known eigenvalues, y;, m = 1 to 4. Since the magnitude of each of the eigenvectors is not uniquely defined by (19), we may arbitrarily fix one of the components. If the last component of each 7; is set to unity, for instance, then there results a 3 x 3 linear system to invert for each unique mindexed eigenvalue, Once the sets of eigenvalue-eigenvector pairs have been found within the nth layer we can represent the state-vector as a weighted summation of the eigenbasis functions which were postulated in (1 8  This expansion can be more succinctly expressed by employing matrix notation as shown below, where V, has its mth column formed from the components of 7; while U,(z) is a 4 x 4 diagonal matrix with nonzero elements given by exp We are now in a position to find the local transition matrix A, as defined in (17) where the thickness of the nth layer is d, =. 2, -z,-An alternative method for obtaining the local transition matrix in each layer is through a finite difference algorithm. Such a procedure is very straightforward to implement and circumvents the need to solve the eigenvalueeigenvector problem which was just considered. This method also has the potential for direct application to the case of continuously stratified anisotropic media.
Consider partitioning the region within the nth homogeneous layer into L, sublayers, each having equal thickness h, = d,/L,. A piecewise linear approximation for each of the four components of the state vector may be employed, as illustrated in Fig. 3. Enforcing the matrix differential equation in (4) The selection of the sublayer interval h, will depend upon the maximum expected rate of change of the state vector fields therein. For the case of the linear interpolate in (25), h, should be small enough to allow accurate tracking of the set of exponential basis functions in (18). A simple analysis shows that h, should be made small compared to the reciprocal of the magnitude of the largest (in modulus) complex eigenvalue, The power method can be used to form an accurate estimate of the magnitude of the largest eigenvalue in an efficient manner, without the need to solve for the smaller eigenvalues Once the local transition matrices are obtained for the various layers, the complete transition matrix can be computed as the product, v41.
This matrix is then used to find the scattering and transmission matrices, using the procedure described in the previous section.
The direct extension of the finite-difference algorithm to continuously stratified media in region 2 embodies the use of an approximate model composed of thin homogeneous layers. The layer thicknesses can be made locally dependent upon the material properties through estimates of maximum wavevectors. A localized recurrence formula, as in (26), will result from the finite-difference method and the' A-matrix can be formed by the product of successive local transition matrices. Enhancements of this approach can be made by implementing more optimized numerical methods such as that of Runge-Kutta [ 141.

COMP~~~ATIONAL VER~;ICATION
The algorithm was implemented on a microcomputer using the finite-difference approach, as described in the previous section, to 'compute the transition matrix. As input to this program, the user must specify the number and thickness of eacti' layer within region 2 followed by the entry of the complex matrix elements within each layer for the permittivity and permeability tensors. The material in region 1 is assumed to be free-space while that of region 3 can be specified as either free-space or a perfect electric conductor (PEC). Further inputs are required as to range and increments of both temporal frequencies and incident aspect angles.
Scattering computations are performed for two incident orthogonal linear polarizations, having the incident E-field either parallel or normal to the xz plane of incidence.
The computer program was checked in three steps. First, simple cases were run to verify that no obvious errors were present, while still exercising all of the computational steps of the algorithm. Next, several test cases involving complicated layered material configurations composed of lossy isotropic media were compared to the results of a previously developed computer program. This separate program, which used the 2 x 2 matrix formulation for isotropic material, had been extensively checked against independent results [ 111. The final validation was for a metal backed (PEC in region 3) homogeneous lossy anisotropic layer whose scattering could be computed by analytic means.
An example of the results from the initial set of computations involving free space in region 2 is shown in the pseudothree-dimensional plot of Fig. 4. In computing the transition matrix, region 2 was partitioned into multiple layers (based on the criterion discussed in the previous section) and the diagonal elements of the permittivity and permeability matrices within each layer were set to their respective free space values, with other elements set to zero. The independent variable axes are normalized frequency, kod and incident aspect angle (0" is normal to the surface). In this plot the case of a PEC in region 3 is considered which, of course, should yield complete reflection, I SI) I = 1. Identical results were obtained for S, while the off-diagonal elements of the scattering matrix were equal to zero. When region 3 was considered as free space, the computed S-matrix was found to be identically zero while the resultant transmission matrix had unit diagonals and zero off-diagonal elements.
Examples of the computational comparisons made between the general anisotropic program and the specialized isotropic program are shown in Figs. 5 and 6 . The first case is that of a single electrically thin layer of lossy dielectric bounded by a vacuum on both sides. Both SI1 and S2, are shown for a range of aspects and frequencies. The magnitude of SI, versus normalized frequency for 0" incidence on a lossy five-layer dielectric structure is shown4n Fig. 6 . In this case, region 3 is free space. Differences of the various validating computations were in the realm of a fraction of one percent in each case.
To check the accuracy of the general program for anisotropic material in region 2, a comparison was made to an A comparison of S-matrix magnitudes is given in Fig. 7 for a particular c&e of magnetoplasma coated metal where the complex permittivity matrix in (3 1) has nonzero elements E , = 100 €0, q, = 50 EO and E , = EO. Accuracy to four decimal places was observed in these computations. Further comparisons are pending, with use of experimental scattering data for multilayered magnetic materials. It is unlikely, however, that such results will be available for publication in the open literature.

CONCLUSION
A general analytic formulation and numerical method has been described for scattering and transmission problems involving generally lossy stratified anisotropic material. The method employs a 4 X 4 transition matrix approach combined with impedance boundary conditions at the entrance and exit planes of the stratified region. A new finite-difference approach is employed for field solutions within homogeneous layers of the structure. Scattering and transmission matrices are obtained which include full polarization information for a range of both frequency and aspect angle.for incident plane wave fields.
Extensive accuracy checks were performed on the general numerical algorithm [IO]. Some of these results are shown ' here for three levels of validation. This included comparisons with results from an earlier program which was designed to handle only isotropic media, using Ruck's technique [ 1 11. An anisotropic case was considered, using an independent analyti-cal solution for scattering by a layer of magnetoplasma covering a planar metallic surface.
A continuing effort is underway for incorporating the algorithm described here into an optimization program for use in the design of transmission and/or reflection structures. Input specifications include bounds on the characteristic behavior of the scattering and/or transmission matrix elements for a range of incident aspect angles and time-harmonic frequencies. Future work will address the analysis of electromagnetic interaction with two-and three-dimensionally inhomogeneous anisotropic material structures using a finite element formulation.

APPENDIX
The elements of the r-matrix found in (4) are given by rll =jk, -