010106-2-T THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING Radiation Laboratory ON ELECTROMAGNETIC FIELD PROBLEMS IN INHOMOGENEOUS MEDIA Technical Report By A. Mohsen Grant NGR 23-005-477 January 1973 'I) 2~ Prepared for: NASA Scientific and Technical Information Facility P. O. BoK 33 College Park, Maryland 20740 10106-2-T = RL-2195 Ann Arbor, Michigan 48105

010106-2-T ABSTRACT Analysis of electromagnetic fields in inhomogeneous media is of practical interest in general scattering and propagation problems and in the study of lenses. For certain types of inhomogeneities, the fields may be represented in terms of two scalars (components of Hertzian potentials). In a general orthogonal coordinate system, these potentials satisfy second order differential equations. Exact solutions of these equations are known only for a few particular cases and in general, an approximate or numerical technique must be employed. The present work reviews and generalizes some of the main methods of attack of the problem. The results are presented in a form appropriate for numerical computation. i

010106-2-T TABLE OF CONTENTS ABSTRACT I INTRODUCTION 1 II INTEGRAL EQUATION APPROACH 3 m DIFFERENTIAL EQUATION APPROACH 6 A: Auxiliary Impedance Function 6 B: Auxiliary Admittance Function 7 C: Auxiliary Reflection Function 7 D: Auxiliary Phase Shift Function 8 IV CONCLUSIONS 11 REFERENCES 12

010106-2-T I. INTRODUCTION The problem of wave propagation in inhomogeneous media is important in many situations in which electromagnetic waves are caused to propagate through media whose index of refraction has arbitrary spatial variations. Typical examples arise in radar signal intercepting ionized wakes, radio waves entering the ionosphere, microwave diagnostics of ionized gas, inhomogeneously-loaded waveguides as well as in the study of lenses. Electromagnetic fields in these inhomogeneous, and more generally anisotropic, media may be represented in terms of two Hertz vectors [1]. Since the solution of Maxwell's equations is uniquely determined by assigning initial values to the components of the electric field (E) and magnetic field (H) at t 0 which in the meantime must satisfy V. e E = V. * p T= 0, where p and c are the permeability and the permittivity, respectively, the degree of freedom in essence is reduced to four arbitrary functions [2]. Employing the gauge transformation, one may represent the field in terms of two components of either the electric (5T ) or magnetic (7 ) e m Hertz vectors or in terms of one component of 5 and another of F. The dife m ferential equations satisfied by these components are in general of higher order than the second, but under certain conditions they reduce to second order [3]. If the medium parameters (p and e) satisfy some further restrictions, the method of separation of variables may be employed to reduce the problem to a set of second order ordinary differential equations. The objective of the present work is to study the main available methods employed in order to solve such equations. The emphasis will be on the method of formulation rather than on different approximate or asymptotic techniques. Exact solutions of the equations encountered in dealing with inhomogeneous media are known for only a few particular profiles. In reference to this the important contributions of Heading [4] and Westcott [5] - [9] are worth men 1

010106-2-T tioning. Since only a few cases yield themselves to exact solutions, numerical or approximate methods are expected to play an important role. This enhances rather than diminishes the importance of those problems for which exact solutions can be found since exact solutions are often useful as starting points for approximate calculations. In addition, they may also help to establish limits of validity for various approximation methods. Even for cases when the exact solution is known, computation of such functions may not be the most economical way to find the solution. The two approaches used to solve the equations are the integral equation and the differential equation approach. Some of the basic computational aspects related to these two approaches are considered next. 2

010106-2-T II. INTEGRAL EQUATION APPROACH The differential equation encountered in the present study has the form 0" + P(x) '+Q(x) = 0 (1) which is valid for the range a < x < b, where the prime indicates derivative with respect to x; P(x) and Q (x) are functions of;p and c. It may be pointed out that the term involving 0' may be eliminated upon making the substitution 1 - P(rM)drj. 0 = 0 e. The function Q (x) may be written in one of the following forms Q(x)= q(x)-T(x) = qav T(x) (2) q (x) ap where q corresponds to the homogeneous case, q is the average value of Q(x) over the interval a <x <b and q (x) represents a convenient approxiap mation of Q(x). Obviously, of all these choices of q(x), the last is the most appropriate as far as the rate of convergence of a perturbation technique is concerned. T (x) may be considered as a perturbation function. The solutions corresponding to T (x) = 0 are assumed to be in terms of known functions which are denoted by f(x) and g(x). By a proper choice of f and g and using Green's function, the solution of (1) and (2) takes the form 0(x)= cf(x)+ XI V(s) K(x,s) 0(s)ds (3) where c is an unknown constant to be determined from the boundary conditions, 3

010106-2-T X (expansion parameter) and V(s) are related to the Wronskian of the homogeneous equation (unperturbed) and the Kernel has the form rf(x) g(s) x>s K(x, s) = (4) f(s) g(x) x < s Usually the Born approximation is used to solve the resultant Fredholm Integral Equation by means of a Neumann series in X [10] - [l]. This perturbation expansion has a limited radius of convergence determined by the lowest eigenvalue of the homogeneous equation [12]. Following Drukarev [12], if the function x W+X f(x) V(s) g(s) 0(s) ds (5) is added to (3) and 0(s) is written in the form 0(x)= N U(x) (6) then U (x) satisfies the Volterra equation U(x)= f(x)- V(s) f (x) g(s)- g(x) f(s) U(s) ds (7) and N= c.1-a V(s)g(s)U(s) ds ol (8) If the Neumann's series is used to solve (7), the solution converges for all x and an expression for 0(x) may be obtained as the ratio of two power series in x [12]. Brysk [12] has shown that Drukarev solution coincides with the 4

010106-2-T determinantal solution of the Fredholm equation. Thus, an iterative (successive approximation) solution of the Volterra equation generates the Fredholm determinants, just as an iterative solution of 0(x) generates the Born series. In other words, the effort normally expended on the Born approximation suffices to obtain the Fredholm solution with the Drukarev transformation. It may be pointed out that a numerical solution of (7) is much simpler than (3) in that U(x) is expressed in terms of its value for s no greater than x, thus starting at x a a, the solution can be developed successively to larger x [13]. 5

010106-2-T mI. DIFFERENTIAL EQUATION APPROACH Equation (1) may be solved, using approximate or numerical techniques, subject to the appropriate boundary conditions which 0 (x) must satisfy. The problem may be reduced to solving a set of coupled first order differential equations satisfied by the field components and an appropriate method, e.g., finite difference, may be used to get the solution [14]. A formal solution to the problem may be found if the original equation is reduced to Hill's equation [15]. In some situations, e.g., the slab problem, step by step integration, starting from an initial value, of the field equation is possible [16]. In general, it is more appropriate to introduce auxiliary functions in order to solve (1). These auxiliary functions are the auxiliary impedance, admittance, reflection and phase shift functions. The purpose of this section is to cast the differential equations satisfied by these functions in an appropriate form which is convenient for computation. A: Auxiliary Impedance Function: Let =0r (9) then, using (1), v satisfies 2 V +v2 +Pv+Q = 0 (10) The term including v may be eliminated and thus reducing the computational effort if we introduce Z- +-P (11) 0 2 Then, Z satisfies the equation z + z2 =- p2+ 1 pt -Q (12) 4 2 6

010106-2-T B: Auxiliary Admittance Function: Following the same ideas as for the auxiliary impedance function, an admittance function may be introduced in the form Y=Q. + 2 T (13) where Y is given by y,-y2=Q 1 T + 1 T2 (14) and T = P + Q (15) Q Due to the appearance of 0 in the denominator of Z, Z = o whenever 0 = 0. Similarly, Y is infinite if 0' is equal to zero. In order to get around this difficulty, Garbacz l17] suggested the use of the equations satisfied by Z whenever Y tends to infinity and vice versa. A more convenient way is through the use of an auxiliary reflection function. -i C: Auxiliary Reflection Function Introducing the function R= ~ 0 +w (16) where i = atv --- exp [-(Q -1) dx] =, the function R then satisfies R' +R2 2(1+P+Q) -U'+ 2 (17) Contrary to Z and Y, R is bounded as may be observed from (16). It may be noted that all the auxiliary functions introduced satisfy similar Ricatti type equation whose right hand side is a known function of position. 7

010106-2-T D: Auxiliary Phase Shift Function: The phase shift concept extensively used in Quantum Mechanics [18], was originally introduced in dealing with potential scattering from cylindrical or spherical regions. If the solution is presented in the form 0(x) = A (x) cos 6(x) J (kx) + sin 6 (x) Y (kx) (18) n n where 6(x) is unknown function, the asymptotic behavior of 0(x) is proportional to cos kr - (n+) - - 6j Hence 6(x) was called the phase shift. Integral and variational expressions [18] as well as upper and lower bounds of this function are known in the literature [19]. Recently, Shafai [20], [21] derived the differential equation satisfied by the phase shift due to scattering from cylindrical as well as spherical regions. The purpose of this section is to extend Shafai's work and to solve equation (1) using the phase-shift approach. Using Lagrange's method [22], a solution of (1) may be written in the form 0(x)= c(x) f(x) + (x) g(x) (19) where a' f + ' g = 0 (20) a'l f' + 1' gt - T (21) Equations (20) and (21) lead to the evaluation of a(x) and 13(x). Introducing the amplitude function A (x) and the phase function 6 (x) given by a (x) A A (x) cos 6 (x) (22) (x)= A (x) sin 6 (x) (23) and using (20) and (21), one may show that 6(x)= - W ( g) os6f+sin6g2 (24) 8

010106-2-T A'(x) = -A(x) T() [cos6f + sin 6 g] [sin6 f - cos6g] (25) A2 W fg) where W(f, g)= f g' -g f is the Wronskian. Hence, A(x) and 6(x) satisfy first order differential equation and notably the differential equation satisfied by 6(x) is independent of A(x). The monotonic characteristic of 6(x), as a function of the perturbation T (x), may be noted from equation (24). In order to solve equation (24) for 6 (x), an initial value is needed which may be evaluated from the boundary condition on O(x). Assuming that such condition requires a0 +b0=0 atx=x (26) it follow, upon using (20), that af' (x1) + bf(x1) tan 6(x1) - (27) 1 ag'(x1)+b g(x1) On the other hand, if g tends to infinity as x -- x1 while the field is assumed to be finite there, an initial value of zero may be taken for 6 (x). Once an initial value of 6(x) is known, equation (24) may be easily solved numerically using, for example, the Predictor-Corrector of the Fourth order RungeKutta method. Once 6 (x) is known, A(x) is then given by A (x) WA exp { - (cos 6f + sin 6g) (sin6f - cos 6g) dW (28) A 2x ) W-f, ge 2 O where A is the value of A (x) at x = x which is assumed to be known from o o the boundary conditions. Assume that we are dealing with two regions whose boundary is at x = x and their corresponding solutions are 0 9

010106-2-T 0= a f(x) + 3 g(x) in region I (x > x ) (29) o 0 0 and =A [cos6f(x)+sin6g(x) in region (x < x ) (30) where region I is assumed to be homogeneous (unperturbed). If the boundary condition at x = x requires the continuity of both 0 and 0', it follows, using 0 (20), that a A(x ) cos6(x ) (31) and /3 =A(x ) sin6(x ) (32) Hence if a is known, usually related to the incident field, as well as 6(x ), 0 0 A (x ) is obtained directly from (31). It may be noted that 1o may be written in terms of 6 (x ) only (P3 = a, tan 6 (x )) which leads to the significant conclusion that only the solution for 6 (x) is needed if we are interested only in the field in region I. 10

010106-2-T IV. CONCLUSIONS The present work deals with some of the computational aspects of the two basic approaches used to solve a second order ordinary differential equation. These approaches are, the integral and the differential equation techniques. In cases when the exact solution is unknown, the integral equation technique is more appropriate if an analytical expression via perturbation is attempted. Generally speaking, the mechanics of solving differential or integral equations by means of digital computers is such that it is preferable to solve non-linear initial value problem as opposed to linear problems with two-point conditions. The reason for this is that the initial value problem can be resolved by means of a simple iteration procedure, which is ideally suited to digital computers, whereas the two-point boundary problem requires the solution of a large system of equations. In addition, a test of the convergence of the solution for the initial value problem is much easier. Another advantage is that a solution may be obtained for intermediate dimensions of the range considered, whereas in two-point boundary formulation the solution must be carried out for each range. For some particular cases, e. g., the slab problem, the reduction to an initial value problem is easy but in general an auxiliary function should be introduced to achieve this purpose. Since they are, in general, well behaved and bounded functions, the reflection coefficient and the phase shift functions are the most appropriate for numerical evaluation. If an initial value is known, it may be convenient to use the reflection coefficient function since it does not require computation of the homogeneous solutions f(x) and g (x) which are needed for the phase shift approach. In general, an initial value for 6 (x) is easier to evaluate since the equation for 6 (x) is independent of the amplitude function. 11

010106-2-T' REFERENCES 1. A. Nisbet, "Electromagnetic potentials in a heterogeneous non-conducting medium", Proc. Roy. Soc., vol. A 240, pp. 375-381, 1957. 2. R. K. Luneberg, "Propoagation of Electromagnetic Waves from an Arbitrary Source Through Inhomogeneous Stratified Atmospheres", New York University, Math. Res. Group, Res. Rept. No. 172-6, Jan. 1948. 3. A. Mohsen, "Electromagnetic Field Representation in Inhomogeneous Anisotropic Media", to be submitted for publication, IEEE Trans. Antennas and Propagation. 4. J. Heading, "Refractive index profiles based on the hypergeometric equation and the confluent hypergeometric equation", Proc. Camb. Phil. Soc., vol. 61, pp. 897-913, 1965. 5. B. S. Westcott, "Exact solutions for electromagnetic wave propagation in spherically stratified isotropic media", ibid., vol. 64, pp. 227-235, 1968. 6. B. S. Westcott, "Exact potential functions in spherically stratified media", ibid., vol. 64, pp. 1089-1098, 1968. 7. B. S. Westcott, "Exact solutions for the electromagnetic wave propagation in a cylindrically stratified axially magnetized plasma", ibid., vol. 66, pp. 129-143, 1969. 8. B. S. Westcott, "Exact solutions for vertically polarized electromagnetic waves in horizontally stratified media", ibid., vol. 66, pp. 675-684, 1969. 9. B. S. Westcott, "Exact solutions for vertically polarized electromagnetic waves in a horizontally stratified magneto plasma", ibid., vol. 67, pp. 491-501, 1970. 10. J. Heading, "Composition of reflection and transmission formulaei', J. Res. NBS, vol. 67D, No. 1, pp. 65-77, 1963. 11. L. A. Pipes, "A perturbation method for the solution of linear matrix differential equations", J. Franklin Inst., vol. 283, No. 5, May 1967. 12

9 010106-2-T 12. H. Brysk, "Determinantal solution of the Fredholm equation with Green's function kernel", J. Math. Phys., vol. 4, no. 12, pp. 1536-1538, Dec. 1963. 13. H. Brysk and M. L. Buchanan, "Scattering by a cylindrical Gaussian potential: exact solution", Can. J. Phys., vol. 43, pp. 28-37, Jan. 1965. 14. L. Shafai, "Electromagnetic fields in the presence of cylindrical objects of arbitrary physical properties", ibid., vol. 48, pp. 89-98, 1970. 15. K. F. Casey, "Application of Hill's functions to problems of propagation in stratified media", IEEE Transactions Antennas and Propagation, vol. AP-20, pp. 368-374, May 1972. 16. J. H. Richmond, "Transmission through inhomogeneous plane layers", IRE Transactions Antennas and Propagation, vol. AP-10, pp. 300 -305, May 1962. 17. R. J. Garbacz, "Electromagnetic scattering by radially inhomogeneous spheres", Antenna Lab., Ohio State University, Rept. No. 1223-3, 1962. 18. L. Schiff, Quantum Mechanics. New York: McGraw-Hill Book Co., Inc., 1965. 19. T. Kato, "Upper and lower bounds of scattering phases", Progress of Theoretical Phys., vol. VI, No. 3, pp. 394-407, May-June, 1951. 20. L. Shafai, "Scattering by cylindrically symmetrical objects, method of phase and amplitude functions", Int. J. Electronics, vol. 31, no. 2, pp. 117-125, 1971. 21. L. Shafai, "Scattering by spherically symmetrical objects", Can. J. Phys., vol. 50, no. 8, pp. 749-753, 1972. 22. E. L. Ince, Ordinary Differential Equations. New York: Dover Publ., Inc., 1956, p. 122. 13