INVESTIGATIONS ON EXCITATION AND PROPAGATION IN IONIZED MEDIA Chiao-M. Chu, John LaRue, and David B. vanHulsteyn 6663-1-F = RL-2143 This document is subject to special export controls and each transmittal to foreign governments or foreign nationals may be made only with prior approval of RADC (EMLI), GAFB, N.Y. 13440. AFLC, GAFB, N.Y., 8 Aug 66-130

FOREWORD This report was prepared by the University of Michigan, Dept. of Electrical Engineering, Ann Arbor, Michigan; under Contract Number AF30( 602) -3381; System No. 760K, Project No. 5579, and Task No. 557902. The RADC project engineer is Leonard Strauss, EMASA. This is a final report covering the period of work from June, 1964 to October, 1965; the originator's report number is 6663-1-F. Release of subject report to the general public is prohibited by the International Traffic In Arms Regulation. This technical report has been reviewed and is approved. Approved/.- v/,f *5h/ LEONARD STRAUSS Contract Engineer Approved: BOND>77 J ~-TIiOMAS S. BOND, JR. Colonel, USAF Ch, Surveillance and Control Division FOR THE COMMANDER: I:I NG JG BELMAN Chief, Acvanced Studies Group ii

ABSTRACT A study of the excitation and propagation of wave like disturbances in an ionized medium, such as the ionosphere, is made based on the linearized Euler's equation and Maxwell's equation. The local propagation constants of the basic modes of propagation are discussed. A computer program for the evaluation of these constants with given ionospheric properties is given. Methods of investigating the propagation of such waves in inhomogeneous and/or bounded media, such as ray tracing, invariant embedding, reflection and refraction, orthogonal expansion, and the use of a general matrix formulation are presented. A unified matrix-operator transform method for investigation of the excitation and propgation of disturbances in an ionized medium is proposed. iii

TABLE OF CONTENTS I INTRODUCTION 1 II THE EQUATIONS OF MOTION 3 2.1 The Boltzmann Equation 3 2.2 The Macroscopic Equations 5 2.3 The Perturbation Equations 10 III PROPAGATION CONSTANTS 15 3.1 The Dispersion Relation 15 3.2 The Basic Modes of Propagation 22 3.3 Propagation Constants in the Ionosphere 41 IV THE OPERATOR TRANSFORM METHOD 48 4.1 The Operator Approach 48 4.2 The Integral Equation 52 V EXCITATION IN A HOMOGENEOUS PLASMA 63 5.1 General Formulation 63 5.2 Collisionless Electron Plasma 66 VI SOURCES IN A BOUNDED PLASMA 73 6.1 Statement of the Problem 73 6.2 The Potential Functions 74 6.3 Vertical Oscillating Dipole Above a Perfectly Conducting Plane 85 6.4 Vertical Oscillating Dipole Above a Perfectly Conducting Sphere 95 6.5 Discussion of Result 107 VII INHOMOGENEOUS MEDIA 108 7.1 Introduction 108 7.2 Ray Tracing 109 7.3 Reflection and Refraction 112 7.4 Invariant Embedding 120 7.5 Operator Transform Method 132 VIII CONCLUSIONS 136 REFERENCES 137 v

APPENDICES A MACROSCOPIC PROPERTIES OF THE IONOSPHERE 139 A. 1 Discussion of the Model 139 A.2 Macroscopic Properties 139 B COMPUTER PROGRAM FOR NUMERICAL EVALUATION OF THE PROPAGATION CONSTANT 153 C MATRIX ANALYSIS 165 C. 1 Discussion of the Method 165 C. 2 General Formulation 165 C.3 Matrix Inversion 169 C.4 Physical Interpretation 170 C.5 Kron's Method 172 C.6 Graph Theory 172 C.7 Electrical Network Theory 175 C. 8 Application to the Excitation Problem 183 vi

EVALUATION The subject report represents about 1 1/2 years effort on excitation and propagation of disturbances in the ionospheric plasma. Towards describing these phenomena, the contractor has applied matrix-operator techniques and developed generalized formulations applicable to a wide range of associated problems. However, particular problems though noted and discussed are not resolved. For instance, excitation is formulated for all possible souces but problems are solved only for extremely simplified media/source conditions. The report significance is based on its scope rather than on detailed calculations. As indicated, the problem types are as broad as possible and the mathematical techniques powerful. However, the formulations are limited to linear macroscopic conditions. The significance of the contractor' s results is that while little new physics is introduced in this report, the mathematical bases for several important problem classes have been developed. LEONARD STRAUSS Contract Engineer vii

I INTRODUCTION This report presents the results of an investigation of the excitation and propagation of wave like disturbances in the ionosphere. Recognizing that the ionosphere is a complicated physical system such that the eventual solution of the problem has to rely heavily on approximate, numerical calculations, the primary efforts of our investigation were devoted to the choice of a model, the investigation of techniques that may yield the solution to the problem directly, and the solution of some idealized problems which may allow some insight into the validity of a particular approach or which may yield results useful to the eventual solution of the complete problem. The model of the ionosphere is taken as a three-fluid plasma containing ions, electrons and neutral particles whose density, temperature and average mass are known. The problem is formulated mathematically by linearizing the coupled hydrodynamic and Maxwell equations. In order to clarify the approximations involved in the use of these equations, they are deduced from the Maxwell-Boltzmann equation in Chapter II. The basic modes of propagation in the three-fluid plasma are discussed in Chapter III. It is found that if all the coupling effects are included, the propagation constants in a homogeneous medium satisfy a tenth order algebraic equation. The basic properties of these waves are introduced through the study of degenerate cases. Based on collected published data on the ionosphere (Appendix A), propagation constants were evaluated by use of a digital computer, for several discrete angular frequencies and directions 1

of propagation, corresponding to an altitude of 100Km. The results of the computations are presented in Section 3. 3. Due to a lack of time and funds, similar calculations were not made at other altitudes. The computer program used in the calculations is discussed in Appendix B. A generalized operator transform approach was formulated as a unified approach to the problem of the excitation and propagation of waves in ionized media. This method, valid in principle for both homogeneous and inhomogeneous media, reduces the mathematical problem to that of solving a generalized integral equation in transform space. For homogeneous media, analytical solutions are possible. Some new results on the solution of the excitation problem in a homogeneous electron plasma have been presented by Wu (1965) in a separate report. A general discussion of the application of this approach to excitation problems in inhomogeneous media is given in Section 7. 3. The excitation problem in a bounded two-fluid plasma is considered in Section VI, using the conventional methods of orthogonal expansions. The formal solution of the fields excited in the plasma by a dipole source above a conducting plane and sphere is obtained in a form suitable for numerical integration. Based on the assumption that the solution to the excitation problem in a homogeneous medium may be used near the source even in an inhomogeneous medium, the effect of the inhomogeneity on wave propagation is investigated in Section VII. Conventional ray theory, the technique of invariant embedding and the reflection and refraction of waves at a discontinuity, is formulated for an electron plasma. A novel approach, using Kron's method of large system analysis, is investigated in Appendix C. All of these methods, in principle, may be extended to more complicated problems. However, the feasibility of obtaining numerical solutions by these methods remains to be tested. 2

II THE EQUATIONS OF MOTION 2.1 The Boltzmann Equation Basically, the study of the excitation and propagation of disturbances in the ionosphere can be formulated in terms of the kinetics of a partially ionized gas mixture. From the microscopic point of view, therefore, one may start from the well known Boltzmann's equation with a source term, which is x a - a- f(r, u, t).+ u f+ f I= I(f., fk)+Sj(r,ut) r j mk a k (2.1) The terms in (2. 1) are defined as follows: 3 th f.(r, u, t)du is the average number of particles of the j component per unit volume in the velocity range (u, u+du); m. is the mass per particle of the j component; J I(fj, fk) Sj(L, t)du is the force (external anl long range) exerted on each particle of the j component; is the rate of change of f. due to short range interaction with the k component, and is known as the collision term; is the source per unit volume per second of particles of jth component externally injected into the mixture in the velocity range (u, u+du). 3

The long range forces of interaction between charged particles may be expressed in terms of the electric and magnetic field in the medium. Thus, the term X may be written as X x. +q(E+uxB) (2.2) where x. is the external mechanical force; q. is the charge per particle of the j component; E, and B are the macroscopic electric and magnetic field respectively. On the other hand, a detailed knowledge of the short range interaction, especially in an ionized media, is not known. Simple models for this interaction (or collision), are generally postulated in the investigation of Boltzmann's equation. The solution of (2. 1), even without the source, the external force and the collision term, is extremely difficult, and has been the subject of a great deal of research which, to date, has yielded only approximate solutions to some idealized problems. Since the pertinent macroscopic properties of the mixture, such as temperature, density, average velocity, etc., are all contained in the first three moments of Eq. (2. 1), the present work is confined to the investigation of the macroscopic equations (magneto-hydrodynamic equations) governing the laws of mass, momentum, and energy conservation. Mathematically, these *. We assume here that the effective electric and magnetic fields acting on each particle are essentially the macroscopic field, see Ginzburg (1961), p. 28. 4

equations, which are interpreted as the first three moments of Eq. (2.1), do not represent a complete description of the system. In general, approximate relations between the higher moments and the first three moments are necessary. These relations are derived either by using perturbation theory (Chapman and Cowlings, 1939), an approximate distribution function (Everett, 1962), or * variational methods (Marshall, 1957). The results obtained by these methods invariably introduce some macroscopic constants (or parameters) into the system of hydro-dynamic equations. For a complicated mixture, such as the ionosphere, some of these "constants" have never been determined experimentally or calculated theoretically with certainty. For an engineering approach, we shall, in this work, simplify these relations so that the "constants" involved in the resulting macroscopic equations are reasonably well known either from experimental measurement or approximate theoretical considerations. 2.2 The Macroscopic Equations The macroscopic equations governing the average motion of each component in an ionized gas mixture can either be considered approximately as the first three moments of (2.1), or deduced directly from the application of the macroscopic conservation laws of mass, momentum, and energy. These laws, and some explanations concerning them, are given below. Due to the vast amount of literature on this subject, a complete list of references cannot be given in this report and thus only representative references are listed.

(i) The conservation of mass and charge. If the mass density of each component in a mixture is p., and the macroscopic average of its velocity is V., one then has -P+ V. (P.V.) =Q (2.3) where Qj is the material source input. The charge density of each component is given by P.q. a: (2.4) j mj J and the current density is given by I =o V (2.5) Thus, the conservation of mass implies the conservation of charge. (ii) The momentum equation. For each component of a mixture, we may define a pressure tensor to account for the transport of the fluctuation part of momentum. This pressure tensor may be written as P.=p.I -T. (2.6) 3 J-= =J where p. is the partial pressure of the j component, T. is the traceless pressure tensor, and I is the identity tensor. Similarly, the kinetic temperature T. of each component can be defined in terms of the partial pressure by j p. =,L kT. (2.7) m. j J 6

where k is Boltzmann's constant. Due to short range interactions (collisions) between different components, the momentum of each component is not conserved. The usual momentum equation for each component is given by a3 (p [p jVjV + I = I =E+I.xB+F.- k. (V -V.) (2.8) at -- - j -- - -j i —j — i i where F. is the external mechanical force per unit volume acting on the jth -J component, and k..(Vj-V.) is the effective force per unit volume acting on th J.th the j component due to short range interactions (collisions) with the i component. The precise calculation of the coefficients k.. requires a deJ1 tailed knowledge of the mechanism of short range interactions. Approximately, it can be shown that 2mm k (number of binary collisions in unit volume per k Xko c _ — ji ii m.+m. sec. between particles of the two components) i j (2.9) (iii) The energy equation. Due to collisions between particles of different components, the energy of each component in a mixture is not conserved. The transfer of energy due to a collision, again depends on the mechanism of the collision, which, of course, is unknown. Assuming a simple model of collisions, such as an inverse fifth power interaction, it can be shown that the energy equation for each component of the mixture is 7

-V. +U. +V P.V2 +UJ)V.+p.VV.- 'T +. = E I.+F. -V.+h. dt J 2 JJ L J2 j j J -j j-j - j -j I iVi I (V- 3k - v { Ml.+- i). (V -V.)+. M(T -T )l (2.10) ->j 1 \ m+mj / -j -1 m.+m. (T- i (2. 10) where i. is the conventional heat flux, h. is the heat input to the j *th component, and U. is the internal energy per unit mass of the j component. Although the transfer of energy due to collisions given by the last term of (2. 10) was derived from a very simple model, its physical implication is quite clear: collisions between particles of different components tend to equalize the average velocity of the components as well as the random energy (the term involving the temperature difference). Besides these conservation laws, obtained from the mechanical point of view, the long range interaction between particles expressed in terms of the macroscopic electric field and magnetic field must also satisfy Maxwell's field equation. These are V (e E)= a.+a (2.11) J V- B= 0 (2.12) aB VxE - -- = 0 (2.13) - at 8

aE VxB- pE I_+I (2.14) VxB 0 oo jt 1i -( i where p and c are, respectively, the free space permeability and pero 0 mittivity (: = c, the velocity of light), a-, and I are the external s - s S oo sources of charge density and current density. The solution of the equations obtained from the conservation laws together with Maxwell's equations, even with some known relations between T J i, v and the other macroscopic variables, is too complicated to =' jk yield any immediate result. Therefore, some approximations, neglecting secondary effects but retaining the basic features of primary physical importance, are necessary. From the system point of view, the system of equations may be considered as several subsystems coupled together, i. e. the systems for the motion of each component, and the electromagnetic system. The motion of each component is coupled through (i) the collision terms and (ii) the electric and magnetic influence via Maxwell's equations. If this coupling is neglected, each component may be treated as a simple gas. For a disturbance in a simple gas, the conventional theory of acoustic waves applies. In acoustic theory we assume that, to a first order of approximation, the entropy remains constant and the acoustic velocity is given by ap constant entropy. 9

Based on the success of using linear theory and the constant entropy approximation in the investigation of acoustic disturbances, even in gas mixtures, our present investigation shall be concerned primarily with the effect of the coupling terms on the acoustic disturbances. Approximate equations that degenerate to the acoustic equations when coupling is neglected are given in the next section. 2. 3 The Perturbation Equations For the investigation of disturbances in the ionosphere, we shall simplify the macroscopic equations deduced above by using approximations consistent with our present day knowledge of the ionosphere. Therefore, the following idealizations are made in this investigation: a) The ionosphere is assumed to contain three distinct components: electrons, ions (singly charged) and neutrals. The average mass per particle of each component is designated by m, m., and m, respectively. b) When undisturbed, the average temperature and number density are known as a function of altitude. The temperatures of the components are denoted by T, T., and T, respectively, while the number densities are dee 1 n noted by N, N., and N, respectively. In general, local electrical neutrality may be assumed, so that N = N. N. (2.15) e 1 o The collision coefficients, k., k, and k., given by Eq. (2. 9) can thereei en in fore be estimated approximately. c) The disturbances are assumed to be weak, so that linearized equations apply. The perturbation of the electron, ion, and neutral density is denoted by n, n., andn, respectively. e i n 10

d) The thermal and viscous effects are considered to be of secondary importance, so that the stress tensor and heat flux can be neglected, and for each component j, 2 V p.- m.U. Vn. (2.16) th where U. may be regarded as the local acoustic velocity of the j component. e) The perturbed electric field is denoted by E, and the perturbed magnetic field by M h. Due to the presence of the earth's magnetic field, the total magnetic field is given by A BB +oh=B b+p h (2.17) - -O- 0 - A where B is the magnitude, and b is the unit vector indicating the direction of the earth's magnetic field. Based on the above assumptions, the equations in Section 2.2 may be linearized by neglecting the second order terms. Taking the Fourier transform (with respect to time) of the resulting linearized equations, the following equations governing the perturbed motion of electrons, ions, and neutral particles are obtained. Electron Motion 2 U (e e r - iwV + Vn = — E+V x B -e N e m L- -e o e k k F ei en e _ (V -V )- en(V -V )+ Nm -e -i Nm -e -n Nm e o e o e (2.18) 11

-iun +V (NV ) =Q e o-e e (2.19) where e is the electron charge (numerical value). Ion Motion 2 o i 1 e -ioV.+ Vr (V.) -i N 1. 1 0 1 k k ei in (V -V.)-V r (V-V )+ F. Nm. e 1 N m - -n -1 0 1 0 1 -i.Z +V, (NV) V Q 1 o-i Neutral Particle Motion (2.20) (2.21) 2 -iwV + Vn = -n N n n k. _- h (V -V.) N m -n - n n k ne (V -V )+ F N m -n e -n n n (2.22) -iJmn +V- (N V ) Q n nn n Maxwell's Equation Vx E- iup h = -K - O- - Vx h + iw E= eN (V-V )+ I -- o- — i e s (2.23) (2.24) (2.25) 12

where I is a current source, and where, for completeness, a magnetic -s current source K has been introduced. The set of Eqs. (2.18) through (2. 25) represents a coupled system of 18 first order partial differential equations relating the 18 field quantities (3 components each of, V V, E, h, and the three scalar n, n., n ) to the sources. The solution of such a system, especially in an inhomogeneous (stratified) medium such as the ionosphere is extremely complicated. Therefore, our investigation of this set of equations is limited to: a) The investigation of the characteristics of the propagation of the disturbances in a locally homogeneous medium. h) The development of techniques that may enable one to solve such large scale systems numerically or approximately. c) The investigation of some idealized problems. For example, in some regions of the ionosphere, if we neglect collisions and ion motion, a relatively simple set of equations involving only the motion of electrons is obtained. These are given by Vx E -- iwu h = -K (2.26) O- - - Vxh+ i? E+eNV= I (2.27) o o- -s -imnN V +mU Vn+eN E+VxB = F (2.28) O- O -- - N V- V+ V (N V)- i = Q (2.29) 0 Owhere the subscript e has been deleted from all variables since no confusion would result. The excitation of disturbances in an electron plasma, 13

based on certain limiting cases of the above equations, has been discussed by many authors but the most comprehensive investigation has been carried out by Y. K. Wu under this contract and the results presented in an earlier report (Wu, 1965). d) Some idealized problems involving the excitation of disturbances in a bounded plasma. The result of these investigations will be discussed in detail in the remainder of this report. 14

III PROPAGATION CONSTANTS 3.1 The Dispersion Relation In principle, in a homogeneous medium, the set of Eqs. (2.18) through (2.24) may be solved by the usual Fourier transform technique. Basically, for each function of the space variables r, we may introduce the Fourier transform defined by 1 ( -is r f(s) =- -e Af(r)dr (3.1) (2r) The transform of the set of Eqs. (2. 18) through (2.24) yields a set of linear algebraic equations for the transformed field quantities, which can be solved by inverting the matrix formed by the coefficients of the algebraic equations (functions of s). To carry out the inverse transform for the field components, it is then necessary to find the roots of the determinant of this matrix. This relation is known as the dispersion relation. The truth of this statement can be seen from the following discussion. Instead of applying a Fourier transform to Eqs. (2.18) through (2.24), assume that the components of the is* r vectors and the three scalar quantities vary as We - where W is independent of r and represent the magnitude of V, etc. This is the usual asex sumption of plane wave propagation where s is the propagation vector given by s 8s (3.2) and r is a position vector in space. Making this substitution into Eqs. (2.18) 15

through (2.24), assuming no sources, a system of linear homogeneous equations are obtained. In order that a non-trivial solution exist for this set it is necessary that the determinant of the coefficients be zero, This result yields the propagation vector s as a function of frequency and, therefore, is the dispersion relation. It is obvious that the determinant obtained in this way is identical to that obtained from the matrix resulting from the application of the Fourier transform. Explicitly, if each field variable such as V (r) is replaced by V (r)eis the following equations are obtained*. e s x E ^ /p h (3.3) s xh x -wE E-ieN (V.-V ) (3.4) o- o -i -e N n =- s V (3.5) e w — e N n.i S- V (3.6) N n D s- V (3.7) n w — n 2 E A Uen V (1+iV.+iv )-i -- — il V xb+s- +i+ iv V.+iv V (3.8) -e i en m e-e -N e - en m ee N e e o No confusion should result if V is used to denote V (a) etc. -e -e 16

2 eE U. n V. (+iv. v. ) -+ i -+iQV.xb+s --- -+iv V +iv V (3.9) -i le in m. w i-i -N t ie-e in-n 1 o 2 U n V (l+iv.+iv )= s n +iv V +iv.V. (3.10) -n ni ne - N t ne-e ni-i o For simplicity, the normalized collision frequency ratio v which accounts for the loss of energy of particles of m kind due to collisions with those of the n kind given by k mn mn: x density of mth component has been used in Eqs. (3. 3) through (3.10). To obtain the dispersion relation from the above systems, we shall, without loss of generality, choose A s = zs and B = (z cos 0 +y sin 0) B -o 0o where 0 is the angle between the direction of propagation and the d. c. magnetic field. By expanding the vector equations and quantities in (3.1) through (3.10) into components, the dispersion relation may be given in the following determinantal form: 11 A12 13 A21 22 A 23 (3.11) A A A 31 32 A33 17

where the sub-determinants (minors) Aij are given by the following equations. (2 -2)( +ive +iv ) - 1 ( 2-s2)s sn 2 sZ,-I I I i (32- s2) [A11]= - 2(2O -s in2a9) o I -iQ (AB2) I - -.3 - - -, - ( -s2)82sin os8 I e (0 -e+liv+vn) -e e O, -81 (o-BiB 8^sin0Co I e I e I +w2s2s nhcosO et en I O I (o-s2)(+1 + ve+ en) -Is2)s c2 2 I 1 (2 2 ) e i-5( -8) coe ) 0 22 + w B Bsin coe e I (3.12) 18

[A121.= 2 2 I 0 el 0 1 2si2192 e.e 0 -~- -- - - — IF- -- e e I -o as2ine9 coe9 e I0 I-s I e I (3. 1:3) [A 13} 0- 0 0 1 0 I 0 I "j-s) en FA21]z.43.2) vie I 0 22 2 2 +-iossi8csin I0 II (:3. 14) 28a 2Cose0 sine0 I I0 I 0 l I ~ W 2(fl-S2 Cos82) l o (3-. ~5) 19)

22_ 0 S '+ve+uim -.-0(f3 W82)Ssisn 6 i I I I I I I I I I (f3-s02)sisn60 coos6 2 2 -B ) 1 2 2 1 +(Oia sin 0coos6 2 22 2 — (Aj39-s siOm )I 1 0- - - - - - - - ~ [A22 21 +1Q( -2 1(2-s2)s sine CosO0 2o0 2 *2 +cw)s sin9O cos6 1 f I I I I I 0 le "in 0 I I I I - I - I I I I I (/3 - )(+li,+iv J ) 0 ie in 1 2 22 2 - p-s )s COB6'81 A2 022 -'(4)2(p-.2Cos26) lo 0 (3. 16) [A:1 2 2 -(/3 -s )vin 0 In 0 0 -(2-2) o in -1- - 0 7 0 0 (3. 17) 20o

l+iv+ / ' I n ne 1 2 n I I -.- I - I. - - - - o 0 I 0 2 I i2 2 ne I 113200i20 (3.18) '~ iv1 0 0' I | - e (3.19) -Al ni = 0 0 0 I ~ - - - I ' - o 0 -iv o I I (3.20) 21

For simplicity we have also introduced the following normalized quantities in Eqs. (4. 12) through (4. 20). 6p ~ (free space propagation constant) 0 o c e, i, n V. (propagation constant of respective e, 1, n acoustic waves) e2N e, i we mei (normalized plasma frequency) e, 1 WE~ M o e, i eB 0 - (normalized gyro-frequency) e,i wm e,i 2 Equation (3. 11), when developed, is a polynomial of 5th order in s, indicating the existence of five basic waves (modes). These modes may be interpreted as the result of collisional and other forms of coupling between the following basic types of waves: neutral acoustic, electron acoustic, ion acoustic, ordinary and extraordinary electromagnetic waves. In order to learn more about the effect of "coupling" on these five basic waves, the characteristics of wave propagation in the subsystem will be studied first. This, of course, implies the investigation of the roots of the subdeterminants of Eq. (3 11). 3. 2 The Basic Modes of Propagation Basically, there are five natural modes of propagation in the model of the ionosphere considered, corresponding to the five basic energy storage mechanisms. The dispersion relation, therefore, yields a polynomial of 5th degree in s 2 corresponding to the five basic tes of waves (or modes of degree in s, corresponding to the five basic types of waves (or modes of 22

propagation). Due to the collisional coupling, and the collective action of charged components through the electrostatic and magnetic interaction, identification of each particular mode of propagation is somewhat ambiguous. In order to understand more clearly the roots of the dispersion relation, some general characteristics of the basic modes of propagation are discussed in this section. Since the thermal energy exchange was neglected in the present first order analysis, the effect of collisions would be to cause attenuation of the waves due to the transfer of ordered energy into random energy. When the appropriate collision terms are small, v << 1, the effect of collisions on the phase velocity of propagation is small and the nature of the basic waves (or modes) of propagation can be determined by neglecting collisions. The effect of collisions can then be estimated by the following consideration. A particle of the ith kind, with ordered motion, will lose some of its energy (ordered motion) when it collides with a different type of particle. This energy loss represents the maximum attenuation due to collisions since some of the ordered energy transferred to the second particle will eventually be returned (fedback) to a particle of the ith kind by means of a series of collisions with other particles. Therefore, the energy lost in the initial collision can be considered as the attenuation due to collisions and the feedback of energy as the coupling effect of collisions. As a result of the above discussion, the nature of the waves (modes) will be studied by first neglecting collisions and then the attenuation estimated by considering only the attenuation due to collisions. a. Neutral Acoustic Waves If collisions are neglected, the motion of neutrals can be considered separately. The equations describing the motion of neutrals are: 23

N n = s- V (3.7) n l - -n 2 U n n n ss V = s ~V (3.10) -n -N w 2 -n o / n These equations yield the basic characteristics of the acoustic type waves, namely they are longitudinal waves in the sense that the velocity is in the direction of propagation. The propagation constant of the acoustic wave is given by s =: (3.21) n U n and is consistent with the definition of acoustic velocity. Collisional coupling would modify the above conclusion somewhat. To estimate the attenuation of the acoustic wave, we may modify (3.21) to s = 0 J i.+W ] n + ni ne k.+k U l+i n- (3.22) U wN m n n n Equation (3. 22) gives the limiting value of the attenuation constant for acoustic waves at high frequency which is k +k k1 ni ne. (3 23) 2U n n n n n 24

b. Electron Plasma (Acoustic) Waves The acoustic type of wave for electron motion is not as simple as that for neutrals. Since the motion of charged particles necessarily causes current and bunching of charges, the mechanical equation of motion of the electrons should be considered together with Maxwell's equations. If we neglect collisions and the d. c. magnetic field, and assume that the motion of the ions is negligible due to their heavy mass, the coupled equations to be considered are sx E= h (3. 3) s x h = -We E+ieN V (3.4') - - - oe N 0 n P-s V (3. 5) e - -e e E s n V =-i — - (3.71) -e m w N i e o The above equations may be used to eliminate h and n, yielding I1- + E - 0V (3.24) L 2 E E o e 0 0 and [1 — V = - m E (3.25) 2J-e Me e 25

The nature of the wave may be made clear by separating E, and V into components in the direction of propagation, (the longitudinal come ponents, EL Ve) and the components transverse to the direction of L' eL propagation, ET and V. In terms of these components, we have T. eT ieN E 0- V (3.26) L w c eL o 2 ieN (1- ) E 0 V (3. 27) 2 T E eT 0 O 2 (1-s ) V - E (3.28) 2 eL m L P e e V i - ie E (3.29) eT m W T e From Eqs. (3. 26) and (3. 28) we find that the longitudinal components of V and E are associated with the propagation constant, -e s = 1-2 (3.30) e e which has the characteristics of an acoustic wave, and is generally called the electron acoustic (plasma) wave. On the other hand, from Eq. (3. 27) and (3. 29), we see that the transverse components are associated with the propagation constant 26

2 s = 1-c (3.31) o e This transverse wave can be identified as the electromagnetic wave. One may also estimate the attenuation of these waves by considering the attenuation effect of collisions and neglecting the coupling effect. For the electron acoustic wave, we have s ~2 s i e +iv +iv. (3.32) e e ei en For the electromagnetic wave, we have, 2 2 e s 1- (3.33) o 1+iv +1iY ei en c. Ion Plasma (acoustic) Wave When the ion motion is not negligible, such as in the case of a low frequency, collionless plasma, the coupled motion of electrons and ions is described by the following equations. 2 BS ieN 1- + -- E -- (V -V ) (3.34) r s L 2 2 - m c ie e3 e e 27

r s s. 1 V-2-. =+ i- E (3.36) 2 B i minu i 1 Again, the longitudinal components are associated with acoustic types of waves. The propagation constants for these waves are determined by the equation s - (1-W )+ (I(-o)I s +3 f(1- -t. ) = 0 (3.37) Le e i i e i e 1 The roots of this equation are given by 2 1 2 (1 2 2 e.2 )- 2 (1- w 2 2 2 22 2 -e e i i Le e i e i e (3.38) Equation (3. 38) shows that the two waves may be considered as the result of the coupling, by electric forces, of an ion acoustic wave and an electron acoustic wave and, in general, cannot be separated into an ion and electron acoustic 2 2 wave. However, at high frequencies where the coupling is weak (W <w <1) 1 e the roots are approximately 1 2l smI. l-U2 (3. 39) and 8sd3 (3.40) e e where (3. 39) is determined by the positive sign in (3. 38) and may be identified as the ion acoustic wave and (3.40) is given by the negative sign and may be identified as the electron acoustic wave. 28

The effect of collisions at high frequencies can be estimated by neglecting the coupling effect of collision only and (3. 39) and (3.40) become s I-=. 2+iV. +iV. (3.41) 1 1 ie n and sxa 1-W +iV e+iv (3.42) e e e ei en At lower frequencies where the coupling effects are large the two waves, of course, cannot be identified as an ion or electron wave; however, since.i >> 3 we can say that an ion type of acoustic wave is given by the positive sign in Eq. (3. 38) and an electron type of acoustic wave by the negative sign. As the coupling effects become larger, the motion of the electrons and ions become nearly equal and only a single wave propagates. This point can be determined from (3. 37) and is 2 2 W 2+(. 1 e 1 22 2 2 When we +W.2 < 1 two waves propagate while for w +W. > 1 only e i e 1 one wave can propagate. Similarly, the transverse component of the velocity and electric field is associated with the electromagnetic wave. This electromagnetic wave has a propagation constant given by 1- 2Jo- - W. (3.43) 29

d. Electromagnetic Waves The transverse characteristics of electromagnetic waves are modified due to the presence of the d. c. magnetic field. Inspection of Eq, (3. 3) reveals that the perturbed magnetic field is always transverse to the direction of propagation, while the electric field, in general, has a component in the direction of propagation. The presence of a d. c. magnetic field causes the electromagnetic field to split into two distinct modes. Based on the propagation characteristics of these waves, they are referred to as ordinary and extraordinary waves. However, in the analysis of a warm plasma problem, these waves are coupled together and, in general, lose their distinguishing characteristics. An appropriate terminology for these waves has been given by Allis et al (1963) and discussed in detail by Wu (1965). This terminology is based on the characteristics of propagation of these waves when the direction of propagation is at a fixed angle with respect to the d. c. magnetic field. In this section, the well-known characteristics of these waves will be deduced * in order to clarify the terminology used to describe them. The basic characteristics of electromagnetic waves in a plasma under the influence of d. c. magnetic field can be brought out by considering the waves in a cold, collisionless plasma. The equations to be solved are therefore, 2 ss ieN ieN (1U- + -)E —V - V. (3.44) 2 2 w -e we -i 0 0 ie E A V -._e -inl V xb (3.45) -e m e -e e This deduction follows the procedure of Stix (1962). 30

and ie E V. + -+i.V.xb (3.46) -1 m o 1ie Since the d. c. magnetic field tends to force the charged particles to have a circular motion in a plane perpendicular to the direction of the d. c. magnetic field, it is convenient to refer all directed quantities to the direction of the d. c. magnetic field. Therefore, without loss of generality, let b x z b=z and A s X s(z cos -ysin ). In terms of the components of E, V and V. Eqs. (3.44), (3.45) and (3.46) — e -1 become 2 ieN ieN (1-s-)E = — V. +. V (3.47) 2 x Oe ix (E ex o o o 2 2 ieN ieN s 2 o o (1- sin 8)E - sin cos E z — V. +- V (3.48) 2 z 2 y w 1z ez o o 2 2 ieN ieN s s 2 o o. ---sincoeO E +(1 Cos )'E -- V. + v (3.49) 2 z 2 y WE LY W ey o o0 31

ie V - E (3.50) ez mo z e E ie x V ---- -i V (3.51) ex m w e ey e E V = -ie-. L+iR V (3.52) ey m W e ex e The equations for the ion velocities are similar to those given in Eqs. (3. 50) to (3.52), and may be obtained from them by replacing e and Q by -e and e., respectively. 1 Although the dispersion relations can be obtained directly by forming the determinants of the above equations, the basic characteristics of the waves can be made clear if we consider the component equations separately. Eliminating the z-component of velocity yields a relation between E and E: z y 2 2 2 (P- sin )E n sin0 cosO E (3.53) z y where, for simplicity, we let A 2 2 P- (1-) -.2 (3.54) e i and 2 2 s n = (3.55) 2 32

Similarly, the transverse components of velocity may be expressed in terms of the transverse components of the electric fields by use of Eqs. (3.51) and (3.52). This is accomplished in simple manner by forming the right hand rotating and left hand rotating field components such as E + iE, E - iE, etc. x y x y Using the rotating field components in (3. 51) and (3. 52) we obtain ie (1- )(V +iV )= - (E +iE ) (3.56) e ex ey m ew x y e (1- )(V -iV ) ~ F(E -iE ) (3.57) e ex ey m eW x y and similar equations for the ions. By substituting these results into (3.47) and (3.40) the following equations are obtained. (2 n )E iE () (3.58) 2 x y and 2 R+L 2 2.R-L -n cossinEos E)E +,i 2 E (3.59) z 2 y 2 x where 2 2 e 1. 1-2 140. e i 2 2 140 1-0f e i From the determinant formed from the coefficients of Eqs. (3. 53), (3. 58) 33

and (3. 59) the equation for s, or n is given by 4 2 R+L 2 2 (R+L) 24Ti 2 n (Pcos 6+- sin 0)-n [p.-(Hcos t)+RLsin 9+RPLx 0 2 2 (3. 60) Equation (3. 60) indicates the anisotropic behavior of the propagation (direction dependence of propagation constant) of the electromagnetic waves, where as the acoustic type of waves discussed previously are isotropic. The directional characteristics, and the identification of these waves can be simply obtained by expressing Eq. (3. 60) in the form 2 2 2 P(n -R)(n -L) tan R+L 2 2RL 2 (n -- )(n -P) 2 R+L (3. 61) Thus, in the direction 0 0, the two waves have the following features. 2 a) n *R, 2 2 we I sxf 1- le - e i (3.62) and from Eq. (3. 58) E * -iE x y which is the characteristic of a right hand circularly polarized wave, and b) n 2L, 2 2 e 1 (3.63) which corresponds to a left hand circularly polarized wave. 34

In the direction 0 = 90, (propagation transverse to the direction of the magnetic field) we have 2 a) n -P s3= \l -u 2 (3.64) o e 1 which is a wave uneffected by the magnetic field and, therefore, is generally referred to as the ordinary wave, and 2 2RL b) nR+ L 2 2 2 2 W U. W W. e 1 e 1 (1- - 1 )(1 - ) s= e 1 e(3.65) 0 2 2 ( e 1 (1\ - e 1 which is generally referred to as the extraordinary wave. In general, for 0 < 0 < e the two waves are coupled and lose their 2 distinguishing characteristics. However, they can be distinguished mathematically by the limiting condition as 0 -- 90. Waves whose propagation constants are obtained from the same branch of the solution of Eq. (3. 61) as the ordinary wave, (3. 64) will be called the modified ordinary waves and the waves whose propagation constant is given by the same branch as the extraordinary wave, (3. 65), will be called modified extraordinary waves. 35

The dispersion relation, given in the form of Eq. (3. 61) can also be utilized to identify the regions where each of the waves can exist along a 2 particular direction (or are heavily attenuated). Since tan 0 is positive 2 along each mathematical branch of the solution, n is a continuous function of 0, and waves are heavily attenuated if n < 0. One may, according to the values of P, R, and L, investigate the possibility of the existence of unattenuated propagation by using the well-known CMA diagram. For the case of a cold plasma, such a diagram is given by Stix (1962). The application of such a diagram to the case of warm plasma was given by Wu (1965). e. Waves in a Warm Plasma The analysis of the dispersion relation pertaining to the propagation of waves in a warm plasma, where the effects of thermal velocities are considered, is generally quite complicated, even for the collisionless case. For an electron plasma, the determinantal equation was laboriously developed and discussed by Wu (1965). A simplified approach to this problem, based on the discussion given in the previous section, is to investigate the proportionality 0 2 constant only in the directions = 0 and = 90. If the values of n in these directions are known, say a1 a2 and 3 for 0 0 and 31, 2 and 3 for 0 X 90, the dispersion relation may be easily written in the form 2 2 2 (n -a 1X -a2 Xn -a) tan O constant2 2 2 (3.66) (n - )(n -3_2Xn -3) 1 2 3 36

where the constant term can be obtained easily from the determinantal equation by a limiting process. For example, consider the case of a collisionless electron plasma. The set of Maxwell's are given by Eqs. (3. 47), (3. 48) and (3. 49) while the equations of motion for the electrons are obtained by adding an acoustic velocity term to Eqs. (3. 50), (3. 51) and (3. 52). The complete set of equations are ieN (1-n )E - V (3.67) x WE ex o 2 ieN 2 2 2 o (1-n sin O)E -n sin 0cosOE V (3.68) z y wE ez O ieN 2 2 2 o -n sin cos E +(1-n cos0 )E = V (3.69) z y we ey 2 E 2 a 2 ie z s (1 — cos 0)V = - - in cos V (3.70) 2 ez m o 2 ey B e e o E V -ie- -iV (3.71) ex m w e ey e 2. E 2 ( — sin)V = — - 8e sinOcos V (3.72) 2 ey m 2 ez Pe e e 37

For 0=0 we have from (3. 68) 2 U. E 2 e ie z (1-n -)V = -v 2 ez m w A e (3.73) While the transverse components of the velocity are identical to those obtained for the case of a cold plasma and are given by Eqs. (3. 64) and (3. 65). The z-component of the system of equations is decoupled and yields one propagation constant as s= ' J - 2' e e (3. 74) which is easily identified as the electron plasma wave. The same procedure used to eliminate V and V for the cold plasma can be applied to the reex ey maining equations yielding the propagation constants,Szf 4Rii'2 1- e 1-e e (3. 75) and s=f3 1L a 0 2 \1 e (3.76) which are easily identified as electromagnetic waves. 38

For 0 90 the components of the equation again decouple yielding the propagation constant I 2 s o\ 1-e (3.31) which is identified as the ordinary electromagnetic wave. For the transverse field, E and E can be eliminated yielding x y 2 W e V (1- )+ i2 V 0 (3.77) ex 2 e ey 1 -o 0 -iF V + - - V =. (3.78) e ex e 2 ey e A second order equation in s can be obtained from (3. 77) and (3. 78) and is S4-2S [ 2 2Xi- 2 2 2)] +2 2 2 [i-w )2- 2] 0 e 0o e e e 0( e - e e 0 (3.79) The roots of (3. 79) can be investigated easily and correspond to waves obtained due to the coupling between the extraordinary electromagnetic wave and the electron acoustic wave. 39

The dispersion relation can then be obtained from the propagation constants and Eq. (3. 66) and is 2 2 e 2 ^0 tan O - 2 (3.80) [n2 -][n4 2 2 e eR+L 2 d) -P -n (1-n -RL O The constant in (3. 66) was evaluated by comparing the limit of (3. 66) as 3 -- ao to the cold plasma equation. A detailed discussion of Eq. (3.80) was e given by Wu (1965) in a previous report. The same approach, when applied to a warm plasma, considering both electron and ion motion, does not yield such a simple result. For the sake of completeness, the result is given below. 2 2 2 2 2 2 2 2 2 2 2 2 2 tan(2 -B12Ls2 -2R)[4- [i2(1 -w)+f (l- 2 s2+P.(1- - 2tan 2 6 3+ a[(s-2 P) s +A +A2s +A3 (3.81) where F 2 2 2 2 2 2 2 2 A, -.(1-c 2-/ )+P (1-o )+- (1-u -cg 1 i1=1 e e e 0 e 1 40

922 r 2 2 9 2 2 29q 2 2 A2=3 [-x -La -Q (1-u)-Q (1- )+Ql] 2 e e 1 1 e e 1 2 2 2.2 2 2 2 2132 2r22 2 _ 2(12 2 + P2. (1-w -u )(-2)-l 2- ) 1- )(1-w - (1 -O i e I I 1 e J e e 1 e e 1 A 2- 2 2 2 22 22_22 22 2)2 2Q22 2 L. [1-W - + ) +2S.2 w. -(1-w - (1- WZ 2 2 3 o e e 1 e 1 e e 1 1 e e i 2 For the sake of comparison, Eq. (3. 81) has been written in terms of s 2 rather than n. It is easy to see, from Eq. (3. 81), that for = 0 two coupled electromagnetic waves are possible and two acoustic waves, due to the coupling o between the electron and ion acoustic waves. For 0 = 90 the first term in the denominator is precisely the same as Eq. (3. 64) and is the ordinary electromagnetic wave. The second term in the denominator yields three waves, an electromagnetic wave and two acoustic waves, which are the result of coupling between the extraordinary electromagnetic wave and the electron and ion acoustic waves. 3. 3 Propagation Constants in the Ionosphere The characteristic waves for the collisionless case discussed in the previous section yield the basic characteristics of propagation, in general, only for the direction of propagation parallel or perpendicular to the direction of the d. c. magnetic field. For intermediate directions of propagation, the basic modes are coupled and lose their identity. Indeed, even for the parallel and perpendicular case, the basic modes are coupled as indicated, for example, 41

by Eq. (3. 38). The effect of collisions, as mentioned in Section 3.1, can probably be considered as composed of two parts, attenuation due to collisions and coupling due to collisions. Examination of the curves in AppendixA shows that for the ionosphere, especially for charged particles, the coupling due to collisions is probably of secondary importance when compared to coupling due to electric forces and due to the d. c. magnetic field as exemplified by the plasma frequencies and gyro frequencies, respectively. Because of this, it is probable that the most important effect of collisions is the attenuation of the waves. The attenuation of the waves can be considered as due to two causes. The coupling due to electric and d. c. magnetic forces and that due to collisions. At frequencies and regions of the ionosphere where the coupling forces are small, the collisions terms (also small) have the greatest effect on the attenuation of the wave and have very little effect on the phase velocity of the waves. As the coupling forces become large, the attenuation due to collisions becomes small compared to electric and magnetic forces and in this region, the collisions terms have their greatest effect on the phase velocity. This point can best be illustrated by examination of Eq. (3. 32). The real and imaginary parts of this equation are 2 2)2+ 2 )2 1-a + (1-u ) +(v.+v 2 e e el en r e 2 -(l-w)2+ l-w 2)2+( i+V )2 2 e, ei en i e 2 42

Since for the regions of the ionosphere considered v.+v <<a it is 2 el en e evident that for W <1 s is a strong function of the plasma frequency and e r 2 s. is largely a function of the collision terms. For o > 1, due to the 1 e factor 1 — w the reverse situation is true. e It should be noted that the above analysis does not apply to the neutral wave since, of course, electric and magnetic forces do not affect the neutral particles directly. Because of this, the propagation constant for the neutral wave should be given to a good approximation, by Eq. (3.2) for the ionosphere. One other fact of significance should be pointed out at this time in regard to the attenuation of the electron and ion types of waves. For loose coupling, 2 << 1, Eqs. (3. 32) and (3.41) are approximately, e (v.+ ) S e+i el en e 2U e (v. +v.) SzA3 +i ie m i1 2U. 1 Thus, the attenuation of the electron and ion acoustic waves is given by the ratio of the collision frequency to the appropriate velocity. Again, referring to Appendix A, this number is relatively large and these waves are rapidly attenuated, in general, and the - point for the magnitude of these waves e is reached in a distance of one kilometer at the most and,in many cases, in a distance of a few meters. A complete analysis of the effects of collision and direction of propagation on the modes propagating in the ionosphere would involve a detailed analysis of the dispersion relation, Eq. (3.11). Because of the complicated nature of this 43

equation and the analytic expressions obtained for the coefficients of the resulting fifth order polynomial in s, the dispersion relation was programmed for analysis by a digital computer. A complete discussion of the program is given in Appendix B. Due to limitations of time and money, only a limited amount of data was obtained and the results tabulated in Table I. The propagation constants were computed for 15 intervals of 0 between 0 and 90 at angular frequencies of 3x10, 3x 10, 3x1O, 3x1O, 3x1O, 3x1O8 and 3x1O radians per second, using parameters appropriate to an altitude of 100 Km in the ionosphere. The roots are presented in groups of five, corresponding to the five basic modes of propagation, the first two roots correspond to the two electromagnetic waves and the last three to the electron, ion and neutral acoustic wave, respectively. The columns in the table correspond to variations in direction of propagations and the rows, in groups of five, to frequency variation, the lowest frequency appearing at the top of the table and the highest frequency at the bottom. The left hand column of the table lists, in addition to frequency, the corresponding values of f30 3e i and 3n while the two extreme right hand columns list the appropriate normalized values of the plasma frequencies, gyro-frequencies and the collision frequencies. Only three significant digits are listed for the real and imaginary parts of the propagation constants; however, many more digits were required to achieve the required accuracy. In this table, the imaginary parts of some of the 4ots are negative. One particular case of importance is that for the neutral wave 3 at 90 and wx 3x10. The negative sign was due to a lack of accuracy in the procedure for determining the roots of the polynomial in that not enough iterations were carried out in the process. The roots of this polynomial were evaluated a second time, using greater accuracy, and the following results obtained for the five roots. 44

TABLE I; VARIATION OF PROPAGATION CONSTANTS ALT 100 Km. WITH ANGLE 9 AND FREQUENCY 4 u..3x 10 30.307010 3 -. 103x 102I 0w=.7o3x1I01 10~31x 1002 3z 10x 10 * 3x1064 l0 3 =.31 x100 7.x 102 0 =3x103 3~0 z~.31x102 3. 7x10 3 o.3x10 710 3: Ax 106 3 =. 10'2 j3=.31x 10 7x 104 w=.3 x 10 -2 =.310'0 ~ =31xl103 e 13.7 x10 5 t=.1x1086 00.139xI0 s-4+i.612 x103. 592 x10 -3+j. 132 x10 -.702 x -IO +i.141 x103.533x10 1+i. 188x 10.103x 10 2j.171 x10 - -4 -2.173xl0 +i.190X10 -2 -4 *190x 10 +i. 175 x10. 566x 10 0+i. 134 x103.528x 10 2+i.234 x10I.103x 10 3+i. OOxlO.519 x10 4+ i. 583 x102. 620x 10- + i.559 x10 -.325 x1 1 +j.l101xIO0 3 *1I.697xl0 +i.286x10.103 x10 4+i. 00x10. 137 x10 -3+i. 13 x10 -.255 xI0- 1+i. 292xl103.26 x100 +i.953 x102.704 x10 4+i. 283x 10.103 x10 5+ i. O~xlO.922 x10- +ji.604 x10 -.958 x10- +i. 174 x10 -.290x 10 3+i.852 x100.704 x10 5+ i.283 x I0.103x 10 6+i. O~xlO 10' +i.295X 106 10 0+i. 263 x106.307 x10 4+i. 806 x100.704 x10 6+i. 283 xI0 lxIO 1+i. 409 x108 1x I 1.19x1-8 l30xIO 5+i.806x 100 7ox1 +i o 10 15 0 *143x 10 -4+i. 622 x103. 602 xIO' + i. 137 x10 -.702 x10 -1+i. 141 x103. 533 x10 I+i. 188 x10 *103 x10 2j.171 x106 *182 x 10 -4 + i. 193 x 10 -2 *194 x 10- + i.184 x10 -.567x 10 0+i. 134 x103 *528x 10 2+i.234 x10 * 103 x 10 + i. 00 x 10 *544 x10 4+ i. 593 x10 -*631 x10 4-2+i.591 x10 - *343x 10- 1+j.l10lxlO0.697x 10 3+i.286x 10 *103 xl0 4+i. 00x 10.141 x10- 3+i.132 x10 -* 262 x10- 1+i.327 x103 *289x10 0+i. 949 xlI02.704 x10 4+i. 283 x10 103 x10 5+ i. 0x10 *923 x10- L+i.595 x10 -.958x 10 + i.178 x10 - *29 x10 3+i. 862 x10 *704 x10 5+i. 283 x10.103 x10 6+i. O~xIO xlI 0+ i. 295 xl0b l xlO 0+i. 263xl106 *307x 10 4+i.806x 10 *704 x10 6+i.283x 10 * 103 x10 7+i. O~xlO xlI 0- i.124x 10 -x 1 01i.14x1-5 *307x 10 5+i. 806x 10 * 704 x 10 7+i. 283 x 10 *103x 10 8+i. OOxlO 300 *636x 10 3+i. 152x1I0-.159x 10 -4+i.657x1I0 *703x -10+i. 141 x103 * 533x 10 I+i.188x 10I.103x 10 2-i. 171 x106 -4 -2.214x10 +i.204xIO. 204 xl02 +i.217xl10-.566x 10 0+i. 134 x103.528 x10 2+i. 235 x10I *103 x 10 3+i. 00x10.634 x10 4+ i. 625 x10-.668 x102 + i.705 x10 -.407 x10 I+i.10lx IO.697x 10 3+i. 286 x10I *103x 10 4+ i. OOxlO.153 x10O 3+i.138 x -10.290 xI0- I+ i.473 x103.408 x 10 + i.933 x102.704 x10 4+ i.283x 10I.103 x 1 0 + i. O~xl.925 x1I +ji. 567x 10 -.956xl0O + i.191 x10 -.287x 10 3+ i.886x 10.704 x10 5+i.283x 10I. 103 x10 6ti. O~xlO l xlO 0+ i.293 x106 l xlO 0+ i.265 x106.307 x 10 4+i.807 x10 *704 x1I06 ~i.283 x10I.103 x 10 7+i. OOxlO 1 x10,+ i. 109 x10-6 l xlO - i.104 x106.307 x 10 5+i.806 x10 7 1.704xlO +i.283x10.103xl0 8+j.O0XlO 450.703x 10 -3+i. 188x10.198x 10 -4+i. 728x10-.704 x I0 1+i.141x10. 533xl0 I+i.188xl0 *lO3xlO 2_i.171 x106.289 xl0 4+ i. 226x 102 *226x10 -2+i. 293 x10 -.567 x100 +i.134 x103 *528x10 2+i. 235x 10 *103 xl0 3+i. O~xlO *839 xI0 4+ i. 688 x102 *743xl102 + i. 980x 10 - *572 x10- 1+i.l10lxIO0 *697x,10 3+j. 286 x10 *IO3xlO 4+.0Ox lo *175 x10- + i.149x -10.364x10- + i.104 x102. 863 xI 0 0+i.886 x10.704 x 0 4+i. 283 x10.103x 10 5+i. O~xIO *928 x10- +ji.525 x10 -.954 x10- + i.213 x10 - * 283 x 10 3+ i. 917 x 10.704 xl0 5+i.283x 10.103 x10 6+ i. O~xlO l xlI 0+ i.290x l0 -xlI 0+ i.267xl106 *307x 10 4+i.807 x10.704 xI0 6+i.283 x10 *lO03xl 17+i.-OOXlO0 xlI 0+ i. 146 x106 xlI 0- i.140x 106 *307x10 5+i. 806x 10 713 08 O 10 600 *832 x10 -3+. 285 x10 *307 x10 -4+i. 867 x103 *706 x10- I+l.141 x103 *533 x10 I+ i.189 x10I 10 02 _.11x 0-6.484 x10 4+ i. 268 x102 * 269 x 10- l. 492 x10 *569x 10 +i-. 134 x103.528 x10 2+i.235xl10 *103 x10 3+ji. O~xIO0.135 x10- +i. 809 x102.896 x102 + i.173 x103.107 x10 0+j.l10lxlO0.697 xl0 3+i. 286 xIO.103x1lo 4+ i. O~xIO0.217x 10- 3+i.165 x 1 0.771 x l0 I+i.118x -10. 640x 10 I+ji.625x 10.704 x10 4+ i.283xl10.103 xl0 5+i. OOxIO.932 xl10 + i. 474 x10 -.951 x10 I+ i.242xl10 -.2 8x 10 3+i. 947 x10.704 x10 5+j.283 xl0. 103x 10 0+j~O~xIO0 l xlIO + i.287 x10 -1xI O + i.271 x 10-.307 x10 4+ji.807 x10.704 x10 6+i.283 xl0. 103 x10 7+ i. OOxIO0 l x 1 0 + i. 982 x I10 l xl1o - i.926xl10.307x 10 5+i. 806 x 0.704 xlO 7+i.283 x10. 103 xl0 8+ji. OOxlO 750.785x 10 -4+i. 122 x102.114 x102 +i. 686xl10 * 534 x 10 I+ i. 191 x 10.715 x10- +i. 141 x103.103 x10 2_i. 171 x106 -3 -2 *128x10 +i.372x10 -2 -3.375xl0 +i.132xl0. 582 x 10 0+ i. 134 x 103 *528x 10 2+ji.237 x10. 103 x10 3+i. O~xIO.319x 10 -3+i. 108XI0 --1 -3.131x1o +i.545xl0 *393 x10 0+ji.lIxlIO.697xl0 3+i.286x10.103xl0 4+j.ooxIO.218 x10Z + i.426 x -10.297 x103 + i.19lx lo-. 104 x10 +ji.l162 xlO0 * 704 x 10 4+i.283x 10I.103x 10 5+j. O~xlO * -4.937xl0 +i.425x10.948xl 10+i.276xlo0 -277x10 3+i.967xlO *704 x 1 5+ji.283 x10 l13x 10 0+ i.283x 1060 0 -6 lX10 +i.283xl0 Ix100i 24 x 06.307 xl0 4+ji.807xlO0.704 x10 6+ i.283 x10. 103 x10 +ji.O00YlO0 1x 1 0 +i.19x1-6 1- -6 lxO-.104 x10.307x 10 5+ji.806 x10.704 x 1 7+ji.283 x10.103 x10 8+ji. O~xlO 900.315.- 10 2+ i.321 x102 -2 -2 *227 x10 +j. 115 x10.346 x10 0+i.141 x10.785 x10 I+i.619 x10.103 x10 2-i. 104 x10 -.899 x 10- + i.109 x1 0.976 x 10- +ji.lOI 1xlO *405 x 10 1+ i. 134 x 103.522 x10 2+ji.772 x10.103 x10 3+ji. O~xlO0.733 x10- + i.298 xI 10.338 x102 + i.236 x1I0.129 xlo 2+ i.129 xl03.697 x10 34i.286 xI0.103 x10 4+ji. O~xIO. 896xI103 + i.3lx lO-.395xl103 +ji.214 x -10.181 xl0 I+ji.131 x10.704xl10 4+ i.283x 10 *103xl10 5+ji. O~xIO.939xl10 -1+ji.402 xIO.945 x10 I+ji.294 x 1O.276x 10 3+ji.974 x10 *704 xlO 5+ i.283xlto *103x 10 0+ji. O~xl 1x10 0+i.49x1-6 1xI0 0+i.18x1-6 *307 x10 +i. 807 xl10.704 x10 6+ i.283 xl0 I xlO -ji. 123 x10 l x 1 0 + i. 124 x 10 -.307 x10 5+i. 806 x 0. 704 xl0 7+ji.283 x10 *103 xl0 8+ i.OOXlO0 2 8..20x 10' Q2. 29 x 104 Q=. 53 xl0 kr 101 L) li11x 106 2.2 x 101 SI. 29 x 10 ~.53 x 102 k -102.2 x1-1 Il29 x 102 SI-. 53 x 103 k10 k =.l 103 2 2 -u.xI 129 x 10I RG.53 x10 - k =10 - 2 0 L)2 x-5 l29 x100 11..53 x105 w x1 Il- 29x1 2) -4 1x1 2 -9xI -2 k 10 2 k. 52x10 kT. 19 x105 k. -.80x 10 k =.29x108o k..23xl108 k -.52x10 -6 k.=. 19x10 i1 k.=. 8xl0 k =.28xl109 k ~23 x10-9 k =.52x1o00 -7 k.=.-19 x10 -2 k.=. 8 x102 -10 k =.29xl 10 k =.23xlI 10 k =.52x10 k.. 19 x10 --3 k..-8 x 103 k =.29xl0 -11 k =. 23 x10 k =.52xl102 -9 k.=. 19x 10O -4 k.=.BxlO -12 k-..29xl0 -12 k =.23xIO -J k -.52x IO k =.19xlO 0 -5 k.=8xl0 in k -.29xlO 1 -13 k -. 23 x 0 k =.52xl0 - k. 19xl0 -6 k =.8x 10 -14 k =.29 xlO -14 k.23xlO

.139x10 +i.612xlO 3.592x103 + i. 132xlO4.702x101 +i.141x10.533xl01 +i. 188x101.103xlO0 +i.251xlO7 It should be noted that all numbers agree except for the imaginary part of the root corresponding to the neutral wave which is now positive, as it should be. Several other examples of roots with negative imaginary parts occur 9 for w= 3xlO 0 for one of the electromagnetic waves. It is obvious that in this case the negative part is due to the numerical calculations and is not a representative number. Thus, it would seem reasonable that the imaginary parts of the roots are probably not reliable when they are much smaller than the corresponding real part. In spite of the discussion, it is felt that the following conclusions are justified, based on Table I and the previous discussion. 1. For very high frequencies, the real part of the propagation constant for all modes is essentially given by 3, 0, P. and j3, respectively. O e 1 n 2. The electron and ion types of acoustic waves are highly attenuated, i.e. the magnitude of this type of wave decreases to 1/e of its initial value, in traveling only a few meters. 3. The attenuation of the electromagnetic waves is dependent on the collision frequencies for ( >~ i 2 but becomes large for ^< 2. e - e 46

4. The attenuation of the neutral acoustic wave is relatively low and it will propagate over relatively long distances without appreciable loss(hundreds of kilometers). 47

IV THE OPERATOR TRANSFORM METHOD 4.1 The Operator Approach The set of linearized equations governing the excitation of waves in a three fluid plasma, (Eq. (2. 3) through (2.14), for a general inhomogeneous media with parameters depending on space coordinates, is close to impossible to solve. In the case of stratified media, one may, perhaps, with tedious algebraic manipulation, reduce these equations to a system of ordinary differential equations of higher order. General solutions for such higher order differential equations with variable coefficients are not known. Therefore, approximate numerical methods must be used to obtain the solution to such problems. A formal procedure for reformulating the systems of equations which may introduce some simplification and offer the possibility of a numerical solution, is the general operator-transform method proposed by Wu (1965). This procedure is an extension of the operator method used by Diamet (1963) to obtain formal solutions for Maxwell's equations. The procedure for the formal reduction of the equations is as follows. FIRST: For the purpose of exhibiting a general solution to a system of basic equations, it is convenient to reformulate them in the following single operator equation v. P (r)= (r) (4.1) where 0(r) is a field vector composed of the field variables such as the electric field E, the velocity field V, etc., 0(r) is the source vector containing various excitation sources such as the electric current source J, the mechanical source F, etc., and 'u/' is the system matrix differential operator relating the field to the sources. 'V contains all the properties of the medium and is a function of 48

the space coordinate r. In general, without loss of generality, the system of basic equations can be rearranged so that some of the submatrices of W are identity matrices. SECOND: Here, we introduce the generalized transform techniques, which amounts to choosing some convenient basis of representation for the solution and transforming the operator differential equation in real space to an operator integral equation in transform space. The generic summation symbol $, such as used in Schiff (1955), wvill be used, which requires that the expression following this symbol be integrated or summed over the entire range of the repeated variable. Formally, for any quantity a(r), we may introduce the following transform pair: Transform A(s) $ $ d(s, r)a(r) Inverse a(r) m c(r, s)A(s) (4.2) with the property that $ c(r, s)d(s, p) "l(r, p) and $ d(u, r)c(r, s) a '(u, s) (4. 3) The idemfactor I(u, s) comprises a Dirac delta function or a Kronecker delta and a unit dyadic, as required. To illustrate the transform pair consider a rectangular coordinate system. The real space variables are coordinates (xj y, z) and the transform space variables may be considered as (s, s2, s 3). The range of the real space and transform space variables is -co to +a). In this case a Fourier transform is appropriate and d(s, r) and c(r, s) are 49

1 -ir~ s d(s, r) - eir (2r) (4.4) i r*s c(r, s) = e Now, we proceed to the transformation of the operator Eq. (4.1). Let 3(s) and (s) be the transforms of the vectors / (r) and p (r), respectively, i. e., (S) ) = d(s, r)O(r) p (r) = c(r, sg((s) (4-) f (s) x d(s, r)0 (r) (r) = c(r, s)(s) (46) Also, we take the transformation law for the matrix operator! as v/ (u, s) = $ d(u, r)-'c(r, s) (4. 7) Premultiplying both sides of Eq. (4. 1) by d(u, r), and then substituting the expansion for O(r) as given by the transform pair in Eq. (4. 5) and summing or integrating over the complete r-space the operator Eq. (4. 1) in the real space becomes the operator integral equation in the transform space $'0 (u, s)((s) ^ (u) (4. 8) This equation has the character of a generalized integral equation of the first kind, with p(s) as the forcing function, (s) as the unknown function, and7LJ(u, s) as the kernel. lI(u, s) is a function of two composite variables of the transform space and retains all the pertinent information about the system. 50

THIRD: Because of the earlier rearrangement and diagonalization, the dyadic kernel iJf(u, s) can be properly partitioned so that the order of the matrices to be manipulated may be reduced, by introducing coupled integral equations of the second kind, which in turn may be recombined into one integral equation of the second kind. For example, we can have for Maxwell's equations l(u, s) -Z(u, s) L -Y(us ) I (u, s) Then, the partitioning of (s) and '(u) into two vectors V(s), W(u) s) (u)= (4.10) produces the following coupled integral equations V(u)=W(u)+$ Z(u, s) I(s) (4.11) I(u)=J(u)+t Y(u, s)V(s) (4. 12) Let V(s) and I(s) correspond, respectively, to the transform of the electric field and the transform of the magnetic field, then these Eqs. (4.11) and (4.12) have the generalized forms of the telegraphist's equations of Schelkunoff (1955) if s is taken to indicate different modes in the waveguide. The elimination of either the field vector V or I in the Eqs. (4. 11)and (4.12) gives the general form of the Fredholm integral equation of the second kind, e.g., V(u)=F(u)+$K(u, s)V(s) (4.13) 51

where F(u)-W(u)+ Z(u, s)J(s) (4.14) K(u, s)=$ Z(u, v)Y(v, s) are both known functions. For homogeneous media the kernel has the ideal form K(u, s)=N(s) 1 (u, s) (4.15, and the integral Eq. (4.13) can be explicitly solved as V(s)= s) F(s). (4.16) For inhomogeneous media, no such idealization is possible, but, in principle, solutions may be obtained by standard techniques of numerical analysis. This generalized operator approach, therefore, may be called a "unified" approach in the sense that within the same mathematical framework, a technique is available which, in principle, is applicable to problems involving either homogeneous or inhomogeneous media. 4.2 The Integral Equation The steps outlined in Section 4. 1 for reducing the system of partial differential Eqs. (2. 3) to (2. 7) into integral equations in the transform domain is straightforward and has been carried out explicitly by Wu (1965). A sketch of this reduction is given in this section. The matrix equation relating the fields E, V, Vi, V n, h, n, nil nn to the sources K, Q e Qi' Qn, I, F, F, is rearranged in such a way that the resulting matrix operator contains two sub-matrices which are identity matrices. This rearrangement is essential in order to simplify the resulting integral equation by eliminating some of the field variables. The resulting 52

matrix equation is given by Eq. (4.17), where S1 20 3 4 are related to the sources I, F, F., and F, while the A.. are source independent operators. - -e f n ij The explicit forms of S., and A.. are given by Wu (1966). — ~1J 53

0 0 0 *vxl (44% 0 0 ~00 0 o 1 0 0 o 0 1 0 iN 0 (V to =? I 1) + i *w 1 0 to 0 0 0 L4L iQ e iQ n to 0 0 iN 0 (V =1 I 1 )+ +VN I to 0 n.I n 1 0 0 0 1 0 0 iN 1. 0 -(V. 1)+LVN *1 to to I U-l A11 A12 A13 A14 A21 A22 A23 A24 A31A32 A33 A34 A41A42 A43 A44 0 0 0 0 1 0 0 0 vi Vn Si -2 0 0 S 3 0 0 0 S -4 (4. 17)

The matrix Eq. (4.17)can be put into an operator form as z/' (r) = P(r) (4.18) where h -iK e n. iQe (0 n iQ. 1 (r) - E ( V n(r) i n (4.19) V, S -1 -1 V S -2 Equation (4..18) can be considered as an abstract relation between the sources and the resultant fields. q (r) is the eighteen-vector representing the field quantities, p (r) is an eighteen-vector representing the source quantities, and 1/ is the system matrix differential operator relating the field to the sources. The generalized Fourier transform, such as defined by Eqs. (4.2) and (4.3), may be used to obtain the transform of Eq. (4.18). The result may be formally expressed in the form of Eq. (4.8) which is repeated below. 55

$ W(u, s) (s) = f (u) (4.20) Equation (4.20) may be put into the generalized forms of the telegraphist's equation by partitioning the transform of the field vector, {i(s), the transform of the source vector, fI (s), and the transform of the matrix differential operator, bA (u, s), as follows: h It(s) n V (s) e e n. V.(s) 1 1 n V (s) 2(s)=$ d(s,r) n n (4.21) E Vt(s) v I (s) e v. I.(s) V I (s) - n.n 56

-3 Jt(s) i&e ) W,(s) iQn 1 W.(s) S W(s) 2 e(s) S3 J.(s) 3 1 S J (s) L 4 n where It(s), Vt(s), I (s), Ii(s), I (s), Jt(s), W (s), J (s), J.(s) and J (s) are t t e 1 n t t e 1 n three by one column matrices, and V (s), V.i(s), V (s), W (s), W.i(s) andW (s) are scalars. In view of the orthonomality of the transformation kernels, c(r, s) and d(s, r), the eighteen-dyadic kernel '/(u, s) can be partitioned as 57

I, - - I I O 0 N 0 O N N Go 0 CD "I0 Co j4 02 "'-4 0 0) 02 0s 02l 0) CD0 - 02 0 0 0 = '-4 02 H4 - 0 0 0 02% I 0 02 0) 0 % ' —4 E —.02 0 co I ) 0 0 H a1) a) H go a) 011%.,-% m m ft -%:3 =1 0 >.4 >-I I I % 02 It — 4 02 0 0 0 02 02 H H$ I I =I 58

The elements in these equations are given explicitly by Wu (1965). By substituting Eqs. (4. 21), (4. 22) and (4. 23) into the integral Eq. (4.20) the following generalized telegraphist's equations can be obtained: I (u) = J (u)+$ Y (u, s) t(s) Vt(u) = Wt(u )+ t(u,) It( s)+$ Tte(U, s)V (s) t t t t Lte e (4.24) (4.25) +$Tti(u, s)V.(s)+$ T (u, s)V (s) ti 1tn n I (u) = J (u)+$ Y (u, s)V (s)+$ T (u,s)It(s) e e e e et t + ei.(u,s)V.(s)+$Y (u,s)V (s) ei i en n V (u) =W (u)+$ Z (u,s)I (s) i(u) = J Y( )V) + TY (u, s)I V(s)+$ Tit(u,s)I(s) (4.26) (4.27) (4.28) +$Y. (u,s)V (s)+$Y. (u,s)V (s) le e in n Vi(u) =Wi(u)+$ Zi(u,s) Ii(s) I (u) = J (u)+$Y (u,s)V (s)+ Tnt (u,s)It(s) n n n nt (4.29) (4.30) +$ Y (u, s)V (s)+$ Yn (U, s) V.(s) ne e 1 Y V (u) = W (u)+$ Zn(u,s)In(s) n nn n (4.31) 59

By properly partitioning (s), f (s) and -A (u, s) as given by Eqs. (4. 21), (4. 22) and (4.23), the basic equations can be reformulated into the general form of the Fredholm Integral Equation of the second kind. First of all, the transform of the field vector is partitioned into three column vectors each with six components as follows:;l(s);2(s) 3(s) (4. 32) where It(s) v (s) e v.(s) V (s) n Vt(s) 02(s)- - (S)!As) I(s)3 (s) - Lns-. L Similarly, the transform of the source vector is partitioned as three six-column vectors (4.33) where 1(s)= J (s) te 1 e w (s) 60 Ji(s) - 3(S) ' J (s) nj

Next, the transform of the matrix differential operator as given by Eq. (4.23) is partitioned in the following form I(u, s) -12(u, s) - 13(U, s) /' (u,s) -UZ 2(u,s) 1(u,s) 0 (4.34) -t1(U, s) 0 1(u, s) Substitution of Eqs. (4. 32), (4..33) and (4. 34) into the integral Eq. (4.20) gives three coupled integral equations which are () - (u)+$ W12(u, s)2(s)+ 24;3(u. s)(3(s) (4.35) 22(u) = f2(u)+ $t 21(u, s)1(s) (4. 36) p3(u) = f3(u)+$ /31(u, s)l1(s) (4.37) where Yt(u, s) 0 IAo' ( z (u, s) 12(u, s) e (4.38) 0 0 0 0 0 0 0 0 IA 13(u' s)- (4.39) Z.(u,s ) 0 0 z (u, s) n 61

Lf (u, S)= 21 Zt(u, s) Tet(u, s) Tit(u,s) nt(U, s) Tte(u, s) Y (u, s) e Y. (u,s) ie Y (u, s) ne Tti(u,s) Tt(u,s) Y i(u,s) Y (u,s) (4en 40) (4.40).1d (u, s)~ 31 Y.(u, s) Y. (u,s) 1 mi Y.(u, s) Y (u, s) ni n (4.41) and (4. 37) into Eq. (4..35) gives Finally, the substitution of Eqs. (4.36) the desired Fredholm integral equation of the second kind for the field variable Jl(s) as 1(u) = F(u)+ K(u, s)!l(s) (4.42) where we have defined F(u) I (u)+$ Z42(u, 8s)2(s)+$ Z13(U s)3 (S) (4.43) and K(u,s) 1$ Z4(u,v) Z1(v,s)+ 13(u ) 31(V, ) (4.44) which are both known functions. Thus, the order of the matrices has been reduced from 18 x 18 to 6 x 6, representing a considerable simplification of the problem. No attempt has been made to obtain numerical solutions to the integral equation for a general inhomogeneous medium, however, application of the results to a homogeneous electron plasma has been illustrated by Wu (1965). 62

V EXCITATION IN A HOMOGENEOUS PLASMA 5.1 General Formulation The general operator formulation given in detail by Wu (1965) and outlined in Section IV, is, in principle, applicable to excitation problems in both homogeneous and inhomogeneous media. For the case of a homogeneous, unbounded media, of course, the set of equations given by (2.18) through (2.25) becomes a system of partial differential equations with constant coefficients. In this case, the direct use of a three-dimensional Fourier transform will reduce the system of differential equations in real space to a system of algebraic equations in transform space. Thus the Fourier transform of the field components can be expressed in terms of the Fourier transform of the sources. Explicitly, if we define the Fourier transforms of the sources and fields by Eq. (4.1), and choose the coordinates such that A iS ZS and B3o (z cos +ysin0)BO the following set of algebraic equations, relating the fields to the sources in transform space, is obtained. Maxwell's Equations; -sE -q( h =iK (5.1) y ox x sE — op h =iK (5.2) y oy y 63

-oAu h =iK 0 z z (5.3) -sh +eoE+ieN (V. -V )=-iI (5.4) y o x o ix ex x sh+wE~E+ieN (V. -V )=-iI (5.5) y o y o ly ey y wEoE+ieN(V. -V )=-iI (5.6) o z o 1z ez z Electron Motion: Ns Q o e n - "V i- (5.7) e w ez o iF (l+i + iv )V +i2 V cos0-in V sin0+ eE -iv.V. - i V = ex ei en ex e ey e ez wm x ei ix en nx wNm e o e (5.8) iF (1+iv.+i )V -ifQV cos+ ie E -.V. -iv V = eY (5.9) ei en ey e ex Wmn y ei ly en ny ON m e o e 2 U iF (1+iv +w )V + iV sin -- en+ ieE -i.V. -iv V - en (5.10) ei en ez e ex w N e wm z ei iz en nz wNm o e o e Ion Motion: N s n.- V =i- (5.11) 1 W 1Z W 64

iF. ie iF. (l+iv. +i. )V. -X.V. cose+ifV. sine — E -iv. V -iv V = ie in ix i ly i iz wm x ie ex in nx wNm. e o 1 (5.12) ie iF (l+iv. +iv. )V. +fi.V. cose- E -iv. V -iv. V -i = (5.13) ie in ly 1 ix om y ie ey n ny wNm. 1 001 2 U2 iF. s i ie lz (1+iw. +iv. )V. -iS.V. sin -- n- E -iv. V. -iv V (5.14) ie in iz 1 ix N i w. z ie iz innz wNm 0 1 0 1 Neutral Motion: Ns Q n - nV =i - (5.15) n u nz u iF (1+ +iv+ )V -iv.V. -iv V = (5.16) ni ne nx nl ix ne ex wNm o n iF (1+i.+iL )V -iv.V. -iv V =- (5.17) ni ne ny nl ly ne ey wNNm n n 2 U iF (+i +i )V -iv.V -i V - n = nz (5.18) ni ne nz ni iz ne ez w N n UN m o n n From the above equations, the transformed field quantities can be expressed in terms of the transformed sources. Formally, the system of Eqs. (5.1) through (5.18) can be represented by the matrix equation [L(s f(s) = S(s), (5.19) 65

where [L(s)] is the matrix of the coefficients of the equations and f(s) and S(s) are column vectors representing the field quantities and the sources respectively. The determinantal equation, L(s)= 0 the condition for the existence of fields in the source free region, has been discussed in detail for the collisional case in Section III. The algebraic procedure involved in finding the inverse of the matrix [L(s)], and then carrying out the inverse Fourier spatial transform to obtain the fields excited by various sources is straightforward, but, nevertheless, extremely tedious. Explicit expressions for the fields are generally obtained only for some ideal problems. The most general case carried out explicitly to this date is, perhaps, the work reported by Wu (1965) for the problem of excitation of waves in a collisionless electron plasma. 5.2 Collisionless Electron Plasma For the case of a collisionless electron plasma, which is, perhaps, an adequate idealized model for the F-region of the ionosphere, the algebraic equations are considerably simplified. In this case, Eqs. (5. 1) through (5.10) form a complete system, so that the inversion of the matrix is relatively simple. If one carries out the tedious algebraic steps, it can be shown that Ei G..K.+H..I.+L..F.+M.Q (5.20) i G.i ij j ij j ij ej iej where i and j have the values x, y and z. Explicitly, the quantities in (5.20) are 66

s6 2 29)4 22 2 2 1 2 2 2 22 G= 2 Cos' es4{ 12-cs+jj x1 -2cos0e-Q sin ij e o o e o 2 1 2 22 21 ~1 2 2 2 2 2 2 2~ -S y2L L4 CosJ +- (- IW- oe-2wPsinO e 0 G, I= ecos0aii1+iai2 0 G -- ia +flcose a +0 sin0 12 (3 Lii e 12 e 13J 0 G13 ~ G~1=j [~cos6 a2l+ia22 0 G22- S i a +SIeco'soa2211 sine a23 =0 67

31 2 e 31+ia32 ] 0 32 2-[ia31+ cosOa32+2 sinOa33] Ge 32 e 33 0 G O G33 0 H11le (iall-2 cos a -0 sin a3) I I (i11 e 12 e 13 o H e sO all iaa2) l3WE e 11 12 0 H13 - si-eal+i(1-2 )al3 ] O e H i1 -52 cosO a -.sin a 21 w-' 21 e 22 e 23 O 1 H 3-(Qcos a ia 22 WE e 21 -a 22) 0 1 2 H 2%sin a2 +i(_ 1s)a3J 23= -we e 1 0 O e 68

H 31 e 32 e H a cos a 21- sina332 H32 we e 31 32 0 2 - s_) H [ sina1 i( 2)33 a33 lj e 1 W1 0 e e = - a1 eo e L2 - 2- a22 22 Wm e 33 e 33 eo L..O, i0 j 13 eQe 2 s a13 e o eQe 2 o 23 69

eQe 3= 2 b 33 13 WE e o where 4 s a11 2 22 o e - s2 [( I ++i~x1-W2)1 +-2)2 0 e a iq cos____ 2 — 12 s f202 o e 1 -V1(- 2 I W e 0 eo3e e 8 4 2 ri 1 2122 2 2 a22o cs L -S+-Xl- w)-jisi9+(1-w ) ) i o e 0 e 0 2 22 azls inG cosG (~- 1) 0 70

a31 2 1e2f 2 a 1- sine -s -+- (1-a)+(1- ) 31 e [22 L2 2 p p eo e o 2 s2 a32 e sine cos (- 1) e 2 o a33 - (- cos 9) — (1-u - cos e)+(1-w ) -n cos e 33 4p e p e o o Using the expressions for the components of E, the components of h may be obtained from Eqs. (5.1) through (5. 3), while the components of V and n may be obtained from Eqs. (5.4) through (5. 7). The Fourier inversion of Eq. (5.20) to obtain expressions for the fields in real space, even for a unit impulse source, is quite complicated. Numerically, the asymptotic expression for the field can be evaluated by using the methods of stationary phase. Examples of such a procedure have been carried out by Wu (1965), where the electric field, due to current sources in the direction of the d. c. magnetic field, has been evaluated for numerical parameters appropriate to various regions of ionosphere. The numerical results obtained by Wu seem to indicate, and confirm the belief of several previous investigators, that a substantial amount of energy is excited in the form of plasma waves. However, from the numerical values of the propagation constants of the waves (Section 3. 3) indicate that the plasma waves, when collisions are not neglected, attenuate quite rapidly. 71

Further detailed calculations are necessary to clear up this point, since this fact is very important in the investigation of radiation from current antennas in the ionosphere. 72

VI SOURCES IN A BOUNDED PLASMA 6.1 Statement of the Problem The following presentation is an attempt to describe the waves that might be generated by moving and stationary sources located in an ionized gas. Only the simplest model of an oscillating dipole is considered in any detail because of the complexity of the equations and the lack of time. It is felt, however, that the methods can eventually be extended to include more realistic and practical situations. For example, a moving vehicle can be considered, as a crude approximation, to resemble a charge moving along a prescribed path. Inclusion of this effect is simply a matter of appropriately selecting the source terms to describe the motion. In addition, the general formulation has been presented in a manner which will allow for spacial variations in the unperturbed charged particle density. An analytic description of the behavior of the fields when this feature is included has not been obtained in terms of known functions. It is felt, however, that the method is definitely amenable to numerical treatment, and this aspect of the problem is one which is definitely worth pursuing. In this section we shall discuss the behavior of a macroscopically neutral plasma which is bounded by a perfect conductor in the interior region and which extends to infinity in the exterior. Since neutral particle effects will be ignored, we shall be concerned only with the ions and electrons, whose motions are governed by the first two moments of the Boltzmann equation together with the equation of state and Maxwell's equations. Certain simplifying assumptions will be made to assist in the analysis, namely that Landau damping and shielding effects may be ignored, that there is no external electric or magnetic field and that collision terms may be neglected. On this basis we can obtain two scalar equations whose solutions can specify completely the behavior of the plasma. 73

The perturbations which take place in the plasma will be considered as being produced by mechanical, electric or magnetic sources imbedded in the plasma. These sources, which are located a finite distance from the conductor cause the particles to move in some manner. A description of their motion is affected by the presence of the conductor which requires, then, a prescription of the boundary conditions to be satisfied. This is, in general, a very delicate matter and has never been well defined to everyone's satisfaction. In most cases, however, the conditions used by Cohen (1962), namely that the tangential electric field and the normal components of the particle velocities, and the magnetic field vanish at the surface of the conductor, are acceptable. Although these assumptions are not rigorously valid, they do have the advantage of simplicity. In addition, any waves stimulated by a source must satisfy some form of the radiation condition. A more exact statement of this behavior will be discussed when we are dealing with the appropriate Green's functions. In the following three sections we shall discuss the manner in which the wave solutions may be obtained. The first section will be fairly general in that it will allow for variations in the particle densities, as well as for arbitrary sources and conductor geometry. In the second section, we shall consider one of the simplest problems, namely the oscillating electric dipole source located above a perfectly conducting plane. In the final section we shall investigate a somewhat more complicated problem involving the electric dipole in the presence of a conducting sphere. In both cases we will observe that the sources give rise to a transverse electromagnetic wave and two longitudinal plasma waves. 6. 2 The Potential Functions The situation we are considering involves the behavior of the particles of a fully ionized gas in an infinite half-space which is bounded by a conductor. Denoting the mechanical, electrical and magnetic sources by F(r), S(r) and J(r), respectively, the Fourier transformed equations of motion and continuity together with Maxwell's equations 74

2 e -iwV +u Vn -N E +F - e e -e m oe 2 e -iwV.+. Vn= N E +F. - i n- i m. o- -1 1 1 1 Vn =-V(V. V )+-VS -e iw - - e w - e Vn.- V(V V.)+ VS. - 1 vL -1 - 1 i10 VxE - H - - C VxH=-E+ (V. -V )+J C- c c -i -e c - (electron motion) (ion motion) (electron continuity) (ion continuity) (6.1) (6.2) (6.3) (6.4) (6.5) (6.6) The subscripts (e) and (i) indicate terms relating to electrons and ions, respectively. The mechanical quantities F and F. represent forces per unit mass, S and S. are the rates at which electrons and ions are introduced into the plasma, e 1 and V - N v and V. =N v. are functions introduced for shorthand notation. -e o-e -1 o-i Substitution of (6. 3) into (6. 1) to eliminate n yields e 2 2 u u e -e e -iwV + V(V- V )- N E+F — VS -e ijt — -e m o- -e iw- e e while similarly, from (6.2) and (6. 4) we have (6.7) 75

2 2 U. U. - iwV.+ + V(V V.) e N E+ F. - VS. -i it - - - m. o- -i li - i 1 From these two equations E may be eliminated to give us (6.8) 2 2 U U. m [-ijV + e V(V' V ) +m -iV. +-i V Vi.) e -e i -- -e L - l - - 2 2 Lm u m.u. e eV + i I VS. i -e it -i = (m F +m.F.) - e-e 1 — (6.9) Taking the divergence of (6. 6) and rearranging terms, we obtain 47r iLO V- V. ^-V -- V.J+- V. E - -i - -e e - - e - - which, when substituted into (6. 9) yields (6.10) 2 u m iwV + e V(V. V )]+ e -e iw - - -e 2 Cu m - i V i+ - i -i iw v(v v )- - -e 47r it eV(V J)+ -V(V e - - e -- = (m F + m.F.) - e-e 1-i 2 2 m u m.u. e e VS + i VS. iw - e io - i Finally, from (6. 5) and (6. 6) we have 2 V. = V + E+- Vx(VxE)- J4, -1 -e e - iwe- - e (6. 11) (6.12) 76

so that V. may be eliminated from (6.11) to give -1 2 2 m u +m.u. -iw(m +m.)V +.... V(V V )= e i -e lw - — e 2 2 m.c { 2 ui c-iI { E-(VxVxE)+ Ve -" - - - 2 -C 1% C c 2 4iriwm. 4rm.u. J+- i i V(VJ) + e iue - - 2,m u (m F + m.F.) - ( m VS + e-e 1-1 \ 1 - e 2 m.u. 1W- - 1 i(M — 1 (6. 13) In (6.13) and (6. 7), then, we have two equations which are similar in form and which may be written in tensor notation as e L. V =g e (6.14) M. V =h. where 2 2 2 m u +m.u. 2 L. -i(m +m i)6 i axax (6.15) 3 a 2 m.c g = e oj e 2 a2 ax.ax J a 2 4jriwm. (VI- + ) 6 E - e- J. + 2 jC a e J c 2 4rm.u. I 1 iWe 2 2 2 a2J. m u as m.u. as. a +. 1 )-( ee e e 1 1 Fe+(m F.i - +F ' ) a x.a x e j i j iw a) x. iw a x. j J j (6.16) 77

2 u 2 e --., Mj= -iw6j6, a d (6.17) j] ]3 i x a x1ax and 2 u aS h e e e e j m o ja a j i ax.6.18) e ] Since L. and M. commute, ja ja Mjg (6.19) From (6.15) through (6.19), then, we obtain the vector equation 2 2 2 2 iL(m +m.) e N (m u +m.u.) 2 u e 1 0 o ee 1i e e- E e-V(V N E) -iw(1 —)V(V- E)+ m m. 2 - m m. 2 - 2 e i c e 1 c c 22 2 2 u u. u i (V2+ -)E e ( V(V2V E)+ i — e V(V- E)+X+VR (6.20) c iwc c In this relationship we have introduced the source terms 2 W iwoe X -4- J ----(m F +m.F.)+ (m +m)F (6.21) - 2 e-e 1-i 2 e i -e c imc - r.c 1 1 and 78

22 47r 2 2 4rue U. R(r) = - (u +u )V. J c w C 2 eu e V2 VJ+i e -- 2 iwm c i m V F +rmV' F + e- -e 1- - e 2 2 (m u S +mur S.)+ 2 e e 11 1 m c 1 2 eu e 2 2 rm.t c 1 2 2 2 2 mu V2S +m.u 2VS. ee e e 1 i (6.22) We now assume that the electric field is composed of a longitudinal and transverse part, to wit E= -Vp +- VxZ -V + -A -C C- - - H= Vx(Vx Z) -- VxA and (6.23) (6.24) With this notation, (6. 20) becomes 2 iw(m +mr) e N e 1 0 o+ IW m m. 2 [-V(N p) +(VN o) + -A m m 2 - o - o c e i c 2 2 m u +m u. 2 e I c 2 2 u +u. iw(1- e 1 0 2 C [- +i A] + 2 2 u u. e V1V4 +X+VR iWe (6.25) 79

The gauge condition is thus 2 2 2 iL(m +m.)e N m u +m.u. 2 t e - o2 e e e i ie N -N + N A mm. 2 im m o o c o z e 1 c e 1 c 2 2 u u. 3 e 4V iW( 2 2 2 (i6 2 - (u +u) 2 +R(r) (6.26) iwc c c where we have assumed that N varies only in the direction normal to the conducting surface. The remaining part of (6. 25) yields iw(m +m.) 2, m +m. 2 2 2 2 e i e_ e m e m e -N n0 - e i w eN A - (V+ )A+X (6.27) mm. 2 o m. 3 o- c - - e 1 c e 1 c c In these two equations we observe that coupling of the longitudinal and transverse waves will occur if N $ O, so that the two types of wave propagation are independent only in a homogeneous plasma. In addition, it can be seen that R(r) and X(r) as defined in(6.22)and(i21)are, essentially, sources of longitudinal and transverse propagations, respectively. It is evident that('6.26) and(6,27) canbe uncoupledwith little difficulty but that solving the resulting equations analytically with variable N can be done only for the simplest of cases. In the following discussion, then, we choose t N r0 and write (6.26) and (6.27) as O 2 4 b 2 c 2 2 2 i2c V4^+b V2+ C + =(V2+a )(V2+ )= - 2 R( (6.28) a a 2l 1 1 u U e 1 and (V2+a2A = X(r) (6.29) LO 80

where 2 b 1 [ b 2 4c /2 02a 2 La a] 1 1 1 (6.30) (6.31) (6.32) 2 2 w2 and a = 2 C m +m. e 1 m m. e i 2 e N 2 C If we assume that the electron kinetic energy is at least as great as the ion kinetic energy, say 2 2 m u Xm.u., ee e 1 1 > 1 (6.33) then from (6. 26) 2 2 m +X)e N b — i(1+ e (-+X) o a 2 1m. 2 m 1 u. 1 u e 1 e 2 2 Ui 1 [ (+x) x 2 m w e pe l m. 2 1 u (6.34) and m 4 c e o - a1 Xm. 4 1 iU. 1 2 2 w +. m 4 pe pi ]_, e ) 2 - m. 4 1 U. 1 2 - I -2 [- 2. (6.35) b c e Since (-) is greater than (-) by the order of - a1 a1 Xm1 (except when 2 e+ e 2 w=i (1+X) --- ) m. pe 1 81

2 2 c/al 2 (2 - w 21 e ](6. 36) e e Ue Xm. 2 1 1 I mi 22 2 b 2 r (l+X)me p 1 a 2 L XM 2I (6.37) 1 u. I ~ 1 2 1 [2 2+ 21 pe (6.38) U -M W 1-( + ) )~(6.38) 2 pe pi J 2 2 (In subsequent sections, we shall delete the subscript (e) on - since 2 2 2 2 pe 1w w +W. =.w.)p We thus have, in (6.36-6.38) three propagation constants p pe pi pe which vary inversely as u, u. and c respectively. 2 e e Turning now to the boundary conditions, we seek to determine what effect the vanishing of n- V. = n- V = 0 will have upon A and p. From(6. 7)and(6. 8) we obtain v(V v)+ v NE+ F -VS (6.39) 2 -e 2 o- 2 -e - e u mu u e ee e and 2 1 i10e i10 V(V- V.)+ V. N E+ F. - VS. (6. 40) 2 N 2 o- 2 -1. u. m.u. u. 1 1 1 1 Subtraction of (6. 39) from (6. 40) yields 82

v[v. (r - ] -1 -ie V1 2 + 2 IL -e ]1 + 1 r L2 +" 22 u m u m.u. 1 e e e I N E + o - -F. F - *iw 1 -e L-'2 2- 2 u. U 1 e - (vs.-vs ) - 1 - e (6.41) From (6. 10), however, V[V (V. -V )] - -i -e = i- V(V E) - - V(V- J) e -- e - - - (6. 42) which, when substituted into (6. 41) yields 2 u -r -_ 1 1 =TW 2Iie N E [|2+ 2 O-2 u. u m u m.u.! e e e 1 1 - i V(V. E) +4- V(V. J) + e -- - e -- - F. F - 2 2 u. u i e - (vs.- vs ) (6.43) Finally, V. from (6. 12) may be used in (6. 43) to obtain — 1 2 1 2((2 1 - =2)V We + u nt.U m u e i ee N E -V(V- E)+ e — Vx(VxE)- - E + o- e -- - e 2 - - 2 - U. c 2 - (V -Ve e 2- (6.44) U. 1 U^. V(VY J)+ilu 1 e - -2 Lk F -^ 2 u e 83

Similarly using (6.12) to eliminate V in (6. 43) results in e 2 2 21 1 1 it it C ( - )V.i ie 2 - N E — V(V-E)+- Vx(VxE)- E + 2 -i e - - - e - - -. u r.u. mu u c 1 e 11 ee e e(- — W -_(VSi_ VS )+ 4_ e- 2 2 - VSe) (6.45) u. ue u 1 e e If the normal components of V. and V are to vanish at the surface, then the - -e two following conditions must hold on the boundary W 2e 2 N E - (V-E)+- - --- -- 2 2 o n e an e an -2-(S.1S ) m.u. mu u ii ee 1 e (6.46) and 2 2 in' [Vx(VxE) - E]+ 4, Jn 0 (6. 47) c c It is a simple matter to use (6. 23) to obtain two mixed boundary conditions for A and 0 at the surface. The third condition, namely that (nxE)= 0 on the n conductor, merely requires that a+iw A O. 0(6.48) at c t while for the final one, n * H 0, we have n- (VxA) = 0 (6.49) 84

6. 3 Vertical Oscillating Dipole Above A Perfectly Conducting Plane For an oscillating vertical dipole we have J i J = erlw sinw t6(x)6(y)6(z-z ) (6.50) Thus, X in(6.21)has only a z- component which we denote by Z(r), so that (6.29) may be written (V+a2 )A = 0 x (V+a2)A - 0 (6.51) y (2+a2) A -Z (r) z 2() The Fourier transform of (.51), assuming that a variable u(r) transforms according to oo u(k, z) u u(r)e-ik- rdx dy (6. 52) gives us A -(k -O)A X 0 (6. 53a) x x A -(k2 -a)A 0 (6. 53b) Y Y A -(k2 2)A c Az -(k - )AZ 2 Z z) (6.53c) Likewise, from (6.28) we obtain rd 22 2- -d 2_ 22 I Lp. [ -(k 2-2 [ -k2- R(k z) (6.54) dz dz Ue U 85

The variables A and A must satisfy the radiation condition. Thus x y A (k, z)- c e x 1 - (k -a z A (k, z)= c 2e Y (6.55) The condition n. H x 0 requires that k c = k c This, however, means that n 'H O, and since this adds nothing to our solution, we may choose c1 c2= 0. It remains, then, for us to solve (6. 53c) and(6. 54) forA and0. Takingthe first of these, we consider the inhomogeneous equation 2 d - (k2-a)g 6(z-e) dz (6.56) which satisfies the condition g( 0, )-=0 and the radiation condition. Thus 2 2 sinh {k - z e Tk2 2 z k -a (6.57) 86

Then 2 c - k-ac z A (z, ) C g(z, ) Z(Q) d + 0(O)e Z ) ( Jo (6.58) where the constant q(0) must be determined from the boundary conditions. For the solution of (6. 54) which satisfies homogeneous boundary conditions, we consider the auxiliary equation 2 2 LG= [ _2 - a-2 -(k2 -2)] G( Z)6-) dz dz which satisfies the radiation condition and the requirements G(0, ~)= 0 (6.59) a (0, ) 0 az (6.60) G(z, g) is constructed from the four solutions of L= [ d (k2-a d (kz_2 2)] 0, dz dz (6.61) namely 1= e 02= e P3= e k2 2 ' -m z (6. 62a) (6. 62b) m z e _ - m2Z me (6. 62c) 87

_1 2 mz 4=e = e (6. 62d) In order to interpret the radiation condition, we note from (6. 54) and (6. 59) that we obtain 2 [CLG-GL] =: 6(z-0)+ 2 2RG (6.63) 2 u. U U. e 1 Integrating from =0 to co, we have Co o O02 O O U U. X0 ( 2 G(e, z)R(f) do + [p (r i)L G-GLg P()] d5 (6. 64) e i o The final term in this equation will vanish if G(Q, z) varies as a linear combination of pl(z) and p3(z) as z - oa. Thus al P(z)+ a202(z)+a3 3(z)+ a404(z) z G(, z) = (6.65) b l(z)+b3P3(z) z> This, then, is the form which satisfies the radiation condition. Using the continuity requirements at z= e and the conditions stated in (6. 60), it can be shown that 88

G(z, )= 2 L(m2-m2 )e 2(m2 -m ) (ml+ m)ml m2 m z -mlz.L -m2(ml+m2)e +2m m2 e -m2z] -m -+ m2z e + (m -m2 )e -ml(ml+m2)e -m z -m z - m2 +2mm2e je (6.66) for z <~. When z> E, the roles of z and E in (6. 66) are interchanged. Hence 2 - WC ~(z) 2 22 U U. e 1 oD 0 (6.67) is the solution of (6. 54) which satisfies the homogeneous boundary conditions. In order to satisfy the inhomogeneous boundary conditions we merely add two functions such that 2 iwc P(z) 2 2 U U. e 1 o where Gl(z) and G2(z) satisfy (6. 60) and Gl(0)= 1 GI(0)= 0 G (0)= 0 G2(0)= 1 2 (6.69) (6.70) 89

The three boundary conditions in (6. 46-6. 48) require that 2 2 p'"(0)- (1+x) 2 +K] '(0)+ - (+X) 2Az(0)= u u e e 2 ( c 2,(o)+ A (0()- - Z(0, k) c 0 O(0)= 0 (6.71) (6. 72) (6. 73) The last of these allows us to set p(0) on the right hand side of (6. 68) equal to zero. The problem remains, then, to determine the other two constants. To do this, we consider 0(z) when z <z. From(6.22), (6. 50)and(6.68)we have 0 -m z 2y erw wf (u+W )-6 (u-u ) 2 m itzeI 2 2 om J u 2 r.fI+ \ e -- L m 2 1 + (1 me 2 + 0) [e- mz m2z m ml 22 J nv- m (6.74) For simplicity, in(6. 74)wehave let G(Q, z), as specified in(6. 66), be denoted by G(, z) = t (z)e -ml +t2(z)e -m2e (6. 75) 90

To evaluate f 1 t (O) we note that t (O)~ 2 1 mI+ m2 mI '2'' (O) =m n12 (6. 76a) (6. 7 6b) and Thus,, It2 2 1 1 2 2 ~ 0 2 U. 1 2 L -m, n~ w2r,+n r - z-r Iz Jl+Xmt em2zo e lo 22 { Me 2 = =(M1+ M IM2+ n~) () ~~ 1 m LO 1 W (6. 77) With the aid of (6. 77), the conditions in (6. 71) and (6. 72) reduce to 2 2 m 1+MI m 2+m2 1 i2 ho *(1+X)co2.4 c 2 U e 2 (0 2 c A (0) z - F~w, k2, z ) 0 -.j; Z(0, k2) 2A;.78) (6 91

For the dipole at z=z, it isclearfrom(6.21)that Z(,k 2)=0, which will 0 simplify matters somewhat. If we denote the determinant of the matrix in (6. 78) by A, we have 2 2 2 _2 _ —a 2 + k2-a)k( -2 + (1+ (6.79) 2+(l+k) (6. 79) 2 2 2 c u. u 1 e Thus, for example, 2 F(w,k, z ) A (0)=' (6.80) A complete analytic solution with this rather complicated expression is very difficult and has not been successfully obtained so far. We proceed, however, to obtain a formal solution which may be evaluated numerically if necessary. For z < z as we are assuming, 0 r - sinh[ (k2- ) z] A (z) = -4re nwo [6(W+- )-6(0i-w x z 0 i o o2 2o \k -a - k2 -a z -k -a2 z e +A (0)e (6.81) z Performing the spacial inverse transform first, we have OD A (w, r,)= 2 A (w,, z)ek -dk dk (6.82) z - x y 47 92

2 Since k appears in A (w, k, z) only as k, it is more convenient to let k = kcos t x k = ksint y x = rcos P y = rsin p (6. 83a) (6. 83b) (6. 83c) (6. 83d) so that k r = kr cos (p -t) and (6.84) (6.85) dk dk =kdkdt x y Thus, (6.82) becomes o00 2r A(, r z)= I kdkA (w, k2, z) ( e co(dt -) 4w 0 0 O0 Az(,k,z)J(kr)kdk (6.86) 93

Substitution of the expression (6. 81) into (6. 86) yields -(z -Z) 2k-a2 0 A (w, r z)= -er") 1 (w*w )- 6(w-w z oi 0 0 -(Z+Z)e k-a, e- X -22 J Mk -a 2 2 ' k _ a o J (kr)kdk+ A (0) e J (kr)kdk: -er W~ - 0 x + 2 2' ia r +(z -z) 6(w+w ) -6(-)]ie r +(z -z)2 ia r +(z +z) e ~r +(z+z) 0 00 2 A (OM)e J (kr) kdk 2r z o 0o (6.87) To obtain the time dependent expression for A (r, z, t), we perform the inverse transformation of (6. 87) which yields 94

ia rr +(z -z) ia r +(z +z) n tis eq ion, o o A (r, z, t)= erlL sin t -+ 21 0 O 2 22 co o iL k r+(zz r+ ( z 0 0 In this equation a. 4 -Vet (6. 89) O C p is the plasma propagation constant for this two-fluid plasma. The first two terms in(6.88), then,correspond to the spherical waves which are generated by the source dipole at z=z and its image dipole located at z= -z. The third term is more difficult to analyze exactly, owing to the complicated nature of the integrand. 6. 4 Vertical Oscillating Dipole Above a Perfectly Conducting Sphere As in (6. 50), we have 6(r-r )6(o) 6(r-r )6(6) J= i erlsinw t = i ersin ot 0 (6.90) -r 2rr sin 0 ~ 2wrr sin 8 where we have assumed the dipole to be oriented along the z-axis. This provides us with cylindrical symmetry and it develops, as in the previous section, 95

that A = A = 0. Thus, we again have to deal only with A. There will, then, be r and 0 components of A, which are given by A = A cos 0 r z r z (6.91) A A sin 0 0= z We introduce now the Legendre transform 7r T = 2r f de sineP (cose ) n o (6.92) such that T: f(r, 0) = f (r). The inverse of T, applied to f (r), is n n 00 -1 1 ^'2n+1 V P (-cose T f (r)= f(r, 0 )=- 2 - P (cose)f (r)=-1 v 1/ n 2w 2 n n 2 cos I v n=O c For the purpose of solving (6. 28) and (6. 29) we note that T: 2fr r 1 d d 2 d n(n+l) f (r) r r and 4 [1 d 2 d n(n+l) [ d 2 d n(n+l (r) T: Vr 2 dr (r) r r r r found to be )f - (r) IV 1/2 (6.93) (6.94) (6.95) 96

We thus have to solve the equations 1 d2, n(n+r d 2 d )+ 2 n(n+l 1 d l) i 2 dr dr 2 2 dr (dr 2 n 2 2-R (r) rr r u e 1 (6.96) and d 2d 2 2 n( + 1c dr(r )+ - (r)= X (r) (6.97) 2dr dr L~ o n W n,r r where we will let /(r, 0) = A (r, )) in order to avoid subscript confusion. z The four independent solutions to the homogeneous equation corresponding to (6. 96) are W h (1) (ar) p2r (2) (ar) P2 n P3(r) hn(1)r) P4(r) n (3r) while for (6. 97) we have ql () n (a r) (2) (6.99) q2 (r) =hn (ar) 97

The constriction of the Green's function for (6. 96) which satisfies the two conditions G (a) = 0 n (6. 100) dG n _ (a) = 0 ar is a straightforward, though complicated, process. It develops, after several pages of algebra that will be omitted, that 2 G (Q, r)= - )Xr p(a)-p3)+(a)p t 3(r)p (6.101) n 4(2 2 )[p(a)p (a)-p3(a)pa (a) I 1 3J for <<r, where Xl(r) [p2(a)p (a))-p3(a)p (a)] p1(r)- [p (a)p5 (a)-p3(a)p(a)] p2(r) + [p(a)p2 (a)-p2(a)pI(a)] p(r) (6.102) and \3(r)= [p3(a)p4(a)-p4(a)p3(a)] p(r)- [pl(a)p4(a)-p4(a)p;(a)] p3(r) + [l(a)p3(a)-p3(a)pI(a)] p4(r) (6.103) When e > r, the roles of e and r in the bracketed expression of (6.101) are interchanged. 98

The corresponding Green's function for (6. 97) satisfying g (a)= 0 is given by g2 q (?) ' r)=i q (a) (a)q1(r)-ql(a)q2(r)J (6. 104) when r<:, while i2 q (r) g, r)=4i q(a) [q2 (a)q1(5)-qI (a )q2( (6.105) for r >. The complete expressions for (r) and $n (r) are obtained in n n manner analogous to the one used in treating the planar case, and it may be shown that 2 OD n(r)=- i U / G (,r)R ()d + (a)dl)(r)+ ' (a)G(2)(r) (6.106) n..2 n n nn n n e 1 a where pI(a)p (r)-pl(a)P3(r) G ( 3r)- (6.107) n pl(a)Pi(a)-p3(a)pl(a) and () P3(a)Pl (r)-p (a)p3(r) n p - (a)(a)-p (a)p(a) (6.108) Likewise, 2 h(l)(ar) n h(1)(a n a n 99

From (6.22), (6. 90) and (6. 92) we note that -22 R()=erw i [6(ww-6(W-w 4ir ri-( 2+- ) ) u. Ue ui r 20r 4w[ ( 1 d 2 d) ndn+l)6 iu 12 d d 2 J 2 J (6. 110) so that. 2 0? WC 2 2 n n U 2 Gn( ~r)Rn( )d^= u J e 1 a ir ieruw O(w O )- 6(J-W) ( ) [a) (a)N(a)-p3(a)pq(a)] "2 2 2 p(o1(r)+ 2+ 2 3- (ro)3(r u u. W u W - e i e 1 c (a, r )X(r) +c(a, r )3(r) Similarly, 2 4 2 56(-r) C Z (j)=-.. ---enw.o [(u+)-6(CO-.o)l 0 w n i o 0 2 whereby (6.111) (6. 112) 2 Cf g(h, r)Zn( ) = e rxo [6 (wo ) - 6 (W ) ii n o 0 o 0 q (r ) 1 oa q la) q2(ah)q(r)-q(a)q2(r] = Kl(a, r()q(r)-K2a, r )O(r) 2 1 0 q 2 0 2 (r (6. 113) 100

Using the definition o0 p(r, ))= 2n+ P (cosO ) (r), (6. 114) 2 n n n=0 we have, from (6.106) and (6. 111) the result 0O 0(r, e)=. 2n+ P(cos) ^lnn(r)+.(r)+ 0(a)G(1)(r)+0'(a)G(2)(r)l (6.115) 2 n I" n nn n ' n=0 Similarly, we obtain from (6. 109) and (6. 113) 00 (1) 2n+1 n n h (r) r, (cos)) (r)+Kq2(r)+ ( (a)} (6.116) 2 2(r ) h(1)(aa) n n=2 n Equations (6.115) and (6.116) represent a formal solution. The problem remains, however, to determine the arbitrary constants (a), 0'(a) and Pn(a). This may be done by considering the boundary conditions stated in (6.46) through (6.48). Substituting in the functional expression for E from (6.23), these requirements reduce to 2 ~i _+ 2 cos0A = 0 at r = a (6.117) car 2z c 2 2 V -(1+)) [1 - os A =0 at r=a (6.118) ar c 2 2]o A U W u e 1 a0 iw 1 - -- a A:0 at r=a (6.119) sin/ 0a c z 101

The term which causes the most difficulty is V 2l. After some algebraic ar ' manipulation, we find from (6.107), (6.108), (6. 111) and (6.115) that 0O a2 1 2n+1P (cos9) -— 3 2 p'(a)+c3p;(a] a 0 ra 2 -2 )p '(a)p (a) 2 1 - 2P3 P) + IP3 (a) r=a (a) (Pi pi- p3p 1 () n 1 P3r- P3P1 n r=a r=a - 2nl P (cos)F (ar)+b b (a) - d 0(a (6.120) 2ir 2 nn nn where F (a, r ), b (a, r ) and c (a, r ) are known quantities. n o n o n o Turning nowtothe boundary conditions and taking (6.119) first, we note that if x = cos0, (2-), n odd n ( -1), n even ap (cos0) aP (x) ( 1 n n n-(4m+l P (6.121) sinO O ax - ( n-(2m+1) m=O 102

Thus OD - 1 al~.,O _ 1A 4n5 0(a) sin9 d O 2 Lr 2n+ n=O n z[n-(4m -3)]P m=O 2n-(2m-1) + n=O 4n+3 (a 2 2n+l ~ m~O [4-(4m -1)] 2nm =..1 4~ zi=O,U= 50 (a) 22n+2 n OD n - (4S+3)P 2n+24n3 () 4S+i)~(x)r = S0 2S+l I2 (anzZ( = n=O S=O (4S+3)P (x) 2S+1 OD... (a) 2 2n+2 n=S OD S0O 0D 'S'4n+3 - (4s+1)P (x)Z~ 2 (a) 2S 2 2n+1 n=S (6. 122) Thus, from (6.119) together with (6.116) and (6.122), >T {4s+3)P (x S0O 2S+ In=S OD 4n+5 (a) +(45+1)P x 2 2n+2 2S n=S 4n+3 (a 2 2n+l + 00 s=O (6. 123) 103

From the orthogonality of the Legendre polynomials, we obtain the first recursion expression, namely a qJ (a)+ C m i) aSm(a)+ o (4n+5)0 (a) O 0 m-12n+2 m-l n= 2 Z (4n+3) (a) 0 2n+l n m/2 m odd (6. 124) m even To show that the two sines expansions do indeed converge, we note that 0o and p(a, 0)+ p(a,Jr) = Z (4n+1) (a) 2 n 2n n=0 oo (a, ) 0p(a, ) (4n+3) 0 (a) nv 2n+1 n=O (6. 125a) (6.125b) Thus, in place of (6.124) we could have written (a)-= - 2[(a., )-P(A.,) (6. 126a) 104

(m-1) iwa {2)0(a,O)-O(a,)]- > (4n+3)p (a) (6.126b) m na 2n+l n=0 for m even, m/ 0 m-1 p (a)= - c [ 2 (a, 0)+(0a,r)] - (4n+l)p (a) (6.126c) m ia 2n na0 for m odd The condition in (6.117) is handled in much the same manner, although the situation is simpler. Omitting the details, we find that 2 W (6. 127a) w o(a) + — 2 l(a)^ 0 (6.127a) C o 2 c 2 u tO c (a) + 2 1 [mi (a) +(m+l)P (a)], m> 0 (6.127b) c m-1 m+I Using the expressions from (6.126) we thus obtain 2 0(a) - - - 2- 2' [,(a, 0)+ (aa, a) (6. 128a) 2 a J 2.(a) — 2 2[r (a, O)- (a.r)] -26 (a) (6.128b) 2w a 105

2 (T2 n00 2n+1 u -1 F (ar)+b- 2 (a,))-d +(a, -(m+l) (a)- (4n + 1) (a) m 20 2 Iao m weehave 22 2 2 F (a,r )+b (a)-rd+ - (1+ X)(1- 1'(a) x o (6.129a) 0m 0 m 2 2 2 0 Ua } e 2 n 2 m o m m 2 6.28 m u to Substitution of the expressions for (a) from (6. 128) into (6. 129) thus m yields a recursion scheme for O (a) in terms of known functions. Once these values are determined, we can obtain 1 (a) and ft1 (a) from (6. 126) and m m (6.128), respectively. 106

We have thus been able to establish a formal solution for 0(r) and W(r) which may be evaluated for particular cases of interest. It is not too difficult to understand that the process for any given situation is not particularly simple. Nevertheless, it should be possible to obtain solutions in the shadow region for example, by using a Watson transformation, or a near zone solution by retaining just the first few terms of the Legendre polynomial expansion. 6. 5 Discussion of Result Two points are worth noting in regard to the preceding discussion. In the first place, it is quite evident that the difficulties encountered in the analysis were due to the fact that ion motion was taken into account. This effect served to give us a fourth order equation in p rather than the second order equation we would have obtained had we considered the ions to be stationary. More devastating from the standpoint of a rigorous treatment, however, was the fact that the boundary conditions were vastly more complicated. It was illustrated in the third part of this section, for example, that even in the case of an oscillating dipole above a perfectly conducting plane, we obtain a solution which cannot be integrated in any simple manner. The second point that should be mentioned regards the matter as to just how realistic the boundary conditions really are. The assumption that the particles undergo elastic reflection when striking a perfectly conducting object has been used by practically everyone treating this type of problem. The fact remains, however, that an electron striking a metallic surface might either be absorbed or undergo some diffuse reflection. An accurate description of the actual phenomenon has never, to our knowledge, been described in a satisfactory manner. This is one facet of the problem, then, that definitely bears further investigation. 107

VII INHOMOGENEOUS MEDIA 7.1 Introduction Although the propagation of waves in an infinite homogeneous ionized medium has been extensively investigated in the literature, the excitation of waves and their subsequent propagation in an inhomogeneous media, such as the ionsphere, has been scarcely touched, due to the complicated nature of the mathematics involved. In this section we shall discuss possible methods, approximate or exact, for treating such problems. If the local dispersion relation (and propagation constants) are known as a function of the spacial variables, the excitation problem may be solved approximately by considering the fields excited in the vicinity of a source and extending these fields to other points in space by ray tracing techniques. A general discussion of ray tracing in ionized media is given in Section 7.2. For drastic changes in the properties of the medium, the ray tracing technique may not give accurate results due to the neglect of intermode coupling. For these regions the reflection and refraction of the waves should be used. A general discussion of the reflection and refraction of waves, due to discontinuities in the medium, is discussed in Section 7.3. For a general variation of the properties of the medium in one direction (i. e., stratified), the recently developed technique of invariant embedding, probably would also apply. The application of this technique to wave propagation in ionized media is discussed in Section 7.4 108

Our primary goal in this work is to develop a unified approach to deal with excitation and propagation problems in both homogeneous and inhomogeneous media. The generalized operator transform technique developed during the course of this investigation is designed to accomplish this goal. The application of this technique to excitation problems in homogeneous media has been demonstrated but due to lack of time and the complexity of the mathematics involved, the application of this method to problems involving inhomogeneous media remains to be investigated. In Section 7. 5 a general discusslon of the application of the operator transform method to problems involving iihomogeneous media is given. 7. 2 Ray Tracing If the propagation constants for different "modes" of propagation are known as a function of space, the propagation of disturbances initiated at any point may be approximately calculated by the well known technique of ray tracing. The application of ray tracing techniques to anisotropic (propagation constant depends on direction) and dissipative media has been discussed in detail by Brandstatter (1963). If the propagation constant does not change appreciably in a wave length, and the inter-mode coupling is negligible, then knowing the intital amplitude and direction of a ray at one point, the ray path and amplitude can be obtained by numerical integration. In terms of the propagation constant written in the form (assume real, for the time being) s= n(r, a) cwhere a is a unit vector indicating the direction of propagation, the standard form of the ray equation may be easily shown to be 109

dr 1 G( ) (7.1) dr 2 3a n - d_ 1 3G(,) (7.2) dT 2 r r n - where a ^ n(r a) (7.3) G(r, ')=2 -a-n (r_,)] (7.4) and T is the time of propagation. Equations (7. 1) and (7.2) are a set of six (6) coupled equations. Knowing the initial position r, and direction r, these equations may be integrated numerically in steps to yield the ray path in the parametric form r (r). The direction of phase propagation at any point along the ray path is also obtained in the parametric form a_ = _(r). For a stratified medium such as the ionosphere, the ray equations can be simplified further. Following our discussion on the propagation constants for each mode, we have a value of n in the form n(r, p) where P =cos0 is the cosine of the angle between the d. c. magnetic field and the direction of propagation. In order to adapt these results to the ray equation, we shall choose a new coordinate system such that the medium is stratified in the z direction, while the d.c. magnetic field is in the direction ^A b=zcos0 +ysinO. (7.5) o o 110

If the direction of phase propagation is designated by A A A a z cosa+x sin a cos 3+ y sin a sin P then cos a cos 0 + sin a sin B sin (7. 6) 0 0 Introducing these equations into (7.1) and (7.2), and simplifying, one finds that = constant (7.7) dc_ 1 an(z, i) 2 az (7. 8) dT 2 3z n dx n sinacos 3+Psin cos P an (7.9) dr 2 J n dr =2 "nsina si ina+(u s inan-sin)- ) (7.10) and dz 1 I an dcr 2 ncost+(Pcosaos- cos0 ) (7.11) With some initial values of a and 0 these equations can be integrated in a straightforward manner. If several rays are traced so that the variation of the cross section of a tube of rays can be calculated, then, by conservation of energy, in the lossless medium, the amplitude of each field component can be 111

computed from the fact that they are inversely proportional to the cross sectional area of the ray tube. When collisions are not neglected, the propagation constant is, in general, complex. We may write the real and imaginary components of s in the form s= [n(z,,)ix(z,1)4 (7.12) Then, approximately, the ray path is obtained from the real part of s, i.e., n(z, P. The fields, however, are now attenuated. The attenuation along a ray may be approximately given by St dr exp a- d(z,)a * -ds where the integration is along a ray path. In general, the ray tracing technique is valid when n does not change drastically, such as near a discontinuity. Near the discontinuity, intermode coupling necessarily exists, and the investigation of such reflection phenomena shall be treated in Section 7.3. The ray tracing technique is also not valid near the regions where n is near zero or infinity and changes sign. Near the infinity of n, the wave is totally absorbed, while near a zero of n, the wave is reflected (Budden, 1961). 7.3 Reflection and Refraction The technique of ray tracing, applicable when the properties of the medium do not vary appreciably in a distance of one wavelength, neglects the coupling of the various modes due to the inhomogeneity of the medium. For strong changes in the properties of the medium, especially for large gradients of temperature 112

or density, such that the medium may be idealized by means of a discontinuity, the intermode coupling may be investigated by the conventional methods of reflection of waves. Assuming that the surface of the discontinuity can be represented by a plane, the calculation of the reflection and refraction of plane waves can be carried out, in principle, by assuming the forms of reflected and refracted waves at the discontinuity, and determining the amplitudes of these waves in terms of the incident wave by use of the boundary conditions. Although such methods are standard for acoustic waves and electromagnetic waves (in dielectric medium), their application to an ionized medium becomes very involved doe to the existence of several modes simultaneously. Existing investigations have been carried out only for the case of a plasma without a d. c. magnetic field, such as reported by Kritz and Mintzer (1960), Haynes and Kahn (1965), etc. In this section, the general formulation of the reflection and refraction problem in the presence of d. c. magnetic field is given. However, due to the complicated nature of the expressions involved, no attempt is made to derive explicit formulas for the reflection coefficients. Consider the discontinuity of a medium to be in the plane z= 0. For clarity, we shall. use the superscript (+) to denote all quantities in the region z > O., and superscript (-) to denote all quantities in the region z< 0. When a wave is incident on the boundary, say from the region z< 0, waves would be partially reflected in the region z< 0, and partially transmitted into the regions z > 0. The amplitude of the reflected and transmitted waves are determuned by the boundary conditions at z=O. Assuming that the velocity components associated with the waves are small (in conformity with the perturbation approximation), so that the effect of the distortion of the boundary is negligible, then the following set of boundary conditions is sufficient to determine 113

the reflected and refracted amplitudes: a-The electromagnetic conditions, E =E x x E =E Y Y + - h =h x x (7.13) (7.14) (7.15) (7.16) h = y h y b-The kinematic conditions, V =V ez ez + V. = V. 1Z 1Z V = V nz nz c-The dynamic conditions, (U +2n = (U -2n e e e e (7.17) (7.18) (7.19) (7.20) 114

+2+ 2 - (Ui )- = (U) n. (7.21) (- + -2 (U n = (U )n. (7.22) n n n n Explicitly, let us consider the case where electron motion is dominant and collisions are negligible. Then, assuming that the magnetic field is in the y - z plane, such that A A b = zcose +ysine, (7.23) 0 0 the propagation constant in any direction A A A A s = z cosa + x sina cos y+ z sina sin ' may be determined through the dispersion relation 2 2_ 2 222 2 2 22 2 2 2 2 o2-[e2( W [- 2(12-w2)] -22 [2 2] s2 [-C 1 -s2 (-2)] 2 +fsi20 [s2-02][s2-o2(1-U2)]82 = 0 (7.24) o o p o e where cos 8 = cos 0 cos a + sin sin a sin y (7.25) 0 0 By algebraic manipulation, the field quantities associated with each mode, i.e. each s satisfying Eq. (7.25), are determined within one constant. 115

This means that for each s, we have E E h h V n x y x y ez e -. = = - = = = - fl f f f f5 f6 where the functions f. are given by 1 (7.26) f (a, ^ys,,3 ) 1 p e sinO cosa-cosO sinasinu'v 2 2 2 0o 0o _ 2- 2 (2 -2.[42 (I-w )2 cose e p 0 p o 0 -isirncos(s2- ) cosO [s22(1 -20 2+sine [s2-2(1 -W)] S 2J (7.27) 222 2 2 22 2 c~Sos i2sinroiy 2 9- 99 +i~(s7- 2) (sineo -sirsin-ycos6) [s - (1- )] 3 -sinrsinsin [s ((7. p 28) (7.28) f3(a,, s, O, 3 p 0 S22 2r2 2 2, = - -iQ(sinO cosa-cosO sinarsin^/)[s-o -0 (1 )1 e w o o 0 0 e p o - sinacoss -e (1-j) - (I - )~ P e p 0 p I (7.29) 116

2222 S 1 2 0 2 2 f4(a, s,, ) -iS cosO sina sins2-1, ) [s -3 (1-w ) P o p sine sinasin-y + — - 2 2 2 1 - r2 (1 2)]J (7.30) 5' p -2 (s-13)ocos e o p e p p 22 22 -il(cose -cosacose)(s2-2)[s 2 (1-w2 )] -ifcosasin8 -/ (1-w )] o o e p p e (7.31) WE f(0, 'v s, 0,~ )= - ns sine222 2 2- 2 fvs,,3 )=- ssino s -0 [s (1 - (7. 32) 6 p -e e o o p e where we have explicitly written down the arguments of the f's to indicate their dependence on the properties of the medium. For any mode corresponding to a propagation constant s., the associated field E, E, h, h, V and 1 x y x y ez n can then be represented by a column matrix e is.(cosa z+ si r.cos.x+ sina.sinY.y) A[i f(ai*,Y i si' p e)]e (7. 33) where A. is the "amplitude" of the wave, and 1 117

f f f [f] 3 (7.34) f 4 f5 are the appropriate components. Across any discontinuity, the total field must be continuous. If a wave is incident from the negative region on the boundary, and the incident field is given by mode i, corresponding to s., o, vi 10 io 10 then we have the following reflected and refracted fields. Incident field: is. (cos. z + sinc. cos '. x+ sinc. sinv. y) afr " ~ ~ o~ 1 1 1 e 1 1 ' 1 0 io p e Reflected fields: corresponding to j = 1, 2, 3: is.(cosa.z + sina.cos'. x+ sina.sinv.y) Rij [f, ^v., s.,I W 3e 1j j J J P e and the transmitted (refracted) fields, corresponding to j= 1, 2, 3, + + + + + + + + + + is. (cos a. z+sin a. cosN. x+ sina. cos'. y) +K + + +] ] ] pe TY, s.,c,1P e ii j~ j p e' 118

The continuity of fields at z=0, therefore, yields the following relation on the direction of the reflected and refracted waves - - - - - - + + s. sina.cos. = s. sina. cos. = s. sina.cos Y. 10 1 1 J J j J J J - - - - + + + s. sina. sinY. = s.sina.siny. = s.sina.cosv. io0 1 i 3 J ] ] ] (7.35) (7.36) where s., a., v. are fixed by the incident wave and j=l, 2, 3, corre10 10 10 sponding to three reflected and refracted waves. From Eqs. (7. 35) and (7. 36), it is evident that y is the same for all waves, while the a's are related by - - - + + s. sina. = s. sina. =.. sina. 10 10 ] ] ] ] (7.37) where on physical grounds, we restrict + ir a. < -2 and - > a 2 For the particular j= i, we have immediately a. = -a. 1 10 ' (7.38) (7.39) S. = S. 1 10 119

for other values of j, s., a. have to be determined by the dispersion relation. After the values of s, and a are determined, the following vector equation enables us to determine the coefficients of reflection and refraction 3 [f(a i,. is. p 13 R [f(.,..,,0)1 lo lo lo p e ij j lo j p e j=l 3 =T f(a,^ -v..1 SJ W.^+ (7. 40) j j lo J p e j=l Since there are six component equations in Eq. (7. 40), the six coefficients R..ij, Ti.. may be solved by inverting a 6x6 matrix. This can be done by a computer with relative ease. Extension of this formulation to a collisionless two-fluid plasma is straightforward. In this case, the dispersion relation has four roots, and the column matrix has eight components. The general form of the result should be the same as that discussed above. 7.4 Invariant Embedding The effect of the intermode coupling on waves propagating in a region where the properties of the medium are changing rapidly, can be investigated by the familiar principle of invariant embedding. This technique has found wide application in the fields of light and neutron scattering and can be applied to the present problem. In this method the region. where the properties of the 120

medium are changing rapidly is considered as a layer rather than as an idealized discontinuity as discussed in Section 7.2. The mathematical formulation of "invariant embedding" has been summarized in an article by Baily and Wing (1963). In the section, the problem of evaluating the reflection and transmission coefficients of a layer of stratified, inhomogeneous electron plasma will be formulated in a manner suitable for numerical computation. The formulation may be extended, in a straightforward manner, to include the case of ion and neutral particle motion, but, of course, the numerical computations are much more involved. For some background on the principle of invariant embedding, consider the set of coupled differential equations given by dz = A(z) U(z)+ B (z) W (z) (7.41) dW(z) - d C(z) U(z)+ D(z) W (z) where U(z) and W(z) are nxl column matrices, while A, B, C and D are nxn matrices. If we interpret U(z) as the field quantities associated with the forward propagating wave, and W(z) as the field quantities associated with the backward propagating wave, then the reflection and transmission characteristics of a layer in a region 0 < z < z may be investigated by solving Eq. (7.41), subject to the following boundary conditions U(0) = U (7.42) W(z ) 0 0 121

where U is the vector, whose components are the amplitude of the 0 incident wave. The overall effect on wave propagation, due to the presence of this layer, can then be expressed in terms of a reflection matrix R and a transmission matrix T. In terms of these matrices, we may evaluate the reflected wave at z = 0 and transmitted wave at z = z by the relation 0 U(z ) = TU 0 0 (7.43) W(o) =RU 0 By using the principle of invariant embedding, it can easily be shown that, in order to obtain T and R, it is not necessary to solve the linear boundary value problem (Eq. (7.41) and Eq. (7.42)). In fact, T(z) and R(z) satisfy the first order nonlinear differential equations given by the following -R(z) = R(z)B(z)R(z) + R(z)A(z) + D(z)R(z) + C(z) (7.44) dz d T(z) = T(z) B(z) R(z)+ T(z) A(z) (7.45) dz This set of nonlinear Eqs. (7.44) and (7.45), can be integrated from z=z 0 to z= 0 to obtain the reflection and transmission matrices. The obvious boundary conditions, for the integration, are R=0 and T=I at z=z 0 To apply the above formulation in the calculation of transmission and reflection matrices for an inhomogeneous slab, let us consider an electron plasma with non-uniform properties in the region defined by z > z> 0. 0- - 122

The regions for z< 0, and z > z, are assumed to be filled with electron 0 plasmas of uniform properties (or free space which can be considered as an electron plasma with zero electron density). Let waves be incident from the region z< 0. From the discussion of the reflection and refraction of waves in Section 7. 3, we note that the propagation constant in the x direction (s sin a cos 'v) and in the y direction ( s cosa sinv) are constants. Therefore, each field component, such as Ex, may be represented in the form is IC is y E -E (z)e e Y x x where the amplitudes of the field components are functions of z only. The source free equations for the collisionless plasma, where the direction of the d. c. magnetic field is the same as given in the Section 7.3, can then be reduced to a set of coupled differential equations. The explicit forms of these equations are: 123

dh x dZ dh dZ dE dZI 0 s s h~I0 2 0 0 ss 0 32 0 wm )28 e p x iujm Li2 e p y 0 0 0 h x hy E x 8 i2O s 2 U2 Faa iw 3coa91- a w 1 {(32 1-02Co2 089 -w I 2 -2Co2aj L 01-2co0 w~m w 2 2sin~cosO 0 e 1-1 22Cos2 8) 0 dE -x dZ dV z dZ dni dZ 2 2 -if 0 - 22- - -2j 0 e (a +ia 1fZcoa6) WMr( I -dCOS29 ) iw2i2 1 QainO p e e (1-fl Cos2 9) 2 [ iua i C.Zeos i U)E.x + P 0 (a - is 1lcos9) e y x WM 1-.1?cos0o c 2 2 dsinftosO e p e 1-12Cos 2 0 2 2 iw ME U) 1Iin8 e( -12 sin2 9) 3 (12 G 2Co2 e (a - Ocosesa 2 1-2 Co2 e,ie2(a2 +2 x ~ IWME Li 2[ 2(l _f22cos29) op e Qsin9(s -jfkosesa (1-02 Cos2 9) 0 E Y 0 aw 2 w2 az 2- 2Co2 (i12os+ p op e _ _ __ L2> e2 p2 ', w2 2 we y 2a2 we xa V z n.1 (7. 4(5)

All the field components in Eq. (7.46) are continuous across any boundary. To reduce the above equations to a form that is adaptable for use in the principle of invariant embedding, let us assume that for the region z<O, the propagation constants of the three modes have been calculated and denote their values by sa, sb and s, respectively. For any fixed direction of the incident wave, s and s are fixed. If we denote x y 2 2 U =+ s -s -s,(7.47) a a x y and similar expressions for Ub and U, then the amplitudes of the incident b c waves may be written as iU iz iz i z P f(U )]e aPf( + P [f(U)) e P U e c (7.48) a a c c where the f's are the column vectors for the field quantities. The explicit forms of the f's are given in Eqs. (7. 26) through (7. 32), and a simple algebraic substitution will yield the expression for f in terms of the U's, instead of s, a and. Similarly, the reflected waves are represented by -iU z -iU z -iU z [f(-U e a + b e b + Q [f-U e c (7.49) The "amplitudes" of the reflected wave can, then, be expressed in terms of the reflection matrix, such that Qa a Q Pb (7.50) Qc Pc 125

Note that [R] in Eq. (7.50) is not the same [R] given in Eq. (7.43). Similarly, for the region z > z, we may denote the values of the propagation constants as s, sb and s, respectively. From the values of s and s, fixed a c x y by the incident field, we can calculate ^: /32 2 2 U = s -s -s. (7.51) a a x y The transmitted fields in the region z > z are, therefore, given by the O general form -i (z-z -i%(z )Z -i (Z-Z) f(U)+ Je e + Qc [+Uc)] e c (7.52) -" where the f's are the appropriate functions for the field components corre - sponding to the properties of the medium in the region z > z. In terms of a O transmission matrix, the "amplitudes" of the transmitted waves can be calculated from 4 p a a Q% = (zo) Pb (7.53) dQ P c c The transmission and reflection matrices may be calculated by the principle of invariant embedding, if we can reformulate the set of Eqs. (7.46) 126

in terms of forward and backward propagating waves. Utilizing the property that the field components must be continuous, we see that the boundary conditions for the differential Eqs. (7.46) are (i) at Z = 0, E x E y h x h y V z Pa(Ua) + Pbf(] + Pf(Uc +af(-Ua +Qbf(-U +Qcf(Uc)] (7.54) n and (ii) at Z = z, o E x E y h x h y V z n - - Qaw ( a)-] Y r^Ub)] Th ^] Au Qu) +.0 UJ+Q f( ) a a- Qb f b- c cl (7.55) We may eliminate T and R from Eqs. (7.50), (7.53), (7.54) and (7.55), by forming linear combinations of E, E, h, h, V and n that may be x y x y z 127

identified as the forward and backward propagating waves. Such algebraic derivations are somewhat complicated, but are, nevertheless straightforward. The procedure is stated here: a) Form the determinant fi(U) f (U) f (U) f (U) f(u) f (U) a 2a a a a fl (U) f(Ub) f3(Ub) f4(Ub) f5(Ub) f (Ub) fl(U) f() f(U) f4() f5(UC) f6(uc) (7. 56) fl(-U) f2(-Ua) f3(-U ) f(-Ua) f5(-Ua) f6(-Ua) fl(-Ub) f2(-Ub) f3(-Ub) (-Ub) f5(-Ub) f6(-Ub) fl(-U) f (-U) f (-U ) f4(-U) f5(-U ) f(-U) 1 c 2 c 3 c 4 c 5 c 6 c and obtain the cofactors corresponding to each of the elements of the first three rows and denote these cofactors by a., b., ci (i = 1, 2, 3, 4, 5, 6). b) Form the determinant 128

f (-u)f (U-) (.U- f() f(-u) f (1 a 2 a 3 a f4 aU 5 a 6U r f(U) f (-U ) F(- f (-U) f (-U) f (-UP 1fc 2 c) 3 C 4 C 5 C 6 c~ (7.57) fr(ii) i) flu T(u f(U) f u) Ila f2Ua 3 a f4a 5 a 6 a f (1) f0) f3 ) f( ) f( ) f( ) 2l3 4 b S 6o I c 2 c, 3' C 4 c 5 c 6 c and obtain -the cofactors of the elements in thie first three rows. Denote these cofactors by a, b, c (,i = 1,9 2, 3, 4, 5, 6) c) DEk.fine LUJ PXy by the transformation,, 1 2 3 a4 5 6 E b b b b b bE I 2 3 4 5 6 y 1 C5 C6 h ~ 6 hy(7. 58) b b b b b b 1 2 3 4 5 6 z L "2 ~3 4 CsC6 129

d) It is easily seen that, at z = 0, P a [U a b A (7.59) P C and A A P Q la lb lc a a wJ AP + A Q (7.60) W = a 2b 2c b b (7.60) 3a 3b 3c c Qc where, la is the determinate formed by replacing the first row of Eq. (7.57) by f(U, f(Ua) f(U ), f(U ), f (U), f (U). The other l a 2 a 3 a 4 a 5 a 6 a determinants are similarly defined. The boundary conditions at z = z are O [] = 0 (7.61) and la Alb A1 Qa la lb Aic a [U] 2a A2b A2 (7.62) 3a A3b 3c Qc where Ala is the determinant formed by replacing the first row of A by fl(U), f(U ).. f(U). The other determinants are similarly defined. 1 a 2 a 6 a 130

e) From the transformation given by Eq. (7. 58), the differential equation for U and W can be obtained. This set of equations is in the standard form of Eq. (7.41), with the same type of boundary conditions as Eq. (7.42). Thus, Eqs. (7.44) and (7.45) may be integrated to obtain [R] and [T] appropriate for [U] and [w. f) For the reflection and transmission of the actual fields, we see that from Eq. (7.62) Ala Alb Alc A2a A2b 4c A3a 3b 3c -__ [F] P a Pb c P a = [T Pb P c (7.63) m Thus, T[] Ala A2a 3a Alb Aic A2b Ac 3b 3c -1 [T] (7.64) Similarly, from Eq. (7. 60), Ala A3a 3a Alb 2b 3b 0 —0 Alc 3a p a Cb pC +A [R P a b P c P a = [T] Pb P c (7.65) 131

Thus, la lb Ic [R] [T]- 2 2b (7.66) A A A2b %0 3a 3b 3c The procedure outlined above seems to be complicated, but numerical results can be obtained in a fairly straightforward manner. 7. 5 Operator Transform Method A general and, perhaps, more basic approach to the investigation of the excitation and propagation problem in an inhomogeneous, ionized medium is probably the generalized operator method outlined in Section IV. It is shown there, that for the most general excitation problem, this formulation reduces to the solution of an integral equation of Fredholm type, i. e., #1(u) = F(u)+ $ k(u, s) (s) (7.67) where the kernel k(u, s) is a 6x6 matrix. For stratified media, this equation involves only a one-dimensional integral, so that a numerical solution of such an equation may be within the reach of present day computer capacity. Additional simplifications obtained by taking advantage of the symmetry of the sources may further simplify the problem. In this section, some possible approximate solutions obtained by utilizing this formulation are discussed. a) The Perturbation Method. If the excitation problem in a medium is solved, for example, for the case of plane waves in a stratified media, then the effect 132

of a slight disturbance of the medium on the resulting wave can be investigated approximately by perturbation methods. Mathematically, this means that if the solution i (u) of 0 O(u) = F(u) + $ k (u, s)7 (s) (7. 68) is known, then for a slight perturbation of the medium, the kernel k(u, s) may be written as k(u, s) = k(u, s)+ Ak(u, s) Assuming a solution of the form (u) =O(u) + A(u) and using the familiar Born approximation Aj(u) is given approximately by ai(u) = $ Ak(u, s) I(s). (7.69) b) Series Solutions. If we rewrite Eq. (7. 76) in the form (u) = F(u) + Xk(u, s)(s) (7. 70) by either considering Eq. (7. 67) as a special case of (7. 70) for X = 1, or by scaling k so that X is a small parameter, then a series approximation for 133

obtaining the resolvent of the integral equation may apply. Formally, we may write the solution of Eq. (7. 70) as, I(s) = F(u) + X $ H(u, s, X) F(u) (7. 71) The formal procedure of obtaining the resolvent H(u, s, X) has been discussed by Diament (1963). Following the derivations for a one-dimensional integral equation, it is suggested that if we write H (us) X)= c(u, S X) (7.72) p(X) where p(X) and c(u, s, X) are entire functions of X, and expand p(X) as oD p E)AP. n (7. 73) n=O and c(u,s,X)as c(u,s,X) = c (us)\n (7.74) n=O then a numerical scheme may be employed to obtain the various coefficients of Eqs. (7. 73) and (7. 74) by the following recurrence relations, 134

p 1, c (u, s) = k(u, s) O -npn= $ Tr cnl(s, s) c (u, s) = k(u, s)pn+ $ k(u, v) cn_1(v, s) where Tr c denotes the trace of the matrix c. The rate of convergence of this scheme, however, has not been investigated. Other possible approaches for obtaining approximate, numerical solutions of integral equations, such as the variational technique, Gelerkin's method, or the approximation of the equation by a linear algebraic equation, which have been used in solving one-dimensional equations, in principle, can be applied to Eq. (7. 67). However, unless the equation is further reduced, the labor of computation is, perhaps, prohibitive. In order to take full advantage of this formulation, perhaps, the most realistic approach is to obtain complete solutions for some ideal kernels, which closely approximate actual physical problems, where the solution may be carried out, in part, analytically. Solutions of a variety of realistic problems can then be treated as perturbations of the solutions of these standard problems. Due to limitations of time, no computations were made for inhomogeneous media. 135

VIII CONCLUSIONS The primary goal in this work was to investigate the excitation and propagation of wave like disturbance in the ionosphere. As a necessary step in carrying out such an investigation, a realistic three-fluid model of the ionosphere was developed. The fundamental modes, which can propagate in an infinite homogeneous media, have been investigated and a computer program developed to compute the propagation constant corresponding to parameters appropriate to the ionosphere. Various approaches to the problem of the excitation and propagation of disturbances in an inhomogeneous medium have been investigated. As a result of this investigation, it is suggested that there are two realistic approaches to this problem. 1. Solution of the excitation problem, if the medium is locally homogeneous near the source. From the known field excited near the source, the fields in other regions may be obtained by using the appropriate formulation for a stratified medium. 2. Solution of a set of standard problems, using the operator transform method, by idealizing the kernel of the integral equation so that the complete solution of a set of standard problems can be carried out, at least in part, by analytical methods. The solution for a real problem can then be obtained by perturbation methods from the known solution of an appropriate standard problem. 136

REFERENCES Allis, W. P., S. J. Buchbaum and A. Bers (1963), Waves in Anisotropic Plasmas, (MIT Press, Cambridge). Arden, B.W. (1963), An Introduction to Digital Computing, (AddisonWesley Publishing Co.). Bailey, P.B. and G.M. Wing (1965), "Some Recent Developments in Invariant Imbedding with Applications, " J. Math. Phys., 6, No. 3, pp. 453-462. Chapman, S. and T.G. Cowling (1939), The Mathematical Theory of Non-Uniform Gases, (Cambridge University Press, Cambridge). Cohen, M.H. (1962), "Radiation in a Plasma-Ill: Metal Boundaries," Phys. Review, 126, No. 2, pp. 338-404. Cowling, T. G. (1945), "The Electrical Conductivity of an Ionized Gas in a Magnetic Field, with Applications to the Solar Atmosphere and the Ionosphere," Proc. Royal Soc., A 183,, 453. Diament, P. (1963), "A Formal Solution to Maxwell's Equations for General Linear Media, " Columbia University Scientific Report 78. Ginzburg, V. L. (1961), "Propagation of Electromagnetic Waves in Plasma," translated from the Russian by Royer and Roger, (Gordon and Breach, New York). Haynes, K.A. and D. Kahn (1965), "Wave Propagation Across a Two-Fluid Plasma Density Discontinuity," The Phys. of Fluids, 8 No. 9, pp. 1681-1688. Kritz, A.H. and D. Mintzer (1960), "Propagation of Plasma Waves Across a Density-Discontinuity," Phys. Review, 117, No. 2, pp. 382-386. Marshall, W. (1957), "Kinetic Theory of an Ionized Gas, Part 2, " United Kingdom Atomic Energy Authority, Atomic Energy Research Establishment, Harwell, Berkshire. Nagy, A.F., L.H. Brace, G.R. Carrignon and M. Kanal (1963), "Direct Measurements Bearing on the Extent of Thermal Nonequilibrium in the Ionosphere," J. Geophys. Res., 68, No. 24, pp. 6401-6412. 137

Schelkunoff, S.A. (1955), "Conversion of Maxwell's Equations into Generalized Telegraphist's Equations, " Bell System Technical Journal, 34. Schiff, L.I. (1955), Quantum Mechanics, (McGraw-Hill, New York). Spitzer, L. (1962), Physics of Fully Ionized Gases, second edition, (Interscience, New York). Stix, T.H. (1962), The Theory of Plasma Waves, (McGraw-Hill, New York). Tanenbaum, B. S. (1962), "Wave Propagation in a Partly Ionized Gas," Raytheon Co. Interim Report No. 9, Contract AF 19(604)-5984. SECRET "U. S. Standard Atmosphere, 1962, " Superintendent of Documents, U.S. Government Printing Office, Washington, D. C. Wu, Y.K. (1965), "Unified Approach to Excitation Problems in Compressible Plasma, " The University of Michigan Radiation Laboratory Report No. 6663-1-T. 138

APPENDIX A MACROSCOPIC PROPERTIES OF THE IONOSPHERE A. 1 Discussion of the Model For the purpose of this report, a three-fluid model of the ionosphere containing an electron gas, positive ion gas and a neutral gas was developed which would contain parameters appropriate to an average ionosphere in the altitude range of 100 Km to 700 Km. Since this model was not intended to represent the real ionosphere at any particular geographic location or time, some assumptions were made which would simplify analytic procedures based on this model but which would maintain the salient features of the real ionosphere. For example, the ratio of the specific heats, v, for the neutral gas was assumed to be constant and equal to the sea level value of 1.4. Similarly, the magnitude of the earth's magnetic field was assumed to be constant and equal to 1/2 Gauss. In addition, thermal equilibrium among the gases was assumed. Both theoretical and experimental evidence is available which indicates that the electron and ion temperature, in some cases, may be greater than the neutral particle temperature and thus, for example, the electron and ion "acoustic velocities" predicted by this model may be lower than corresponding values for the real ionosphere (Nagy et al, 1963). A brief description of source of the parameters used in the model ionosphere follows. Since the acoustic velocities and the collision frequencies were computed from the parameters of the model ionosphere, these quantities will be discussed separately. A.2 Macroscopic Properties a. Neutral Gas The properties of the neutral gas were taken from the U. S. Standard Atmosphere (1962). Since the ionosphere is weakly ionized, it was 139

assumed that these parameters would provide a reasonable representation of the values for the neutral gas. Curves of the temperature, molecular weight, neutral particle density and collision frequency are shown in Fig. 1 through 4, respectively. b. Electron and Ion Gas Single ionization was assumed and thus, the number densities of the electrons and ions are equal. In addition, the molecular weight of the positive ions was assumed to be equal to that of the neutral particles. Also, thermal equilibrium was assumed and thus, the electron and ion temperatures are taken to be the same as the neutral gas temperature. The curve of electron density versus altitude is a composite curve taken from recent data appearing in the literature and is the same curve as used by Wu (1965). This curve is shown in Fig. 5. c. Acoustic Velocities The acoustic velocities for the r gas were computed from the usual equation VRT ' r r U = M r M r where y is the ratio of the specific heats, R is the gas constant with the r -1 value 8314 joules (OK )(hg-mol), T is the temperature in degrees Kelvin and M is the molecular weight. r The value of y for the neutral gas was taken as the value given by the U. S. Standard Atmosphere (1962) for the altitude between 0 and 90 Km., i.e., 1.4. The motion in the electron and ion gases is supported by electric forces involving 1 degree of freedom only, and thus, y for the charged particle gases was taken as 3 (Spitzer, 1962). The acoustic velocities in the electron, ion and neutral gas are shown in Fig. 6. 140

700 -600 -500 -E w 400 -Q I I - - 300 -< / 200 -100 -0 --- 0 100 200 300 400 500 600 700 800 900 100011001200130014001500 TEMPERATURE ~K FIG. 1: NEUTRAL GAS TEMPERATURE AS A FUNCTION OF GEOMETRIC ALTITUDE. 141

700 -600 -500 -- 400 -LJ 2 300 -200 -100 -0 - i T I II+ I 14 16 16 I 2 18 20 22 24 26 28 30 MOLECULAR WEIGHT FIG. 2: MOLECULAR WEIGHT AS A FUNCTION OF GEOMETRIC ALTITUDE. 142

700 -600- \ 500 400 Q - $ 300 <1: 200 100 O0- I 1 I I I I 1067 10 08 09 110 1011 1012 1013 NEUTRAL PARTICLE DENSITY (No. cm'1) FIG. 3: NEUTRAL PARTICLE DENSITY AS A FUNCTION OF GEOMETRIC ALTITUDE. 143

700 - Ven -/ei 600 - vin- \ 500 w400 - t300 200 - 100 -0- I 1I I I I 1 I id3 1o2 10 10~ o10 Io2 o13 I1 10 NO. OF COLLISIONS (SEC.') FIG. 4: ELECTRON-ION, ELECTRON NEUTRAL AND ION NEUTRAL COLLISION FREQUENCY AS A FUNCTION OF GEOMETRIC ALTITUDE. 144

700 - 600 - 500 --- 400 -300 -j 300 - 200 - 100 -0 - i I I I — I 106 ELECTRON DENSITY (No. cm3) FIG. 5: ELECTRON DENSITY AS A FUNCTION OF GEOMETRIC ALTITUDE. 145

700 -UN —r i VUo600 -500 -400-.. 300 -200 -100 0 -2 10 5 102 103 104 105 VELOCITY (M/SEC) FIG. 6: ACOUSTIC VELOCITY IN ELECTRON, ION AND NEUTRAL GAS AS A FUNCTION OF GEOMETRIC ALTITUDE. 146

d. Collision Frequency Since the ionosphere is weakly ionized, the value of the ion neutral collision frequency, v. was taken to be the same as the collision frequency in given in the U. S. Standard Atmosphere (1962). The electron-neutral and electron-ion collision frequency we-re computed from the following equations due to Cowling (1945). 1 -8 T - 'V 1.8x10-8T 2N en 300 n 3 -3 T -3 v.=6.1x10 (3- 2N. ei 300 where v and v are the number of collisions per second of an electron en ei with neutral particles and positive ions, respectively. T is the temperature in degrees Kelvin and N and Ni are the number density per cubic centimeter of neutrals and positive ions, respectively. Since collisions are a momentum transfer process, conservation of momentum requires that the following equation be satisfied N mvs =N mvr r rr s ss Thus, the ion-electron, neutral-electron and neutral-ion collision frequencies are given by the following equations, respectively. N m e e. = - y ie N. m. ei 1 1 147

N m e e V = — v ne N m en n n N. m. v = —v ni N m in n n The values of the ion-neutral, electron-neutral and electron-ion collision frequencies are shown in Fig. 4 and the neutral-electron, neutralion and ion-electron collision frequencies are shown in Fig. 7. e. The Electron and Ion Plasma Frequencies The electron and ion plasma frequencies were computed from the usual relation 2 2> N q ir (U) = pr m E V r o where N is the particle density per cubic meter, q is the electronic r charge in coulombs, m is the mass of the particle in K and E the dir g o electric constant of free space in farads per meter. Values of the electron and ion plasma frequency are shown in Fig. 8. f. The Electron and Ion Gyro Frequency The electron and ion gyro frequencies were computed from the equation (a0 gr mr 148

where q is the numerical value of the electronic charge, 0 is the magnitude of the earth's magnetic field and m is the mass of the particle. Since the magnitude of the earth's magnetic field was taken to be constant, the electron gyro frequency is a constant with a value of 8. 76x 10 cycles/second. The ion gyro frequency which varies with the mass of the positive ions is shown in Fig. 9 149

700-, 600\ 500 E 400 -0 -J 5 300 200e - Svneni 100- " ne 1 0" I 10 1i9 Id8 id6 165 104 NO. OF COLLISIONS (SEC:') FIG. 7: ION-ELECTRON, NEUTRAL-ELECTRON AND NEUTRAL-ION COLLISION FREQUENCY AS A FUNCTION OF GEOMETRIC ALTITUDE. 150

700' 600 -500 -v 400 300 \ \ -J 200- / / 100 1105 107 108 ANGULAR FREQUENCY (SEC') FIG. 8: ELECTRON AND ION PLASMA FREQUENCY AS A FUNCTION OF GEOMETRIC ALTITUDE. 151

700 600 500 E - L / 400 -w 300 - 100 0 - I I I I 0 100 200 300 ANGULAR FREQUENCY (SEC') FIG. 9: ION GYRO FREQUENCY AS A FUNCTION OF GEOMETRIC ALTITUDE. 152

APPENDIX B COMPUTER PROGRAM FOR NUMERICAL EVALUATION OF THE PROPAGATION CONSTANT Due to the complicated nature of the dispersion relation, Eq. (3.11), the problem of evaluating the propagation constant was programmed for solution by a digital computer. The data read in was taken from the ionospheric model in Appendix A. Calculation is then done in single precision to find the elements for two distinct three-dimensional matrices called A and B. The third dimension is to allow for each element being complex and in double precision. Matrix B has two complete columns of zeros. Each element of a determinant is of the form a+bs, where a is an element of the matrix A, and b is an element of the matrix B. When this determinant is evaluated it will give a fifth-degree equa2 tion in s To obtain the constant term of the polynomial the matrix A is evaluated as a determinant. All determinant evaluations are done using a Gauss-Jordan pivotal matrix inversion with the diagonal elements used as the pivotal elements and the value of the determinant found by taking the product of all the pivotal elements. The method is further described by Arden (1963). The computation is done in double-precision complex arithmetic and the value returned is in 2 double-precision complex form. To introduce s into the determinant, one column of the A matrix is replaced by the corresponding column of the B matrix 2n and then used as the determinant. The coefficient of the s term of the polynomial is the sum of the values of all possible determinants with n columns of A replaced by the corresponding columns of B. 153

The roots of the polynomial are found by a specialized version of Newton's method. All computation is done using double-precision complex arithmetic. Synthetic division was used to find the values of the polynomial and its first derivative. As each root is foundthe equation is reduced. When it is reduced to a quadratic, the quadratic formula is used to find the last two roots. The initial approximations used for roots, the order in which the roots were extracted and the decision in some cases to extract the inverse of the root from the inverted polynomial, were all obtained as the result of previous analysis and hand computation done by Harold Hunter. The external routines used were: DCXPNl previously described which evaluates the determinant; PROP, which does double-precision complex multiplication; DIV which does double-precision complex division; ADD, which form the double-precision complex sum of two products of a real single-precision number and a double-precision complex number; DPMOIV, which raises a doubleprecision complex number to a real single-precision power; NEWT, an almost separate program linked to the main program by program common which, as previously described, finds the roots of the polynomial; QUADwhich finds the roots of a double-precision complex quadratic equation. System routines used, but not displayed are: DPFA, DPFM, DPFDV, which perform, respectively, double-precision addition, multiplication and division; SIN and COS, which find the single-precision sine and cosine; ZERO, which zeros out a specified block of storage; MOVER, which copies a block of storage into a second block of storage; FTRAP, which sets to zero floating point underflows; DPSIN, DPCOS, DPEXP, DPELOG, DPSQRT, ATANi, which find, respectively, the double-precision sine, cosine, power of e, logarithm to the base e, square-root and single valued arctangent.

The coding is done in MAD, a compiler similar to FORTRAN. A manual for the MAD can be obtained from The University of Michigan Computing Center. 155

REF L:RTNC LS OJ *001 60OLLAN -iS *001 PRO&"',AM LOMkMON JF,R1,42,43. 0 0 Vu-Cf(]' VALUJE-S H7AD=tIH1S2O1-0q2HP40'F %UMBER 13/14+H FUR ALTITUDE=13,S5,6HU.0O0 G~iA,5 OAUANGLrL-=F6.2/41l WE=E14.3,S2,3HWI=E14.8,S2,3 * 0 0 1 HB3=E-14i.~,S2, 3HHI=El4.3,S2,3iHBEFl4.o,S2,3HfAN=E14.8/5H KIN=C1 *003 1 4.rl,S2,41,-EI=Fl4.8S2,4A<EN=E14.S,S2,4HKNI=El4.8,S2,4HiKI[E=4 *003 1.8,S2,4H1K'E=7r14. -~/4,H UE=L-14.8,S2,3HAI=E14.8*$ * 0 0 INT ERIF' L,M *0O0 vECrOR V4LULS INPUT =s7F10.2.$ *00~ INTFU'ER ~1f*O V ECrTo~, V\AL UE-S 0=0R.0Cc F T RA P..00C S fART RE~Al) FORMA'~T IN\0urT, V I qVF I, VEN VN I, VI E, VNE, UN UI, UE W~LPHWI B 0 1 H 0 1 K — Al)L) T A. 011 U E=JE l! ~,.. 0 1 VECTJ ) VALUES L)FOI=.0SdC7.0 1 ) I M L T4 IA ( RI, 42,3(4).0 1 THq~t) II 1-4O, FO; 4ALUE S OF W=3E 3.0 1 w~I~ I 0/W/v 0 o 2= * ~' 0 2 AN: Wi 02~ R2 ( I =. /r I / 'iI 0 Z [R). (1(2 IR (4) 4 2 1)... 12 (4),R 3(2)... R3 (4)) 027 K I4 =V I 1.41 *02 KFNI = V., I /W.032 K I E VI I Fi 4.0 3 KNLJ E L/ I 0 3 L)F J uU'~/ W.Q3~ uI 1O(/.0 3 n H ' VC t h S *Q3~ KAI =) 03 K I 0E = 03k 0 *Q3 L N 1) JF ClIv1,40DI I IONA L.04; v - C T)( /JILUES) P1=3.141592'7 *04 T~~,Ikh'lu, rND,FI)R VALUES IF A NV-=0.. 04 \tJ)V = JltJM +1.04 PA IN T FtC MAT H[EA D, NUM,9ALI, W,ANJG, wF, wI B 0, BI E, N,K IK CI,K FN,KNI,K I E.04~ Il M- JS I 01 (4A, ~i, I) (71.4),C F (6.4).04 I r 3 L I 5 7 3 L i 5 7 I I 3' 5 5 7 8 9 D 1 2 3 e+ 5 6 5 7

ZEKnO. A...A(7,7,4)t... B7,7, 4)) *048 ZLER).(CF...CF(6,4) ) 049 [HETA=A^J/180.*PI *050 SINT=SIN. (THETA) *051 COSr=COS. (THETA) *052 A(l, 1,1 )=H02*WE *053 A(1, 1,3)=-B02*(KEI-KNI ) 054 Al1,2,1)=802*(.-WI) *055 A(1,2,3) =B02*( KIE+KIN+K I ) *056 OEC=O)E*COST *057 A( 1, 3, 1)=-BO2*OEC*KNI *058 OIC=OI C3ST *059 OES=OE*S INT *060 UIS=f)1 S INT *061 A( 1,4, 1) =B02* IC*KNI *062 A( 1,,3)=-B02*OIC *063 A(2,, 1 )=-. *064 A(2, 1 3)=-KNI-KME-KEN *065 A( 2,2, 1) =-l. *066 A (2,2,3)=-KNI-KNE-KIN *067 A(2,3, 1)=OEC*(KJI+KNE) *068 A(2,3, 3) =-OEC *069 A(2,4,1)=-OIC*(KNI+KNE) *070 A(2,4,3)=OIC *071 A( 3,1,1 )=R02*3EC*KNI *072 A(3,2, )=-B02*OIC*KNI *073 A( 3,2, 3) =B02*0OIC *074 A(3,3,1)=802*WE *075 A ( 3,3, 3)=-B02 ( KE.I-KNI ) *076 A( 3,4, 1 )=02*( 1.-WI) *077 Q1 A(3,4, 3)='02*(KIE+KIN+KlI) *078 A(3,5, 1 )=-B02*ES*KNI *079 A(3,6, 1)=B02*OIS*KNI *080 A( 3,6, 3)=-B02*[)IS *081 A( 4, 1,1 ) =-OEC*( KNI +KNE) *082 A( 4,1,3)=OEC *083 A(4,2 1)=UIC*((NI+KNE) *084 A(4,2, 3) =-O IC *085 A( 4, 3, 1 )=-1. *086 A (4, 3, 3)=-KM I-KNE-KFN *087 A(4,4, l.)=-1. 088 A{4,4,3)=-KNI-KNE-KIN *089 A(4,5,1)=OES*(KNI+KNE) *090 A (4, 5, 3 ) =-OES *091 A(4,6, 1 ) =-OIS(KNI+KNE) *092 A(4,6, 3) =OI S *093 A(5, 3,3)=-OES *094 A( 5,5, 1) =1.-WE *095 A( 5,5, 3) =KEI +KEM *096 A( 5,6,1) =WI *097 A 5,6, 3) = -KIE *098 A (5,7, 3) =-KNE *099 A( 6,4, 3) =UIS * 100 4(6,5,1 )=WE *101 A(6,5,3)=-KEI *102 A(6,6, 1)=1.-WI *103 A( 6,6, 3) =K IE+KIN * 104 A(6, t,3)=-KNI *105 A( 7,5,3)=-KEN *106 A( 7,t6 ))=-KIN *107

Al 7,7,1)-lj *. io A l 7,97,93) zKNI +KNE ti 1#1, 3)zKEI-K%41 eic b (I 12,1 ) =- I * 11 9]I,2,3)=-KIE-KIN-KNI *11; B 1 t3 1 )=ODE C KN I. MlI 4, 1)=-OIC*K4 Bi 3,1,1 ) =-OE1C.C1*I Bt3q132I1)=-OEC*KNI* 1 ii(3 2 3)=-O I C 1 C(3, 3 3) =KE I-KN I. 31 3, 4, 3)=-K IE-KI1N-KN I.12] 83 3t 5~ 1)= OE S *KN I 1 12 ti 3,b61)=-OIS*KNI *12: bi 3,b,3)=OIS.12' t~ (5, 5, 1 ) =-1*. /BE/iE.12' btbs,65 1): - 1. /i I / Ii 1 12 [( 7,7,91) =-1./8/N1 21 S ET U P. *12~ DCA I \JV. ( 7, Dr )ET) 1 I2 STORE. (0) *1 3 T IROU(;(-H *.1 F OR T L=I1,TL. G. 5 3 L zGS ( T L) S E TU P. 13 MOVL> C t(A ( CL 1, 1).A L, 7,4), DL l1). D.DL,7, 4)'1 UL XI '44. (7,D,tDE[T)13 STf3RE.C 1) * 1 3 MUVr-.(A 411,I,)...417,7,4),( (1,1,1)....D(7,7,4) ) * 131 C.JMVE:~ * I CL L, I,1)..B(3 L, 7,4),D IL, 1,1)...DU L, 7,4)) '1 32 0UC X.I IV/. C 'IL), JET ii 3 S TO0R E * (4) * A4 Al LON TI'IUE: '14] r HIrO;H 4?, F oF T L= 1,1, TL * G. 5.1 4 L=GS I TL) *141 T HROUGcH A2,FDIR IM=TL+1,1,TM.G.5 1 4 M=GS ( TM) *14' S E TU P. * 14 MOVER. (A(L,1, 1)... A(L,7,4),DCL,1,1)...DCL,7,4) ) *14 M)V'. ( A IM, A,).( M, 7,4) MDC'.i1, 1 ) D..DMt7,4) ) 14 UCX I IJV.( 1, D, DET) 14~ s To 0.L.. * (2 ) 1 SC MflV',.(A(1,1,1)... A7,7,4),DC 1,1,1)...D(7,7,4) ).15 MOVE R * ( r IL, 1 t1 )... 8(L,7,4),DCL, 1,1)... U( CL, 7,4)).5 MOV E I * C I i '1,,)...li M, 7,4),DCM, 1,1)...0CM, 7,4) ) *15 UL X!V. 7,O, DE T)15 STORE.I 3).5 A2 CONT P4Ur A'-JV r R. I I 1,1,91 )..A( 1,7,4),DI1,1q1)...DC7,7,4) ) I)CX INV. 1,tD, OCT) 5 StO*AE.( 5) *5 rHA0U'IJ,H A 3,oFOR L =Of 1 #L.5 1 PRI'JT FJR'4AT SIHOS2fl,1H~I1,1H=UE22.16,5H +IDE22.16*$,L,CFCL~l,1)...C F L 1,4)'6 I F (LI 94) '6 A3 cnJT I AUE *6 f w T. 6 E'LW LONTIAUE16 T RANSFER TO START.6 UIMEMISIU UErI4) 16 3 4 D I 3 5 5 7 3 L 2 3 5 5 7 3 4 L 3 5 5 7 3 4 D 1 3 5 5 7 5 D L 2 3 4 5 5 7 3 4 L L 7 3 5 5

INTEGER TMTL *162 INTERNAL FUNCTION SETUP. *16( MOVER.(CAC 2,1,1)... A( 2, 7,4), DC 2, 1,1 )..0( 2,7,4)) '16 MOVR.( C (4, 1, 1)... AC 4, 7,4),0C4, 1 1 )..D( 4,7,4)) 'ti MOVER.(C ri C1,1, 1t )...B 1, 79,4),D (1,1, 1) D.(1 97,4)) ' 17 MOVER. ~(H3C5,19*1)... BC37, 74)tDC(5,19,1) *...037974)) '17 FUNCTIOJN R~FURN 17 ENO OF FUNCT ION i INTERNAL FUNCTION STORE.(X) '7 INT EGER I,9K'7 I = 4 *X 1 UPFA.(CFC I+1),CFCI+2),DETCI),DETC2),CFCI+1),CFCI+2))'1 D PF A.(CF(I1+3),CFC 1+4),DETC3 ),pDETC4),CFC 1+3),CFC 1+4)) '8 FUNCFrInN RETURN'1 INTEGO'ER GSqX18 ENI'4) OF FUiNCT ION 1 VECTOR VALUES GSC1)=1,3,596,7'1 END (iF PROGRA4'1 I 3 5 2 3 4 5

EXTERNAL FUNCTION (N,BDET).001 ENTRY TO CCXINV..002 ZERO.(DET(1)...OET(4) ).003 CET(1)zl..004 INTEGER IJKN,L.005 THROUGH SI, FOR K=1,19K.G.N *006 PROD.(DETB(KK,0),DET).007 F2.(B1,KtKK,K).008 THROUGH S2, FORJ=1,IJ.G.N.009 WHENEVER J.NE.K.010 F2.(82, K,JK,K).011 OPFDV.(B2(lB2()2(2)tBI(I),8(2),B(KJ,1),B(KJ,2)).012 CPFDV.(82(3),B2(4),Bl(1),Bl(2),B(K,J,3),B(K,J,4)).013 S2 END OF CONDITIONAL.014 CPFDV.)1(KK,1),B(K,K,2),Bl(l),B112),B(KK,1),B(KK,2)).015 CPFDV.(B(K,K,3),B(K,K,4),BI(1),Hl(2),B(KK,3 B(KK,4)) 1016 B(K,K,3)=-B(KK,3) 3).017 B(K,K,4)=-B)K,K,4) *018 THROUGH SI, FOR I=1,1,I.G.N.019 WHENEVER I.E.K,TRANSFER TOS1.020 THROUGH S3, FOR J=1,1,J.G.N.021 WHENEVER J.NE.K.022 F1.(B2,1,K,K J).023 THROUGHS5, FCR L=1,1,L.G.4.024 55 B3(L)=-82(L).025 CPFA.(8(1,J,1),B(IIJ,2),3(L) 83(2),B(I J, ),B(IJ,2) ) 026 DPFA.(B(I,J,3),B(I,J,4),B3(3),B3(4),B(IJ,3),B(IJ,4)).027 S3 END OF CONDITIONAL.028 Fl.(82, I,KKK).029 THROUGH S4,FORL=, 1,L. G.4.030 S4 B(I,KL)=-B2(L).031 SI CONTINUE.032 DIMENSION B1(4)2(4),B31 4).033 FUNCTION RETURN.034 INTERNAL FUNCTION (Z,P,Q,R,S).035 CIMENSION TA(4), TBI4), TCI4).036 INTEGER PQ,R,S,T.037 ENTRY TO Fl..038 TB( 3)=8(R, St 3) 039 TB)4)=B(RS,4) *040 TRANSFER TO CONT.041 ENTRY TO F2..042 TB(3)=-B(R,S,3).043 TB(4)=-B(R,S,4).044 CONT TB) 1)=B(R,St 1).045 TB(2)=B(R,S52).046 THROUGH SETA, FUR T=1,1,T.G.4.047 SETA TA)T)=B(PQT).048 EXECUTE PROD. )TATBTC).049 THROUGH STOC, FOR T=It1,T.G.4.050 STCC LIT)=TC(T).051 FUNCTION RETURN.052 END OF FUNCTION *054 END OF FUNCTION *054

-,EXTFRNAL FUNICT ION NEWT *.*001 REFERENC-S ON *002 PROGRAM COMMON A,R1,R2,R3.003 VECTJR VALUES UNF( 1)=l.,O.,0.,O..004 INTEGER UEG,K, J,L *005 UIMENSIO10 (C,A)((0...5).4),(T,R1,R2,R3,R4,R5,F,OEL) (4).006 VECTOR VALUES CHEK=&1H 6(S5,0E15.9)*&.007 D E G= 5)*008 WOUT. ( 41 I.009 U IV. ( 01NIE,k1, RT).010 'hRINT C(h N O THE SQUARE OFc THE ROOTS.1 PRINT FORMAT RES,4T(1)...RT(4) *012 UPMOIV.( *5,RT,RT).0 13 PRINT CJMMENT $ THE ROOTS.0 14 PRINT FORMAT RES,RT(1 )...RT(4).0 15 VECTO)R VALUES RES=$TH S20,UE22.16,S5,5H+I DE22.16.S *016 ROOT.(14 2) *01?7 U IV.( (I'NEL,R2,r\ ).01 8 PRINT CDMMENT $0 THE SQUARE OF THE ROOTS.0 19 PRIN4T FOR,\MAT -<ES,RT( 1)...RT(4) *020 V)PMI1 IV. (.5,RT1,RT ) *021 P') I T COM)%MEN T $ THE ROOTS.022 PR INT FO(<MAT RE S, T ( 1).RT (4) *023 ROOT. I R3).024 1)I V. I lNE, <43, R) T 025 PRINTi COMM~,%ENT 11C THE 'SQUARE OF THE ROOTS *026 PRINT FORliMAT RFS,RT(l)...RT(4).027 UPMOIV.(.~,RT,Rr).028 PRINT CYAMENT$ THE ROOTS.029 PR INT FJIMA T 'S, R f(1).RT (4) *030 (K= 191,K..2,1)1 V.*(AlK,0), A, A(K, 0))).031 uWIAD. ( A( 4 ), (86), OUM) *032 iUIMEN\J10IU DUM(8) *033 (K=1,1,K.G8.4, <4 (K) =UUM (K)).034 PRINT COJMMENT 1.0 THE SQUARE OF THE ROOTIS.035 PRIN4T Fn4MAT 'NES,R,4( 1..-R4(4) *036 tUPM(IV.(.5,R4,RI) *037 PRINT CO-JMMENT $ THE ROOTS.038 PRINT Fr)T~MAT IE-S,RT(1 )...RT(4).039 PR I'41 CH'AMENT $0 THE SQUARE OF THE ROOTS.041 PRIN4T FOI<MAT kES,RT(l)...RT(4) *042 uPMnIV.(.5,R[,RT) *043 PRI'4T CWMMENT $ THE ROOTS.044 uIMF4JSIO4 RT(4).045 PRINT FniRMAT RES,RT(1)..RT(4).046 FUNC II O' R ETIURN.047 1 N TC~t4 AL F UNC TIODN ROOUT. (R).048 AUA IJ (K=1,1,K.,,.4,T(K)=A(U)EG,K)).049 TIIROJJH A1,FOR K~OEG-1,-1,K.L.0.050 P4UD. (T,l,T) *05os AU DD I *, A( K, 0),1 T, T )*052 (J=1,1,J.6.4,C(K,J)=T(J)) *053

Al CONTINUE *054 (K=l 1,K.G.4,F(K)=f(K),T(K)=A(DEG,K)) *055 r HROUGH A2 FR K=DEG- 1,-1,K.L.1 *056 PROD. (T. RT T) *057 ADD. (l.,AI(K,O),l.T,T) *058 42 CONTINUE *059 UIV.( F, r,)EL ) *060 WHENEVER.ABS.DEL(1).LE..ABS.(lE-12*R(l)).AND..ABS.DEL(3).LE..ABS..061 1 (1E-1?.R(3)),TRANSFER T] OUT *061 ADI).(1.,R,-l.,OEL,R) *062 TRANSFER TO AGAIN O063 OUT (K=1,1,K.G.4,A(DEG-1,K)=A(DEGtK) ) *064 THROUGH 43,FUR K=DEG-1,-1,K.L.1 *065 (L=, 1,L.G.4,A(K-1,L)=C(KL)) *066 A3 LONTI NU E 067 DE-=UEG-1 *068 FUNCTI'IN RETURN *069 cND JOF FJNCTION *070 END OF FJICTIInN 071

EXTERNAL FUNCTION QUAD.(F3,CR).001 REFFERENCES ON.002 PROD. (RR,82) *003 DIMENSION (B2,SQ)(4).004 ADD.(l.,82,-4.,C,-SQ).005 DPMOIV. (.5,S(JSQ) 0006 Afln.( —.5,F3,-.59SQ,R).007 ADO (-.5,8,.51,SQR(4)).008 FUNCTION RETURN.009 END OF FUNCTION.010 EXT[R"14L FUNCTION PRDO!.(AOC).001 DIM'AE lis If: fl(q) *003 DP Ftl. CAC I), A42),R C1),(2),D( I),D(2)) *G C4 OPFNk * (AC A (4) Ml 3),9 R (4), Dnl 3) it (4)) D. O5 DPF-H. ( A( ), 42'),f3),3(4 ),D( 7) r)(C4) *007 0(3) =-DC3) )( 4) =-(4) DOFFA. ( ( 1) o( 2),( 3),D(4) 1,C( C (2C)3 J1' nPFA. ( rl(~) ( 6),fl( 7),O C (3),C (4) @01 1 FUN~C TI 1:j - F TURN F"JO flF FI'IJCTIflN.1 3

EXTERNAL FUNCTION ADD.(X,A,Y,B,C) *001 VECTOR VAI.UE' Z=O. *002 OPF-r,. t X,rA - ), (ZR f1 RIZ) '00 DPFM.(X,Z,A(3),A(4), IT1,IT2) *004 DPFM. (Y,Z,B (1) tB(2),RSl,RS2) *005 DPFM. (v.Z,B(3),B(4), IS1, IS2) *006 DP-A.(R11,K tlRS1,RS,2CI i)tC(2), *007 DPFA.( IT1 IT2,IS1,IS2,C(3),C(4)) *008 FUNCTION RLrURN *^.0 FND OF FUNCTION *010 EXTERNAL FUNCTION DIV. (ABC) *01 T(1 )=R (1) *003 1 (2)=R(2) (00 T ( 3 )=-B( 3 ) *00 T(4)=-(4) *006 PRln. ( A,T,N) *006 PRD[). (R, T) *0 nPF )V. ( l( 1 ),N( 2), T( ),T( 2),C (1) C (2) *00 [l) IH)V.("( ),N(4),T(ll),T(2),C(3),C(4)) *009 DIMENSIOnl (T,N) (4) *1 FIUNCTION) RETURN *01 END nF FUNCTION FXTERNAL FUNCTION DPMnIV. (N,SO,RS) *001 DIMENSION TA(4),(R,ANG,COS,SIN)(1) *002 t^(1 )=SO( 1) *003 TAI( )=SQ(2) *004 TA(C)=-SC(3).05 TA( 4 ) =-S( 4 ) 006 PROD. ( SQ, TA, TA) *007 DPSt.RT. (TA(), R) *J08 ATAN. (SQ'(3),SO(1),ANG) *009 OPFM. (,., ANG, AN ( 1 ), ANG,ANG( 1 ) *010 OPFLOG. ( R, R 011 OPFM. (N,-.,R,R( 1),R R( I) ) 012 DPEXP. (R, R) *013 npcns. ( AN, cns) 014 OPS I N. 4G,SIN) *015 OPFM. (R,R ( L),CS,COS(1),RS(l),RS(2)) )016 L)PFM. (R,R ( 1 ), SIN,SIN( 1),RS(3),RS(4) ) 017 FUNCT IO' RETURN o018 EN'D) F FUINCTION *019

APPENDIX C MATRIX ANALYSIS C. 1 Discussion of the Method In this appendix a method for the piecewise inversion of matrices of large order will be developed. The result obtained will be similar to the method of "tearing" large scale systems developed by Kron (1963) but the work presented here will begin with a system of linear algebraic equations which, it will be assumed, have been obtained as the result of the analysis of a physical system where as Kron's technique begins with a consideration of the system itself and the applicable basic physical laws (Kron, 1963, Branin, 1959). The reason for the different approach is that in some cases of interest, such as the analysis of the plasma source problem discussed in Section II, a direct approach, such as Kron's, may not be fruitful. A simple example of Kron's approach to the problem of piecewise analysis of large systems will be given in terms of electrical network theory and results of the two methods compared. The technique discussed in this section will be applied to the excitation of disturbances in a homogeneous plasma and the extension to a stratified media discussed. C. 2 General Formulation Consider the system of linear equations given by y1 allxl +a12x2 + * +ainxn (C. 1) 165

where a.. 0 The systems of Eqs. (C. 1) can be written in matrix form as Y=A X. s The matrix A of order n can be expressed as a product of matrices as follows A =PBQ s (C.2) (C. 3) where B is a diagonal matrix and p..i and qij have values of +1, -1 or 0, and the indices i and j take all values from 1 to n. From Eq. (C. 3) an element of A is given by 5 a..= p. b q.. ij Z ir rs sj (C. 4) r s Due to the fact that B is diagonal matrix, Eq. (C.4) reduces to a.j Pirb q. i] ir rr rJ (C. 5) A particularly useful form of the matrices P and Q can be obtained by arranging the first n elements of B such that b.. a. 11 ii (C. 6) 166

where a.=f..+... f../ 0 (C. 7) aii = fiii f~i (C. 7) 11 11 11 P, B and Q can now be partitioned so that Eq. (C. 3) has the form A [U G] F +L 0 s n (C. 8) 0 B K where the elements of F, and L are f.. and I.. defined by (C. 7) and o 11 11 the elements of B2 are the aij. U is a unit matrix of order n and 2 n 0 is the zero matrix. Equation (C. 8) can be written in terms of P and Q as A:F +PMQ (C.9) M= L 0 -0 B2 -At this point it must be recognized that some of the entries in M may be zero. Since, in later analysis, it will be necessary that the inverse of M, or some submatrices of M, exist it will be convenient to eliminate the zero entries at this time. This can be done easily, since M is a diagonal matrix, by partitioning M in such a fashion that all zero entries on the main diagonal of M appear in a submatrix of M. The resulting expression is as follows 167

M2 Q2 [1' 2 'm] |. | (C. 10) M Q m m Equation (C. 10) can be written as a sum as m P.M.iQi. (C.11) i=1 Performing the multiplication indicated in (C. 10) or the summation in (C. 11) the zero entries contribute nothing to the product and the following result is obtained. m GHK= PMi Qi (C.12) i=1 Where H is obtained from M by deleting the appropriate rows and columns of M corresponding to zero elements on the main diagonal. Thus, H is a diagonal matrix with no zeros appearing on the main diagonal. G and K, of course, are obtained from P and Q by deleting the columns of P and the rows of Q corresponding to the columns or rows deleted from M. Equation (C. 12) can now be substituted into Eq. (C 9) and the following equation obtained 168

A =F +GHK s o (C. 13) H in Eq. (C. 13) can now be partitioned into s submatrices, all of which are diagonal, and the matrix product written as a sum as follows S A = F+ S O i=l G.H.K. 1 1 1 (C.14) Summing over j of s terms (j< s) of Eq. (C. 14) an expression for A. is obtained as follows J j A.= F + > G.H.K. o 1 1 1 i=l (C.15) From Eq. (C. 15) it follows immediately that A =A. +G.H.K. J J-1 J J J (C.16) C. 3 Matrix Inversion One method of inverting the matrix A is to make use of the following matrix identity (Branin, 1959) (F+GHK)-1 F-1- F-1G(KF-1 G+H -1)-1KE1 (C.17) 169

Applying the identity (C. 17) to Eq. (C. 16) the following result is obtained. 1-l -1 -1 -1 -1 -1K -1 A. =(A. +G.H.K.) =A_-A G.(KA G+H. ) KA._ (C. 18) j A1 i j] J - j-1 j- 1 j j jj j J-1 Equation (C. 18) expresses the inverse of the matrix Aj in terms of the inverse of A. and H., G., K. Therefore, repeated application of Eq. (C. 18) j-1 J J J will yield the desired inverse of A s By means of Eq. (C. 18) the problem of inverting a matrix of order n has been changed to that of inverting a matrix of order r, the rank of H. and 1 multiplication and addition of matrices. Two extreme cases are of interest. If r= 1, the inversion indicated on the right side of Eq. (C. 1) amounts to taking the reciprocal of a number,while if none of the aij in Eq. (C. 1) are zero, Eq. (C. 12) with no partitioning becomes A -(A +GHK) 1=A -A -G(KA G+H-1 )KA (C. 19) S O 0 0 0 0 2 and thus involves finding the inverse of a matrix of order n. Thus, as is usual when the inverse of a matrix is desired, the technique to be used for the inversion process is dependent on the form of the particular matrix to be inverted. C. 4 Physical Interpretation Assuming that the system of Eq. (C. 1) is related to a physical system, the aii may be considered as the sum of the self impedance or admittance of the ith portion of the system plus the coupling impedance or admittance from other parts of the system. Thus, f.. refers to the self impedances and iii to 170

sum of the coupling impedances. The aij, of course, refer to the coupling imnpedances or admittances between the various parts of the system. This interpretation can be easily extended to the matrix Eq. (C. 14). If the systems of Eq. (C. 1) are related to a coupled mechanical-electrical system, for example, and the coupling effects between the systems are to be studied, F 0 of Eq. (C. 14) can be partitioned into two parts F =F +F o ol01 o2 Where F refers to the mechanical system and F refers to the electrical ol o2 system. Similarly, the summation on the right side of Eq. (C. 14) can be considered to be composed of three parts, that is 3 L i iKi=IHIK2 H K+ H 2K2+ G3H 3K3 il 1 where H1 contains all the mechanical coupling terms, H contains all the electrical coupling terms and H3 contains the electro-mechanical coupling terms. Thus, if G H K is not included in the summation in Eq. (C. 14) the mechanical and electrical systems are not coupled and their behavior in the absence of coupling can be determined by repeated application of Eq. (C. 16) or any standard matrix inversion technique. The effect of coupling the two systems can then be determined by adding G3H3K3 by means of Eq. (C. 16). It should be noted that H can also be partitioned and the effect of any coupling term or any group of these terms determined by repeated applications of Eq. (C. 16). 171

C. 5 Kron's Method The result of the analysis of the preceding sections is similar to that obtained by Kron (1963). However, Kron's method of dealing with large systems begins with fundamental physical laws and yields the equations of the system in the form of (C. 13) as a natural result. Because of the basic importance of Kron's technique from the standpoint of the physical laws involved, a brief discussion of this method will be given in this section. The development is couched in the terminology of electric networks, primarily because of the ease with which the final result may be obtained. This in no way implies a limitation on Kron's methods and for a more detailed discussion of the application of this method the reader is referred to references given above. The development of electric network theory is, of necessity, brief and includes only material necessary to achieve the final result. C.6 Graph Theory Several excellent treatments of linear graph theory as applied to electrical networks exist in the literature and thus only a brief summary will be given in this report. The development given here will follow roughly that given by Reed (1961) or Reed and Seshu (1961). It is assumed that the reader is familiar with such terms as oriented graph, tree, tree complement, node, etc. The number of elements in the graph will be denoted by e and the number of nodes by n. Elements in a tree of the graph will be called branches and elements in the tree complement will be called chords. Thus, for any connected graph there will be n-1 brances and e-(n-l) chords. It is well known that the connectivity relations of a linear oriented graph can be specified by various matrices formed according to different rules (Branin, 1959). For purposes of this report, three of these matrices will be of particular interest and will be defined in the following paragraphs. 172

In all cases the columns of the matrices will correspond to elements of the graph, and the columns will be arranged such that reading from left to right the first n-1 columns will correspond to branches and the last e-(n-l) columns will correspond to chords. Thus, all of the matrices will contain e columns. The incidence matrix a with elements a... The rows of this matrix 1J correspond to n-1 of the nodes of the graph, one of the nodes, called the reference node, being omitted. The elements a.. are defined as follows: 1J a..= (+1, -1, 0) if the j element is positively, negatively, not in13.th cident upon the i node. The matrix A can be partitioned according to trees and chords as As [ATA ] where AT has an inverse (Reed and Seshu, 1961). The segregate matrix S with elements s... This matrix contains n-1 ij rows corresponding to the branches of a tree. Each row of this matrix corresponds to a set of segregate elements defined as follows. If the i branch is removed the tree is divided into two subgraphs denoted by G1 and G. Let th 1 2 the i branch be positively incident on G1 (directed toward G1) contain the nodes a1, b1, cl... k and G2 contain the nodes a, b2, c2... k2. Then the th chodel o thth j chord belongs to the i segregate set if it is connected between nodes and j2. th th Thus, s..=(+1, -1) if the j element belongs to the i segregate set 1j and is positively, negatively incident on G and zero if does not belong to tne th segregate set. i segregate set. 173

Because of the manner in which S is formed, it may be partitioned S [U T Sc] (C.20) where U is a unit matrix containing n-1 rows or columns. T The mesh matrix B with elements b.. has e-(n-l) rows corresponding 1Jth to the chords of the graph. When all chords but the i are removed from the graph a single closed path is formed, the i path, and the orientation of this path is determined by the i chord. Thus, b.. is b..=(+l. ) if the j th 1i U J branch belongs to the i path and is positively (negatively) oriented with respect to this path and is zero otherwise. Because of the way it is formed B can be partitioned B= [ U]u (C.21) where U is a unit matrix containing e-(n-l) rows. c Although other connectivity relations may be defined (Branin, 1961), the above are sufficient for the purposes of this report. Indeed, the matrices A and S are closely related, and it can be shown that one can be obtained from the other by simple transformations which can be obtained from the graph. A very important relation which exists between S(A) and B is as follows: SB'- 0 (C.22) A B' 0. (C.23) The proof of this relation is well known and given in many references in the literature (Branin, 1959) (Reed and Seshu, 1961). 174

C. 7 Electrical Network Theory The mesh and node equations of electrical network theory can be developed in terms of the graph theory discussed in Section C. 6. In order to do so, however, it is necessary first to define the elements of an electrical network. In general, three elements are necessary to describe the behavior of an electrical network. These are: 1. Arbitrary voltage source denoted by the symbol e. For this element there is no relation between the terminal voltage of the element and the current through the element. The terminal voltage is an arbitrary function. A one column matrix containing B of these elements will be denoted by the symbol E; = e 2 ek 2. Arbitrary current source h. In this element the current is an arbitrary function and there is no relationship between the terminal voltage and current. A one column matrix containing k of these elements is represented by the symbol H, i.e. hk 3. Impedance or admittance elements denoted by z. or For these elements the voltage v. and current i for the i element are related by 175

or v = z. i i i i i.= y.v 1 1 1 iiy=1yi where ixi 1 A matrix representation of these elements is given by V = ZI (C.24) or I= YV (C.25) where ZY = U and where U is a unit matrix. It should be noted that the use of the terms impedance or admittance implies "steady state" or the use of a Fourier or Laplace transform. It is evident from the definition of the incidence matrix A that Kirchhoff's current law may be expressed as: AI = 0 e (C. 26) where I is a column matrix representing the currents in the elements of the e network. From the fact that A= [A A J and AT has an inverse, we can -1 multiply (C. 26) by A to obtain T -1 -1 AT AI= O = U A A I = U S = SI = 0. T e ce T c e (C. 27) 176

Thus the segregate equations may be considered as a generalized form of current law which may be stated as follows: the algebraic sum of the currents in a segregate set must be zero. If the matrix of the element voltage is represented by I, then Kirchhoff's voltage law is given by B V = 0. (C.28) e It is convenient to petition the element voltage and current matrices as follows: VIT ve V Ie = (C.29) e V 1 I c c V H H where the symbols E and H have been defined previously. When these symbols are used as subscripts, they refer to the matrix of the currents of an arbitrary voltage source and the voltage of an arbitrary current source, respectively. In addition, the subscripts T and c refer to impedance or admittance elements located in the tree or tree complement, respectively. Thus, ohms law, EqS. (C. 24) and (C. 25) may conveniently be written as V [ ] [T 0 (C. 30) V 0 Z I c. c T.. - - 177

and IT YT T 1 (C.31) I 0 Y V The segregate matrix S can be partitioned in the following manner: 0O S13 S14 (C. 32) T S23 S24 where the first n rows correspond to arbitrary voltage sources. e Also, the mesh matrix B can be partitioned B B U 0 B= 11 12 c 3) c (C. 33) B 0 -B21 B22 ~ Ull where the last nh rows correspond to arbitrary current sources. The incidence matrix A can be partitioned in a similar fashion except that unit matrices do not necessarily appear. Since the S matrix can be considered as associated with a generalized form of Kirchhoff's current low the A matrix can be obtained as a special case of a particular S matrix and thus will not be discussed further in this report. The mesh equations can be obtained from (C. 28) as follows. The partitioned form of this equation is 178

B11 B21 B12 B22 U C 0 U u ii E VT c VH = 0 (C.34) From Eq. (C. 34), we obtain B11 + [B12 Uc rv T = 0 c (C. 35) Using ohm laws in Eq. (C. 35) we obtain B116 + tB12 Uc] ZT O IT = 0 0 Z I bL c.M L. C.J (C. 36) At this point we assume I =B'I e a where the currents I are as yet unknown. a I B'll B'21 I B'11 B'21 IT B'12 22 I U 0 c c H 0 U4 (C. 37) Writing Eq. (C. 37) in detail [Ia] (C. 38) 179

From Eq. (C. 38) it is obvious that I = 1 thus (C.39) H thus IT- B11 B 22 'i (C. 40) I U 0 H c_ L CJ Using this result in Eq. (C. 36) the following equation is obtained T [12 F 01[ + Z B' H 0 [B12 C[ + 12ZTB22H ( 0 0 Z J U (C.41) C C The only unknowns in Eq. (C. 41) are the chord currents I of which C there are e-(n-l)-nh in number and thus the I may be determined if L12 C T 12 [] 0 Z U c c -1 exists. The above sets are the so called mesh equations of ordinary network theory. In order that the equations be valid, it must be shown that Kirchhoff's current low is satisfied. This can be shown easily by premultiplying (C.23) by S thus SI = SB' I e a Since SB' = 0, Kirchhoff's current law is satisfied. 180

Segregate equations. The segregate (and thus node) equations may be developed as follows. The partitioned form of Eq. (C.27) is c S13 S14 I T = 0 (C. 43) c C S24 I C H From Eq. (C.44) we obtain [c 23 LT +S24 H= 0 (C. 44) Using Eq. (C.25) in Eq. (C.44) [U 23 +S24H = 0 l YC (C.45) Assume V = S'V in detail e a U 0 V 0 U V' V (C.46) c 13 S23 V S4' S I H 14 224 181

From the above it is obvious that V = 6 V[ Thus VT 0 UT Vci L13 23. Using Eq. (C.48) in Eq. (C.45) (C.47) (C.48) V T - T c s23] LO 01 UT c 23 LI - VT+S23YcS3 + SH=0 T 23(C.13 49) (C.49) Equation (C.49) can be solved for the VT providing EUc -1 T ~0 UT Yc 223 -1 LT 23 c 2U (C.50) This result is the usual form of the segregate equations. The form of the inverse required for the mesh and segregate (node) equations, (C. 42) and (C. 50) respectively, are of the same form as Eq. (C. 13) and, thus, the same inversion technique will apply. There is one significant 182

difference between the network equations developed in this section and the development in Section C. 2. This difference lies in the coupling terms a.i and F.i(Y..). In the case of the electrical network theory it has been assumed that F..= F... Thus, if the matrix product B F TB is partitioned ij ji 12 T 12 sich that FT1 contains only Fij and the matrix product taken Fij will appear in four positions. It will appear as the coupling term between the i and jth mesh, as the coupling term between the jth and ith mesh and also in such a position that it is added to the self impedance of the ith and jth meshes. The condition the F. = F.. is not always true in network theory. Consider, for example, the mesh equations for linear analysis of a vacuum tube or transistor. (Lo et al, 1955). This situation is usually accounted for by the addition of a voltage source at an appropriate place in the equivalent circuit. In contrast to this situation, the general development presented in Section II adds each element in separately, regardless of whether it appears as a self impedance or admittance or as a coupling term. This allows somewhat more flexibility in the way the terms are added to obtain the final result. C. 8 Application to the Excitation Problem The method for the piecewise inversion of matrices can be applied to thle excitation problem discussed in Section V. The system of Eqs. (5.1) through (5.18) written in matrix form is given by Eq. (5.19) which is [L (s) f(s) = S (s) (5. 19) where 183

K y iK z x h y x z h Q ye h 1 z ( E Q. x i E Li y Q E n z (1 n F e ex 1 o e n F n e V F ey ez ez to m V. O ly 0 1m v. F. 1z l nx ONm.x 0 1 vl F. nz l F nix nfl F nfl F N m nfl and 184

V In [I'-, I 'I I I I I c 0 0 0 0 0 0 0 0 0 0 r. 0 0 0 0 0 0 0 0 0 0 0 a a 0 0 1 so Io el 0 W 0 0 1 0 -0 "I 0 6 0 0 0 P a Ets 0 a 0 a 0 a 0 a 0 0 a 0 0 1 0 9 It I 3It 0 D 0 0 1 0 to V I — b5It 0 0 0 0 co 0 0 0 0 0 1 1 I.0 B.0 le 0 Z I P C'.'ta a C) a I 0.,:tI.-. m I D.-.-P m 57 QD 0 so 0 - I I — a II — 0 0 0 0 0 0 I PD 0 Z 1(,Dclto Go 0 5I.(D 0 5 r 0 (D 0 a 0 0 a a 0 0 1 I — (D:3 to 0 M 0 0 Q) I.(D:: ta 0 (D 0 CbIM I (D:0O.0 an 5, CD (D 0 I — I (D "Ib.I — 0 1 I — (D It 0 cl I I — 0 0 "I Cs 0 0 0 0 0 0 0 0 0 0 0 0 0 B., 0 0 0 0 C) 0 a CD 0 0 0 I.0 0 0 0 0 0 a 0 0 0 0 0 0 0 c A 0 0 to 0 0 E 01% 0 0 0 0 0 a 0 a 0 1 1 E loz 0 oz 0 0 0 0 ba e 0 z 0 a 0 0 I Im 0 E Of% 0 0 0 a 0 0 0 a I I — 0 m 0 z 4I (D 0 0 z 0 0 I — a 0 z 0 I — 0 z 0 0 0 a 0 c1 't 0 0 0 0 0 0 0 0 0 a 0 0 0 0 0 0 0 a 0 0 0 0 0 a 0 C) 0 I c r. 0 0 0 to 0 0 0 0 w1 0 0 0 8 I —. v PI 0 a %-p 0 1 a c oz Io a 0 0 0 a 0 0 a 0 s 0 0 0 0 a a a 0 a 0 0 0 0 0 -A 00 ol 'I-, 00 m ts0 I b0 0 I 0 0 E 0z 0 0 0 0 0 0 0 a 0 I

lr

where a = 1 + iv + iv e ei en a.= 1 + iv. + iv. 1 ie in a = 1 + iv.+ iv n ni ne Denoting the matrix L(s)], for simplicity, bv L, the matrix L can be written as L= L + + + L L (C.51) o T E '2 v where L is a diagonal matrix composed of the elements on the main diagonal of L with the exception of the collision terms contained in a, a. and a e 1 n L is a matrix containing the thermal velocities U, U. and U, L is T e i n E the matrix containing the electric force terms, L is the matrix containing the magnetic force terms and L is a matrix containing the collision terms. These matrices are obtained in straightforward manner from L Each of the matrices LT, LE, L0 and L can now be written as a product in the form LT= GTLTdK T LE GELTd E LQ= G2L~dKQ L =G L K v v vd v where the subscript d has been used to denote a diagonal matrix. These matrix products are easily formed, although somewhat tedious to write in detail. 187

In order to illustrate the procedure, consider the matrix L. Examination of the matrix L shows that there are eight entries which contain the electron and ion gyro frequency, thus, the matrix L-t will be a diagonal matrix of order 8. Therefore, GQ will be an 18 x 8 matrix and KQ will be 8 x 18. The arrangement of the elements in LQd is usually determined by the results desired but is not critical since LQi can be partitioned and thus, L- expressed as a sum. In this case, LQd will be arranged as follows i2 cosl e -in sinO e -i2 cosO e i2 sinO Ld= e -i. cose 1 is ine 1 D2 cos i -iS sin0 i The algebraic signs of the entries in L have been preserved in L d and thus the entries in GQ and K will be +1 or zero. The matrices G and K can now be determined by inspection as follows. Consider the element in the third row and columnar of L,, i.e., -i RcosO. This term was obtained from the eleventh row and tenth column of L, hence, gll, 3=1 and k3, 10= 1. Similarly, the non-zero entries in G and K are Q 6 188

1 1 Tv 1 1 1 g1o =1 I i. 1 =1 ^0 2=1 Ko, 12=1 g11, 3=1 K3 10=1 g12. 4=1 K, 10=1 g1 3, 5=1 K5 14=1 g1 3, 6=1 K6, 15=1 g4, 7=1 K7, 13=1 g15, 13=1 K15, 13=1 The matrix L in the form of (C. 51) can now be inverted by application of Eq. (C. 18) or the terms in (C. 51) can be partitioned further. For example, again considering the L, since ~.>> i the matrix can be partitioned such that the Qe and i are separated. In this case, a solution could be obtained e neglecting Qi and then this solution modified by the addition of Q~. and, thus, the effect of the ion gyro-frequency determined. Another useful result can be obtained from (C. 51) by neglecting the collision matrix L. The matrices LT, L and L can then be partitioned in such a manner that all terms belonging to the electron plasma are in the first part of the partitioned matrix, all other terms are in the second part of the matrix. The matrix for the collisionless electron plasma has been inverted by Wu and, thus, is known. The inverse of this matrix can then be modified by use of Eq. (C. 18) to yield the general case of a three- fluid plasma including collisions. One additional remark should be made at this time about the collision matrix L. Since a, a. and a contain collision terms it is evident that the V e 1 n 189

matrix L also contains entries on the main diagonal. These main diagonal terms are, of course, coupling terms which appear in the "self impedance" and correspond to the.ii of Section C. 1. More specifically, these terms are the attenuation collision terms discussed in Section 3. 3 while the off diagonal terms are the coupling portion of the collisions. Thus, the method of piecewise inversion of matrices allows these two effects to be easily separated. The technique for the piecewise inversion of matrices is, perhaps, most useful for numerical work. Assuming that appropriate numerical values, including c and 0, have been inserted in L, this matrix is still a function of the Fourier space transform variable s. Because of this, the work in-l volved in finding L has been magnified considerably and the process of inversion can be carried out by one of two different methods. a. The inverse of L can be evaluated for a sufficient number of discrete values of s to achieve the desired result. b. The coefficients of the powers of s can be computed and thus -1 L obtained as a function of s. The method selected would depend, of course, on the form in which the result is desired and the number of mathematical operations involved. A relatively simple example of the latter procedure is discussed in Appendix B in connection with the computer program for evaluating the propagation constants. In this case, the coefficients of the polynomial were evaluated directly from the determinant defining the dispersion relation. The procedure for the piecewise inversion of matrices can be formally extended to the case of a stratified media as follows. Assume the medium to be stratified into n layers in the Z direction and the boundaries of the layers located at the points Z, < Z1, < Z2.... < Z. The solution to the set of O i n 190

Eqs. (C. 51) can be obtained for the i layer (Z. < Z< Z.) in the form of a matrix equation. This equation may then be evaluated at the edges of the layer, i. e., at Z= Zi 1 and Z.. If this procedure is carried out for each layer and the boundary conditions matched at Z= Z, Z... Z a matrix equation in the form of ( 1 ) would result. The matrix inversion technique discussed could then, in principle, be applied and the solution obtained. However, the complexity of the system of equations obtained in this mnanner is such that, in all probability, this type of analysis may be prohibitively expensive. 191

UNCLASSIFIED Security Classification DOCUMENT CONTROL DATA- R&D (Security classification of title, body of abstract and indexing annotation must be entered when the overall report is classified) 1. ORIGINATIN G ACTIVITY (Corporate author) 2a. REPORT SECURITY C LASSIFICATION University cf Michigan UNCLASSIFIED Dept of Electrical Engineering 2b. GROUP Ann Arbor, Michigan NA 3. REPORT TITLE Investigations on Excitation and Propagation in Ionized Media, 4. DESCRIPTIVE NOTES (Type of report and inclusive dates) Final Report, June 1964 through October 1965 5. AUTHOR(S) (Last name, first name, initial) Chu, Chiao-M. LaRue, John J. vanHulsteyn, David B. 6. REPORT DATE 7a. TOTAL NO. OF PAGES 7b. NO. OF REFS June 1966 208 19 8a. CONTRACT OR GRANT NO. 9a. ORIGINATOR'S REPORT NUMBER(S) AF30(602)-3381 b. PROJECT NO. 6663-1-F 5579 c. 9b. OTHER REPORT NO(S) (Any other numbers that may be assigned 557902 this report) d. RADC-TR-65-484 10. AVA IL ABILITY/LIMITATION NOTICES This document is subject to special export controls and each transmittal to foreign governments or foreign nationals may be made only with prior approval of RADC (EMLI), GAFB, N.Y. 13440. 11. SUPPLEMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY Rome Air Development Center(EMASA) _riffiss AFB NY 13440 13. ABSTRACT A study of the excitation and propagation of wave like disturbances in an ionized medium, such as the ionosphere, is made based on the linearized Euler's equation and Maxwell's equation. The local propagation constants of the basic modes of propagation are discussed. A computer program for the evaluation of these constants with given ionospheric properties is given. Methods of investigating the propagation of such waves in inhomogeneous and/or bounded media, such as ray tracing, invariant embedding, reflection and refraction, orthogonal expansion, and the use of a general matrix formulation are presented. A unified matrix-operator transform method for investigation of the excitation and propagation of disturbances in an ionized medium is proposed. 1 JAN 64 Imj4J UNCLASSI FlED DD R JAN 64 1473 UNCLAS S IFIEc D Security Classification

UN CLAS S I FIED Security Classification 14. LINK A LINK B LINK C KEY WORDS ROLE WT ROLE WT ROLE WT Electromagnetic Waves Propagation Ionosphere Ionospheric Disturbance Plasma, Homogeneous g Inhomogeneous Bounded Plasma INSTRUCTIONS 1. ORIGINATING ACTIVITY: Enter the name and address of the contractor, subcontractor, grantee, Department of Defense activity or other organization (corporate author) issuing the report. 2a. REPORT SECURITY CLASSIFICATION: Enter the overall security classification of the report. Indicate whether "Restricted Data" is included. Marking is to be in accordance with appropriate security regulations. 2b. GROUP: Automatic downgrading is specified in DoD Directive 5200.10 and Armed Forces Industrial Manual. Enter the group number. Also, when applicable, show that optional markings have been used for Group 3 and Group 4 as authorized. 3. REPORT TITLE: Enter the complete report title in all capital letters. Titles in all cases should be unclassified. If a meaningful title cannot be selected without classification, show title classification in all capitals in parenthesis immediately following the title. 4. DESCRIPTIVE NOTES: If appropriate, enter the type of report, e.g., interim, progress, summary, annual, or final. Give the inclusive dates when a specific reporting period is covered. 5. AUTHOR(S): Enter the name(s) of author(s) as shown on or in the report. Enter last name, first name, middle initial. If military, show rank and branch of service. The name of the principal,author is an absolute minimum requirement. 6. REPORT DATAE Enter the date of the report as day, month, year; or month, year. If more than one date appears on the report, use date of publication. 7a. TOTAL NUMBER OF PAGES: The total page count should follow normal pagination procedures, i.e., enter the number of pages containing information. 7b. NUMBER OF REFERENCES Enter the total number of references cited in the report. 8a. CONTRACT OR GRANT NUMBER: If appropriate, enter the applicable number of the contract or grant under which the report was written. 8b, 8c, & 8d. PROJECT NUMBER: Enter the appropriate military department identification, such as project number, subproject number, system numbers, task number, etc. 9a. ORIGINATOR'S REPORT NUMBER(S): Enter the official report number by which the document will be identified and controlled by the originating activity. This number must be unique to thi. report. 9b. OTHER REPORT NUMBER(S): If the report has been assigned any other report numbers (either by the originator or by the sponsor), also enter this number(s). 10. AVAILABILITY/LIMITATION NOTICES: Enter any limitations on further dissemination of the report, other than those imposed by security classification, using standard statements such as: (1) "Qualified requesters may obtain copies of this report from DDC." (2) "Foreign announcement and dissemination of this report by DDC is not authorized." (3) "U. S. Government agencies may obtain copies of this report directly from DDC. Other qualified DDC users shall request through (4) "U. S. military agencies may obtain copies of this report directly from DDC. Other qualified users shall request through (5) "All distribution of this report is controlled. Qualified DDC users shall request through w, b If the report has been furnished to the Office of Technical Services, Department of Commerce, for sale to the public, indicate this fact and enter the price, if known. 11. SUPPLEMENTARY NOTES: Use for additional explanatory notes. 12. SPONSORING MILITARY ACTIVITY: Enter the name of the departmental project office or laboratory sponsoring (par ing for) the research and development. Include address. 13. ABSTRACT: Enter an abstract giving a brief and factual summary of the document indicative of the report, even though it may also appear elsewhere in the body of the technical report. If additional space is required, a continuation sheet shall be attached. It is highly desirable that the abstract of classified reports be unclassified. Each paragraph of the abstract shall end with an indication of the military security classification of the information in the paragraph, represented as (TS), (S), (C), or (U). There is no limitation on the length of the abstract. However, the suggested length is from 150 to 225 words. 14. KEY WORDS: Key words are technically meaningful terms or short phrases that characterize a report and may be used as index entries for cataloging the report. Key words must be selected so that no security classification is required. Identifiers, such as equipment model designation, trade name, military project code name, geographic location, may be used as key words but will be followed by an indication of technical context. The assignment of links, rules, and weights is optional UNCLASSIFIED Security Classification