RL-878 PROPAGATION OF ELECTROMAGNETIC PULSES THROUGH PLANAR STRATIFIED MEDIA - A FINITE DIFFERENCE APPROACH Kasra Barkeshli April 1990 RL-878 = RL-878

RL — 8 7 8 Propagation of Electromagnetic Pulses through Planar Stratified MediaA Finite Difference Approach Kasra Barkeshli Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109 April 25, 1990

Abstract The propagation of electromagnetic fields through a one dimensional inhomogeneous medium is studied by solving the Maxwell's equations with a numerically nondissipative finite difference scheme. The medium is assumed to be stratified in layers of specified permittivity, permeability, and conductivity. In particular, the reflection from and transmission through composite slabs of finite thicknesses is studied for a normally incident Gaussian pulse.

Propagation of Electromagnetic Pulses through Planar Stratified MediaA Finite Difference Approach Kasra Barkeshli Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109 April 25, 1990 Abstract The propagation of electromagnetic fields through a one dimensional inhomogeneous medium is studied by solving the Maxwell's equations with a numerically nondissipative finite difference scheme. The medium is assumed to be stratified in layers of specified per.mittivity, permeability, and conductivity. In particular, the reflection from and transmission through composite slabs of finite thicknesses is studied for a normally incident Gaussian pulse. 1 Introduction Although a direct solution of Maxwell's equations via finite differences was proposed as early as 1966[1], a serious application of the method to general electromagnetic scattering problems was not considered until 1975[2]. Traditionally, the scattering problems are solved by the complex time harmonic approach which is based upon the assumption of an eiwt time dependence in the fundamental Maxwell's equations. The solutions are obtained at 1

single frequencies which in principle may be combined using spectral methods to yield the final result. On the other hand, in the study of the transient response of scatterers, it may be more efficient and instructive to solve the time dependent form of the equations directly, using finite differences. This approach is known as Finite-Difference Time-Domain method(FDTD) in the electromagnetics community[2, 3]. In this paper, an example of this method is given for the case of propagation of an electromagnetic pulse in planar stratified medium (Figure 1). Starting from Maxwell's equations in one dimension, a consistent leapfrog scheme is developed for the solution of field quantities in the medium. The stability of the scheme and its compliance with the Courant-Friedrichs-Lewy criterion is discussed next. Since the scattering problem under study represents an open problem with an unbounded physical domain, a method has to be devised to enable a proper truncation of the domain of computation. This is done by the incorporation of a first order absorbing boundary condition which accurately simulates the medium surrounding the outer surface of the computational domain. Finally, several numerical examples are presented to display the performance of the method. The treatment presented here is similar to that given in [4] although we will take a different approach. In particular, the stability criterion and the absorbing boundary conditions are discussed in some detail. 2 Maxwell's Equations The propagation of electromagnetic waves in space is described by the Maxwell's equations which form a system of coupled differential equations for electric and magnetic field vector intensities, E and H. The Maxwell's equations in a source free medium are[5] V xE -- V D= 0 at (1) OD VxH= -+J V B=O at where D and B are the electric and magnetic flux densities, respectively and J is the volumetric conduction current density. They are related to the field 2

quantities through the constitutive relations and for an isotropic medium D = E B = H (2) J = E where c, /, a are assumed to be time-independent functions of space denoting permittivity, permeability, and conductivity of the medium, respectively. From Maxwell's equations, it may be shown that the fields satisfy the inhomogeneous wave equation in the lossy medium. Considering the electric field for example, we have after combining (1) and (2) 1 a2 a (V2- - 2)E + V(ln i) x (V x E) + V [E V(ln e)] - at-E = O (3) where c = 1/v/if is the speed of light in the medium. Assuming the medium is stratified in the x direction, the Maxwell's equations may be specialized to the one dimensional case. Retaining only the Ey and Hz components, we have 9Ey I H (9xHz at - (x) ax + -()E (4) aH, 1 aEy At Pi(x) Ox and the wave equation (3) reduces to 02 12 p 2 ( o p)EY -)E (- + _ 2 C2 at 2 Itax 0t where ' denotes differentiation with respect to x. 3 The Finite Difference Scheme In this section we give a suitable finite difference scheme for the solution of (4). We will first treat the propagation in a lossless medium (cr = 0). It 3

is noted that in this particular case, (4) represents a hyperbolic system of flux-conservative equations with variable coefficients as t( Ey/ 0 O- da EY (6) at (Hz - -1/p 0 ) Ax Hz) which may be put into the standard vector form a a -(u) = A(x)x (u) (7) Since we are interested in the propagation of the waves for a relatively long periods of time, a numerically non-dissipative scheme is preferable in this case. Therefore, considering a leapfrog difference scheme1, we have n+l _ En-l -(Hn -H n-) (8) H 1 = Hn ---(En+- En1) where as usual Fn = Fy. Ax, n At) and we have set A = At/Ax. 3.1 Consistency The leapfrog method is second order accurate in both space and time andunlike first order methods (such as Lax-Friedrichs) which require significantly small time steps-is limited only by the stability considerations. 3.2 Stability In order to establish the stability of the method for the general case, we may discuss the stability of the constant coefficient form and extend the result to that of variable coefficients. 1In the leapfrog scheme, the space and time derivatives are approximated by central differences. 4

3.2.1 Constant Coefficient Case From (7), the numerical scheme in this case takes the form Vn+'(x) = vn-1(x) + AADov"(x) (9) where the matrix A is given in (6) and Do denotes the central difference operator and we may write ( vn+1) = 2tADo I )( Vn ) (10) V - n I 0 ( n0-1 (0) n+l = QoWn. (11) Applying the von Neumann analysis, we take the discrete Fourier transform of both sides of (11) to obtain Wn+l = QWn (12) where Qo denotes the symbol of the scheme given by o 0 -2iA/c sin I Qo0 ( -2iA/sin 0 (13) I 0/ where w = wAx and w is the wave number. In order to check the condition on the power-boundedness of Qo we form the characteristic equation X4- 2[1 - 2(cA sin)2]X2 + 1 = 0 (14) and find that X2C- i 1 - (15) where = 1 - 2(cAsin )2. It is clear from (15) that IXl = 1 if I11 <' 1 or, equivalently, if cA < coA< 1 (16) where co is the speed of light in vaccum. As expected, condition (16) is compatible with the Courant-Friedrichs-Lewy (CFL) criterion for stability. 5

3.2.2 Variable Coefficient Case The numerical scheme is now of the form Wn = Qo(x)W7 (17) where Qo() = (2AtA(x)Do ) (18) To establish the stability of the scheme in this case, we employ the following theorem: Theorem: The stability of the leapfrog scheme with variable coefficients vn +(x) = v-l(x) + 2AtPo(x, t)vn(x) (19) is compatible with the CFL criterion for the corresponding frozen coefficient case if Po is anti-symmetric, i.e. if (Por, s) = -(r, Po) V r, s (20) Clearly, the difference operator Po = A(x)Do is not anti-symmetric. However, it may be easily shown that it is located within a bounded distance from an anti-symmetric operator given by Po(X) = 2[A(x)Do + Do (A(x).)]. (21) Incorporation of this new operator into our formulation leads to the following scheme E -+ = E - ( H1 + H-1 - (H+l-H )- 1 H ) (22) + Hnl ( 1 (E EP) 1 /En V 2 v+1V+1 YVV+ V YV1 - 6

The nonzero conductivity case (a $ 0) is considered next. Referring to (4), a simple extension of the previous scheme is to replace the first of equations (8) by +1 n- A _2At E+1 En-1 — (H+1 - Hun_)...E (23) and leaving the second equation unchanged. Unfortunately, it was found that for nontrivial conductivity values(cr t 1 S/m) the stability of this scheme is questionable for moderate cA values. This is because the stability condition (16) was derived for negligible losses and the introduction of the new term clearly violates the condition. A solution to this difficulty is to write the loss term in (23) implicitly as E+l = 1 E - A (Hn - H ) - 2At + so that EV+1 [eEV - A(H+1 - H-1)] /(e + 2Atac) (24) Hn+1 Hn-1 A (Ev+ - n) A stability analysis of the new scheme shows that the new term in (24) introduces numerical dissipation, thus making the scheme inherently more stable for nonzero or. It is, of course, clear that this modification does not affect the results obtained from (8) for the lossless case. 4 Incident Field For the incident field we may consider a Gaussian electric pulse of the form Ey(x, t) = Eoe-P( —ct)2 (25) where Eo is the peak amplitude and Xp denotes the initial position of the pulse peak. The coefficient / specifies the shape (and equivalently, the bandwidth) 7

of the pulse. For a pulse of amplitude Eo that decays to aEO in t, seconds (a < 1), d is given by -Ina (26) (ctw)2 (6 Writing (25) as a difference equation we obtain for the two field components E = Eoeln a[(v - p - cAn)/k]2 (27) n Eo eln a[(v - p - cAn)/k]2 (28) where p = xp/Ax, k = ctw/Ax and Z = a/C is the local intrinsic impedance of the medium. The above expressions for the incident fields are used to obtain the initial data at n = 0, 1 as called for in (8) and (22). For this purpose, the pulse is initially positioned at a node where it is not disturbed by the presence of layers whose properties are different from that of the background. 5 The Absorbing Boundary Conditions The problem under consideration is an example of open problems in which the domain of computation is not bounded by physical boundaries. For this reason, an artificial boundary has to be introduced in order to simulate the propagation of the outgoing waves in the medium surrounding the outer surface of the domain of the computations. The new conditions imposed at these boundaries are called absorbing boundary conditions2. Such a boundary condition not only has to minimize the artificial reflections which occur at the boundaries but also has to guarantee stable difference approximations. Engquist and Majda[6] give a systematic method for obtaining a hierarchy of such boundary conditions. Also, Bayliss and Turkel[7] derive another type of conditions based on the Wilcox expansion. Here, we employ the boundary conditions proposed by Engquist and Majda and further discussed by Mur[8] for electromagnetic problems. 2Also referred to as the radiation conditions. 8

Considering the governing wave equation (5) and assuming, for the moment, that the losses are negligible in the vicinity of the boundaries, we write 02 1 a2 ) E y ( OsX2 C2 9t2 t aX We now seek conditions on the boundary, d0 such that the solution to the above equation is consistent with an outgoing wave propagating in the extended medium of the same electrical properties as that of the boundary neighborhood. This implies,' = 0 in the vicinity of the boundary. For the family of solutions W = f(x z ct) to the (homogeneous) wave equation traveling toward the boundaries, the first order condition (f- +-)W = (30) OX c't Of would determine a W on the outer surface which is outgoing, is not reflected, and is thus absorbed. Although this is a first order condition, it actually provides zero reflection for the waves normally incident on the boundary surface[6]. It is also interesting to note that neglecting the effect of the conductivity term o is of second order and does not deteriorate the order of accuracy of (30) as discussed in [6]. A finite difference approximation to the above boundary condition is given by 1 Dr(Wo ) —D(Won) = 0 C 1 C yielding Eo+1 = Eo +cA(E1 - E) (31) E+1 = EM-cA(E -EM_1) (32) where v = 0 and v = M denote the left and right boundaries, respectively. 6 Numerical Results An interactive computer program was developed which implements the preceding formulation. The user specifies the mesh size, pulse width, the number 9

of slabs present and their respective locations and electrical properties. The program calculates the field quantities at specified time and space increments and displays the results in both static and animated forms. Figure 2 shows the reflection from and transmission through a 9 cm thick purely dielectric ([t = 1, cr = 0) slab of relative permittivity er == 4 by an Ey-polarized 400ps (to 0.1% amplitude, alpha = 0.001) Gaussian pulse propagating in the free space. The plots correspond to different time steps. Figure 3 shows the same results in a space-time plot. The reflection and transmission of waves are clearly visible. Also, it is seen that the waves reaching the outer boundaries are absorbed with minimal reflection. As a further example, the same configuration is examined with a slab of permittivity 4 and conductivity 0.8 S/m. The results are shown in Figures 4 and 5 where the dispersion of the pulse due to the presence of losses are evident. Finally, Figure 6 shows the propagation of a pulse through a lossless nonreflective slab (Er = ulr) with a normalized intrinsic impedance of unity. 7 Summary The Finite-Difference Time-Domain technique is an efficient method of solving transient electromagnetic problems. In this paper, the Maxwell's equations were directly solved by a second order leapfrog scheme and successfully applied to the problem of propagation of electromagnetic pulses in an stratified medium. A first order absorbing boundary condition was employed to minimize reflections from the artificial boundary introduced by truncation of the domain of computation. Further study of this problem would include a spectral analysis of the solutions obtained when using various waveforms of different bandwidths as well as employment of higher order absorbing boundary conditions to improve the accuracy of the solutions. Also, the more complicated problem of electromagnetic propagation in time-varying inhomogeneous media whose electrical properties are functions of both space and time may prove amenable to finite difference approaches in an attempt to model transient responses of targets in plasma. 10

References [1] K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. Antenna Propagat., AP-14(3), May 1966. [2] A. Taflove, M. E. Brodwin, "Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell's equations," IEEE Trans. Microwave Theory & Tech., MTT-23(8), Aug. 1975. [3] A. Taflove, "Application of the finite-difference time-domain method to sinusoidal steady state electromagnetic-penetration problems," IEEE Trans. Electromag. Compat., EMC-22(3), Aug. 1980. [4] R. J. Luebbers, K. S. Kunz, and K. A. Chamberlin, "An interactive demonstration of electromagnetic wave propagation using time-domain finite Differences," IEEE Trans. Education, Vol. 33, No. 1, pp. 60-68, Feb. 1990. [5] J. Stratton, Electromagnetic Theory, McGraw-Hill, New Jersy, 1941. [6] B. Engquist and A. Majda, "Absorbing boundary conditions for the numerical simulation of waves," Math. Comp., Vol. 31, pp. 629-651, July 1977. [7] A. Bayliss and E. Turkel, "Radiation boundary conditions for the numerical simulation of waves," Comm. Pure Appl. Math., Vol. 33, pp. 707-725,1977. [8] G. Mur, "Absorbing boundary conditions for the finite-difference approximation of time-dependent electromagnetic-field equations," IEEE Trans. Electromag. Compat., EMC-23(4), Nov. 1981. 11

List of Figures Figure 1: Geometry of a stratified medium. Figure 2: Reflection from and transmission through a dielectric slab by a 400ps unit amplitude Gaussian electromagnetic pulse. Slab thickness 9cm, or = 4, ir = 1, a = 0. Spatial increment Ax = 1.5mm, coA = 0.5. Figure 3: Propagation of the pulse of Figure 2 in space and time. Figure 4: Reflection from and transmission through a dielectric slab by a 400ps unit amplitude Gaussian electromagnetic pulse. Slab thickness 9cm, er = 4, jri = 1, a = 0.8S/m. Spatial increment Ax = 1.5mm, co0 = 0.5. Figure 5: Propagation of the pulse of Figure 4 in space and time. Figure 6: Transmission of a Gaussian pulse through a lossless nonreflective dielectric slab.

(E, ) 0 Figure 1: Geometry of a stratified medium.

- 1.0 0.5 r0.0 -0.5 -0.5. n=210 t=525 ps 0.5 > 0.0 -Pi -0.5 -1.0- '. i --. 0 30 60 90 x, cm -1.0 X, cm 1. - 0.5. 0.0 -0.5 --1.0- I. I I 0 30 60 90 x, cm 1.0-. --.- - --. 0.0 ft -0.5 --1.0- I ' I - 0 30 60 x, cm 90 -- Figure 2: Reflection from and transmission through a dielectric slab 400ps unit amplitude Gaussian electromagnetic pulse. Slab thickness r = 4, Tp = 1, a = 0. Spatial increment Ax = 1.5mm, CoA = 0.5. by a 9cm,

0.5 Ea -0.5 -1.0 0.5 0.0~ -0.5, -1.0 XI cm X. cm 1.0- n=630 t=1.-575 0.5 ->0.0 - -0. 5 -1.0 x. cm Figure 2. (ontinued

FDTD SOLUTION OF EM WAVES Propagation of EM waves ( 1, 120, 2) 22.0 17.6 13.2 t 8.8 4.4 0.0 0.0 18.0 36.0 54.0 72.0 90.0 x 04/25/90 12:41:57 Figure 3: Propagation of the pulse of Figure 2 in space and time.

- - 0.s5 0.0...... -0.5 -1.0 - ', I ' - - ' 0 30 60 90 X, cm _ __ 1.0 0.5. 0.0 -0.5 -1.0 1.0 0.5. 0.0 w -0.5 -1.0 n=540 t=1.350 us DI",, I, - ' 0 30 60 90 x, cm x, cm Figure 4: Reflection from and transmission through a dielectric slab by a 400ps unit amplitude Gaussian electromagnetic pulse. Slab thickness 9cm, e, = 4,,r = 1, o = 0.8S/m. Spatial increment Ax = 1.5mm, coA = 0.5.

14.0 11.2 8.4 t 5.6 2.8 0.0 FDTD SOLUTION OF EM WAVES Propagation of EM waves ( 1, 34, 1) - - ~,s I~ I,, I I, I - I 0.0 18.0 36.0 54.0 72.0 90.0 i x 04/25/90 15:43:34 Figure 5: Propagation of the pulse of Figure 4 in space and time.

FDTD SOLUTION OF EM WAVES Propagation of EM waves ( 1, 40, i) 20.0 16.0 12.0 t 8.0 4.0 0.0 0.0 18.0 36.0 54.0 72.0 x 90.0 04/25/90 15:52:25 Figure 6: Transmission of a Gaussian pulse through a lossless nonreflective dielectric slab.