© 2010 Marco Antonio Milla STUDY OF COULOMB COLLISIONS AND MAGNETO-IONIC PROPAGATION EFFECTS ON INCOHERENT SCATTER RADAR MEASUREMENTS AT JICAMARCA BY MARCO ANTONIO MILLA DISSERTATION Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Electrical and Computer Engineering in the Graduate College of the University of Illinois at Urbana-Champaign, 2010 Urbana, Illinois Doctoral Committee: Professor Erhan Kudeki, Chair Professor Steven J. Franke Professor Jianming Jin Assistant Professor Jonathan J. Makela ABSTRACT In this dissertation, Coulomb collisions and magneto-ionic propagation effects on the incoherent scatter radar measurements have been studied and analyzed in detail. The present study aims at modeling radar observations of the equatorial ionosphere carried out at the Jicamarca Radio Observatory (Lima, Peru) using antenna beams pointed perpen- dicular to the Earth’s magnetic field B. A Monte Carlo procedure based on the simulation of charged particle trajectories in a magnetized plasma (with suppressed collective interactions) was developed to account for the effects of Coulomb collisions on the shape of the incoherent scatter spectrum. Statistics of simulated electron and ion trajectories, single-particle ACF’s, and associated Gordeyev integrals are utilized in the general framework of incoherent scatter spectrum models [e.g., Kudeki and Milla, 2006] to produce theoretical spectra for different plasma configurations. Our simulation method effectively extends the procedure of Sulzer and González [1999] into three dimensions and is valid for all magnetic aspect angles including the direction perpendicular to B. The 3D trajectories, randomized by Coulomb collisions, are described by a generalized Langevin equation with velocity-dependent friction and diffusion coefficients taken from the standard Fokker-Planck collision model of Rosenbluth et al. [1957]. A statistical analysis of the simulated trajectories shows that the ion motion is well modeled as a Brownian-motion process with Gaussian displacement distributions (and constant friction and diffusion coefficients), in which case, an analytical expression for the single-ion ACF can be obtained [e.g., Woodman, 1967]. However, the simulated electron motions do not fit a Brownian model because the electron displacement distributions in the direction parallel to B are sharper than a Gaussian. To account for these effects on our incoherent scatter spectrum model, a numerical library of electron statistics in an oxygen ii plasma (single-electron ACF’s) had to be developed. The library spans a set of densities, temperatures, and magnetic fields as needed for Jicamarca F-region applications. The antenna beams used in perpendicular-to-B radar observations at Jicamarca have angular widths of the order of a degree. Within this range of small magnetic aspect angles, different modes of magneto-ionic wave propagation are excited. These characteristic modes vary from linearly polarized in the direction perpendicular to B (Cotton-Mouton regime) to circularly polarized at aspect angles greater than 0.5◦ (Faraday rotation regime). In order to model the magneto-ionic propagation effects on incoherent scatter radar measure- ments, a computer algorithm based on the Appleton-Hartree equation for electromagnetic wave propagation in a magnetized plasma was developed. Simulation studies show that magneto-ionic propagation effectively modifies the shapes of the radar beams and does have an impact on the incoherent scatter radar measurements because the polarization of the incident and backscattered fields vary as they propagate through the ionosphere. A soft-target radar equation, which incorporates our collisional incoherent scatter spec- trum and magneto-ionic propagation models, is formulated to model the radar measure- ments collected at Jicamarca. Voltages detected by the radar antenna are represented as the beam-weighted sum of ionospheric backscattered signals corresponding to the range of magnetic aspect angle directions illuminated by the antenna beam. This integration is carried out numerically using a finite-element-like integration method that takes advantage of the slow variation of physical parameters in the direction transverse to the geomagnetic field. The resultant radar model is utilized in the inversion of ionospheric parameters in a three-beam radar experiment conducted at Jicamarca. The experiment interleaves radar observations with perpendicular-to-B and off-perpendicular antenna beams. The data model matches very closely the different features of the measured data; for instance, it predicts the enhancement of the measured power in the direction perpendicular-to-B at ionospheric altitudes where the electron temperature is greater than the ion temperature. F-region electron density and temperature ratio (Te/Ti) estimates were obtained using a least-squares inversion algorithm. The inversion results show a good agreement with ionosonde data, validating our model for incoherent scatter radar measurements. iii To Silvanita iv ACKNOWLEDGMENTS I would like to express my deepest gratitude to my adviser, Professor Erhan Kudeki, for his guidance during the development of this dissertation. I greatly appreciate the time he spent teaching me the “mysteries” of incoherent scatter. I very much enjoyed our conversations and discussions on various topics of science; they were always insightful and challenged my curiosity. I would also like to thank Professors Steven Franke, Jianming Jin, and Jonathan Makela for agreeing to be members of my doctoral committee and for their helpful suggestions and comments. I am very thankful to the staff at the Jicamarca Radio Observatory, who always put in great effort and show dedication and commitment to their work. Their professional assistance was crucial during the experiments we conducted at the Observatory. I am very much indebted to Dr. Jorge Chau for his support during my visits to Jicamarca. I would also like to thank Dr. Ronald Woodman for our conversations on topics related to this study. I would like to thank my family for all their support throughout my years at the University of Illinois. I have dedicated this dissertation to my wife, Silvanita; without her love and patience, I would not have finished this work. I acknowledge the use of the Turing Cluster maintained and operated by the Computa- tional Science and Engineering Program at the University of Illinois. The Turing Cluster is a 1536-processor Apple G5 X-serve cluster devoted to high-performance computing in engineering and science. This work was funded by the NSF grant ATM-0215246 to the University of Illinois. v TABLE OF CONTENTS Page LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Perpendicular-to-B Incoherent Scatter Radar Measurements at Jicamarca . 1 1.2 Studying and Modeling Coulomb Collisions and Magneto-Ionic Propagation Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 CHAPTER 2 INCOHERENT SCATTER SPECTRUM, LANGEVIN EQUATION, AND FOKKER-PLANCK COLLISION MODEL . . . . . 10 2.1 Soft-Target Radar Equation and Ambiguity Function . . . . . . . . . . . . . 10 2.2 General Framework of IS Spectral Models . . . . . . . . . . . . . . . . . . . 12 2.3 Single-Particle ACF and Langevin Equation . . . . . . . . . . . . . . . . . . 14 2.4 The Link between Fokker-Planck and Langevin Equations . . . . . . . . . . 17 2.5 Spitzer Friction and Diffusion Coefficients . . . . . . . . . . . . . . . . . . . 18 2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 CHAPTER 3 MONTE CARLO COMPUTATION OF PARTICLE TRAJECTORIES AND SINGLE-PARTICLE ACF’S . . . . . . . . . . . 21 3.1 Computer Simulations of Plasma Particle Trajectories . . . . . . . . . . . . 21 3.2 Estimation of the Single-Particle ACF’s and Gordeyev Integrals . . . . . . . 26 3.3 Statistics of the Ion and Electron Trajectories . . . . . . . . . . . . . . . . . 31 3.3.1 Statistics of the Ion Displacements . . . . . . . . . . . . . . . . . . . 33 3.3.2 Statistics of the Electron Displacements . . . . . . . . . . . . . . . . 37 CHAPTER 4 COLLISIONAL INCOHERENT SCATTER SPECTRUM: RESULTS AND COMPARISONS . . . . . . . . . . . . . . . . . . . . . . . 43 4.1 Simulated Electron Gordeyev Integrals . . . . . . . . . . . . . . . . . . . . . 43 4.2 Simulated Collisional IS Spectra . . . . . . . . . . . . . . . . . . . . . . . . . 45 CHAPTER 5 MAGNETO-IONIC PROPAGATION EFFECTS ON THE INCOHERENT SCATTER RADAR MEASUREMENTS . . . . . . . . . 53 5.1 Propagation of Electromagnetic Waves in a Homogeneous Magnetoplasma . 53 5.2 Model for Radiowave Propagation in an Inhomogeneous Ionosphere . . . . . 62 5.3 Soft-Target Radar Equation and Magneto-Ionic Propagation . . . . . . . . . 71 vi CHAPTER 6 MAGNETO-IONIC EFFECTS ON SIMULATED AND MEASURED ISR DATA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.1 The Differential-Phase Radar Configuration . . . . . . . . . . . . . . . . . . 74 6.2 Magneto-Ionic Effects on the Radiation Patterns . . . . . . . . . . . . . . . 79 6.3 Simulated and Measured Power and Cross-Correlation Profiles . . . . . . . . 85 6.4 Magneto-Ionic Effects on Spectrum and Cross-Spectrum Measurements . . . 89 CHAPTER 7 THREE-BEAM EXPERIMENTS AT JICAMARCA . . . 93 7.1 Experiment Description and Antenna Configurations . . . . . . . . . . . . . 94 7.2 Radar Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 7.3 Data Model and Inversion Technique . . . . . . . . . . . . . . . . . . . . . . 107 7.4 Results and Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 CHAPTER 8 CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . 123 APPENDIX A PARTICLE DYNAMICS DESCRIPTION OF “BGK COLLISIONS” AS A POISSON PROCESS . . . . . . . . . . . . . . . . . . 127 A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 A.2 Derivation of the BGK Gordeyev Integral following a Particle Dynamics Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 A.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 APPENDIX B FRICTION AND DIFFUSION COEFFICIENTS FOR “TEST PARTICLES” IN NON-ISOTHERMAL PLASMAS . . . . . . . 137 B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 B.2 Steady-State Solution of the Fokker-Planck Equation with Spitzer Collision Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 B.3 Modified Spitzer Friction Coefficient . . . . . . . . . . . . . . . . . . . . . . 145 B.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 AUTHOR’S BIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 vii LIST OF TABLES Table Page 3.1 Typical plasma parameters for an equatorial F-region ionosphere. . . . . . . 22 3.2 Representative physical parameters of a plasma. Definitions and sample values are provided. The values were computed for an O+ plasma with parameters given in Table 3.1. Note that the plasma and gyro frequency formulas are given for angular frequencies, but the sample values are given in Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3 Plasma parameters spanned by the library of single-electron ACF’s. . . . . 30 7.1 Radar operating parameters considered in the 3Ba and 3Bb experiments. . . 99 viii LIST OF FIGURES Figure Page 3.1 Sample trajectory of an electron moving in an O+ plasma with density Ne = 1012 m−3 and temperatures Te = Ti = 1000 K. The presence of a uniform magnetic field B parallel to the z-axis was considered in the simulation. The intensity of the field was assumed to be 25µT. On the top, the left panel shows the particle trajectory in 3D space, while the panel on the right is the projection of the trajectory on the xy-plane (i.e., the plane perpendicular to B). On the bottom, the displacement of the electron in the direction parallel to B is displayed as a function of time. The trajectory was sampled every ∆t = 0.1µs. A total of 4×103 samples are displayed corresponding to a time interval of 400µs. In all three panels, the red dots show the starting location of the particle (the origin). In addition, the red curve depicts the first 1.4µs of the simulated trajectory (about one gyroperiod). . . . . . . . . 25 3.2 Probability distributions of (a) ion and (b) electron velocities resulting from the simulations. In each plot, the pdf’s of the velocity components (dis- played in colors) are compared to a Gaussian distribution (black line). No- tice that the velocity axes are normalized by the thermal speeds of the corresponding particle species. The plasma parameters considered in the simulations are given in Table 3.1. . . . . . . . . . . . . . . . . . . . . . . . 32 3.3 Probability distributions of the displacements of a simulated ion in the di- rections perpendicular (top panels) and parallel (bottom panels) to the mag- netic field. On the left, the displacement pdf’s are displayed as functions of time delay τ . On the right, sample cuts of the pdf’s are compared to a Gaussian distribution. Note that all distributions at all time delays are normalized to unit variance. The displacement axis of each distribution at every delay τ is scaled with the corresponding standard deviation of the simulated displacements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4 Variances of the displacements of a simulated ion in the directions perpen- dicular (top panels) and parallel (bottom panels) to the magnetic field. The simulation results (color lines) are displayed for two time intervals: 5 ms (left panels) and 500 ms (right panels). The dashed lines correspond to the fits of the results to the theoretical expressions for 〈∆r2⊥〉 and 〈∆r2‖〉 of the Brownian-motion model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.5 Simulated single-ion ACF’s at different magnetic aspect angles α for two radar Bragg wavelengths: (a) λB = 3 m and (b) λB = 0.3 m. The simulation results (color lines) are compared to theoretical ACF’s computed using ex- pression (3.22) of the Brownian-motion approximation (dashed lines). Note that there is effectively no dependence on aspect angle α. . . . . . . . . . . 36 ix 3.6 Same as Figure 3.3 but for the case of a simulated electron. All distributions at all time delays are normalized to unit variance. Note that the distribu- tions of the displacements parallel to B become narrower than a Gaussian distribution in less than a millisecond. . . . . . . . . . . . . . . . . . . . . . 38 3.7 Same as Figure 3.4 but for the case of a simulated electron. In this case, the results are displayed for time intervals of 0.02 ms (left panels) and 10 ms (right panels). Note the different scales for the displacement variances. . . . 39 3.8 Simulated electron ACF’s for λB = 3 m at different magnetic aspect angles: (a) α = 0◦, (b) α = 0.01◦, (c) α = 0.05◦, (d) α = 0.1◦, (e) α = 0.5◦, and (f) α = 1◦. Note the different time scales used in each plot. . . . . . . . . . . . 40 4.1 Simulated electron Gordeyev integrals as functions of Doppler frequency and magnetic aspect angle for two radar Bragg wavelengths: (a) λB = 3 m and (b) λB = 0.3 m. An O+ plasma with electron density Ne = 1012 m−3, temperatures Te = Ti = 1000 K, and magnetic field Bo = 25µT is considered. 44 4.2 Simulated incoherent scatter spectra as a function of magnetic aspect angle and Doppler frequency for λB = 3 m (e.g., Jicamarca radar Bragg wave- length). An O+ plasma with physical parameters given in Table 3.1 is considered. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.3 Sample cuts of the simulated spectra of Figure 4.2 taken at different mag- netic aspect angles: (a) α = 0◦, (b) α = 0.01◦, (c) α = 0.05◦, (d) α = 0.1◦, (e) α = 0.5◦, and (f) α = 1◦. Our simulation results (blue lines) are com- pared to IS spectra computed using the Brownian-motion model (green lines), the collisionless electron model (red lines), the simulation results of Sulzer and González [1999] (light-blue lines), and the aspect angle depen- dent collision frequency model of Woodman [2004] (magenta lines). Note the different frequency scales used in each plot. . . . . . . . . . . . . . . . . 47 4.4 Electron density dependence of the simulated IS spectra for λB = 3 m at aspect angles α = 0◦ (left panel), α = 0.1◦ (central panel), and α = 1◦ (right panel). An O+ plasma in thermal equilibrium with Te = Ti = 1000 K and Bo = 25µT is considered. Note the different frequency scales used in each plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.5 Te/Ti dependence of the simulated IS spectra for λ = 3 m at aspect angles α = 0◦ (left panels), α = 0.02◦ (central panels), and α = 0.1◦ (right pan- els). The same spectra are plotted on the top and bottom panels; however, on the bottom, every spectrum has been normalized to its peak value. An O+ plasma with Ne = 1012 m−3, Ti = 1000 K, and Bo = 25µT is consid- ered. Each curve corresponds to different values of electron temperatures Te (1000 K, 1500 K, 2000 K, 2500 K, and 3000 K). Note the different frequency scales used in each column. . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.6 Electron scattering efficiency factor η(k) resulting from the frequency inte- gration of the simulated incoherent scatter spectra as a function of electron- to-ion temperature ratio Te/Ti and magnetic aspect angle α. An O+ plasma with Ne = 1012 m−3, Ti = 1000 K, and Bo = 25µT is considered. . . . . . . 52 x 5.1 Coordinate system used for analyzing wave propagation in a magnetized plasma. The magnetic field Bo is on the yz-plane at an angle θ from the propagation vector k which is parallel to the zˆ-direction. The wave field E has three mutually orthogonal components Ek, Eθ, and Eφ in directions kˆ = zˆ, θˆ = −yˆ, and φˆ = xˆ, respectively. . . . . . . . . . . . . . . . . . . . 56 5.2 Geometry of wave propagation in an inhomogeneous magnetized ionosphere. 63 5.3 Electron density and Te/Ti profiles as functions of height. . . . . . . . . . . 66 5.4 Polarization coefficients for the mean square voltages detected by a pair of orthogonal linearly polarized antennas placed at Jicamarca. The antennas have very narrow beams and scan the ionosphere from north to south prob- ing different magnetic aspect angle directions. Note that, for most pointing directions, the polarization of the detected fields rotates (Faraday rotation effect), except in the direction where the beam is pointed perpendicular to B, in which case, the type of polarization changes from linear to circular (Cotton-Mouton effect). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.5 Co-polarized (left panel), cross-polarized (middle panel), and total (right panel) backscattered power detected by a pair of orthogonal linearly po- larized antennas (see caption of Figure 5.4). Power levels are displayed in units of electron density. In each plot, the dashed white lines indicate the directions half a degree away from perpendicular to B, while the continuous lines correspond to the directions one degree off. . . . . . . . . . . . . . . . 68 5.6 An example of the ALTAIR scan measurements collected during the 2004 EQUIS 2 campaign. In colors, levels of backscattered power calibrated to match electron densities in the probed region are displayed. . . . . . . . . . 70 6.1 Antenna diagram of the differential-phase technique developed for the Ji- camarca incoherent scatter radar system. In this mode of operation, the Jicamarca antenna is phased to point perpendicular to B. In transmission, the southeast polarized down antenna quarters are excited. In reception, backscattered fields detected by the up and down polarizations of the north and south quarters are combined to synthesize meridional and zonal polar- ized field components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.2 Two-way antenna radiation pattern of the differential-phase radar configu- ration (in free-space). In the left panel, the radiation pattern (in dB units relative to its peak value) is displayed as a function of direction cosines with respect to the x- and y-axes of the Jicamarca frame of reference. In the right panel, the pattern is displayed as a contour plot in a rotated coordi- nate system in which the x- and y-axes are now aligned with the east and north directions (the Jicamarca frame of reference is rotated 51◦ positive azimuth). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.3 Meridional and zonal effective radiation patterns for r = 0, 100, 150, 200 km. 82 6.4 Same as Figure 6.3 but for r = 250, 300, 350, 400 km. . . . . . . . . . . . . 83 6.5 Same as Figure 6.3 but for r = 450, 500, 550, 600 km. . . . . . . . . . . . . 84 6.6 Simulated incoherent scatter signal power and coherence (magnitude and phase) profiles for the differential-phase radar configuration of June 8, 2004. 87 xi 6.7 Measured incoherent scatter signal power and coherence (magnitude and phase) profiles that correspond to the differential-phase experiment con- ducted on June 8, 2004. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.8 Incoherent scatter Doppler spectra (left) and north-south baseline coher- ence spectra (right) measured with the differential-phase radar technique at Jicamarca on June 8, 2004. . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.1 Antenna configurations of the 3Ba (left panel) and 3Bb (right panel) exper- iments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 7.2 Two-way radiation patterns of the west (top row) and east (bottom row) beams used in the 3Ba (left) and 3Bb (right) experiments. The patterns are displayed as functions of direction cosines in a coordinate system in which x- and y-axes are aligned with the east and north directions at Jicamarca. Note that in the 3Ba mode, the patterns of the west beams are identical, but in the 3Bb mode, the west beams have different shapes. The same can be observed in the case of the east beams. . . . . . . . . . . . . . . . . . . . 97 7.3 Two-way radiation patterns of the south beams used in the 3Ba (left) and 3Bb (right) experiments. Note that, in the north-south direction, the 3Ba south beam is narrower than the 3Bb beam. . . . . . . . . . . . . . . . . . . 98 7.4 Sample self- and cross-spectra measured in the 3Ba experiment of June 19, 2008 (top rows), and in the 3Bb experiment of January 9, 2009 (bottom rows). The sample data were obtained after five minutes of integration. The self- and cross-spectra are plotted as functions of Doppler velocity and radar range. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.5 Range vs. time plots of the signal power and cross-correlation data measured in the 3Ba experiment of June 19, 2008. On the left, the power data collected by each of the radar beams are displayed in linear scale. On the right, the magnitudes of the cross-correlation data are also plotted in linear scale, while the phase data are plotted in degrees. . . . . . . . . . . . . . . . . . . 103 7.6 Same plots as described in Figure 7.5, but for the signal power and cross- correlation data measured in the 3Bb experiment of January 9, 2009. . . . . 106 7.7 Sample plots of the results obtained from the inversion of the 3Ba and 3Bb data (top and bottom panels, respectively). The estimated Ne and Te/Ti profiles are plotted in the left column. In the other columns, the power and cross-correlation measurements (dotted points) are compared with the modeled profiles (continuous lines). . . . . . . . . . . . . . . . . . . . . . . . 117 7.8 Range vs. time plots of the electron density and Te/Ti estimates obtained from the inversion of incoherent scatter signal power and cross-correlation data collected in the 3Ba experiment of June 19, 2008 (left column), and also in the 3Bb experiment of January 9, 2009 (right column). . . . . . . . . 118 7.9 Magnitudes and phases of the radar calibration constants obtained as part of the inversion of the 3Ba and 3Bb radar data. . . . . . . . . . . . . . . . 120 7.10 Goodness-of-fit values resulting from the inversion of the 3Ba and 3Bb radar data. The blue curves correspond to the values obtained by fitting only for the electron density Ne and assuming that Te/Ti = 1 at all altitudes. The green curves correspond to the values obtained by fitting for both parame- ters, Ne and Te/Ti. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 xii 7.11 Comparison of the peak values of plasma frequencies (foF2) measured with the 3Ba and 3Bb radar techniques (blue lines) and with an ionosonde system also located at Jicamarca (red dots). . . . . . . . . . . . . . . . . . . . . . . 122 7.12 Comparison of the altitudes corresponding to the peak values of plasma frequency presented in Figure 7.11. . . . . . . . . . . . . . . . . . . . . . . . 122 B.1 Relative difference between the effective electron temperature T effe computed using (B.31) and the value of Te considered in the calculation of the Spitzer coefficients for an oxygen plasma. The difference is displayed as a function of Te/Ti. Each curve corresponds to different constant values of Te. . . . . . 144 B.2 Same as Figure B.1 but for the case of the effective ion temperature T effi . . . 145 xiii CHAPTER 1 INTRODUCTION 1.1 Perpendicular-to-B Incoherent Scatter Radar Measurements at Jicamarca The incoherent scatter radars are the most powerful ground-based facilities built for the study of the Earth’s ionosphere. Using these instruments, ionospheric plasma pa- rameters such as densities and temperatures can be measured remotely as functions of altitude. Only a small group of these radars is distributed around the World, and among them, the Jicamarca radar is special by its location. Situated near Lima, Peru, the 50 MHz Jicamarca radar is the only incoherent scatter (IS) facility that is routinely used for probing the equatorial ionosphere. Since it began operating in 1961, the Jicamarca measurements have become the primary resource for scientists to conduct studies of the equatorial electrodynamics. One of the most important applications of the Jicamarca radar is the measurement of ionospheric plasma drifts. For this purpose, the Jicamarca antenna array is phased to point perpendicular to the Earth’s magnetic field B. These measurements are possible on a pulse-to-pulse basis (as opposed to ACF measurements conducted with multi-pulse techniques) because the bandwidth of the incoherent scatter radar (ISR) spectrum narrows considerably when the radar beam is pointed perpendicular to B. The original technique of Woodman and Hagfors [1969] has been recently improved by Kudeki et al. [1999]. Very ac- curate plasma drifts determined from the Doppler shifts of the sharp peaks of the measured spectra are obtained using the new technique. The procedure developed by Kudeki et al. [1999] consists in fitting the measurements to a simplified ISR spectrum formula. This 1 equation results from the beam-weighted integration of theoretical spectra corresponding to the range of magnetic aspect angles illuminated by the radar beam. The magnetic aspect angle is the complement of the angle between the scattered field wavevector k and the ambient magnetic field B. Kudeki et al. [1999] developed their model based on the well-established collisionless IS theory [Fejer , 1961; Farley et al., 1961, and others] and the collisional IS theory of Woodman [1967]. Analysis and comparisons of these theories led Kudeki and his team to assume that electron Coulomb collisions could be neglected in the formulation of the simplified model (though ion Coulomb collisions should still be con- sidered). The simplified model predicted that the frequency bandwidth of the measured spectrum would be proportional to the plasma temperature. Because of this, the initial expectation of Kudeki and his team was to measure ionospheric temperatures (as well as drifts) by applying their spectral estimation technique. To their surprise, however, the temperature estimates they obtained were about half of what is expected for an F-region ionosphere. Bhattacharyya [1998] carried out additional theoretical and experimental work in an attempt to attribute the low temperature estimates to a modification of the spectrum caused by magneto-ionic propagation effects. Despite these efforts, the temperatures ob- tained were still underestimated. All this work gave Kudeki and his team enough evidence to conclude that the actual measured ISR spectrum was in fact narrower than what the theories available at that time were able to predict, and in consequence, a revision of the IS theory for modes propagating perpendicular to B was needed. The very same year, Sulzer and González [1999] revised the IS theory and showed that electron Coulomb collisions do have a significant impact on determining the shape of the IS spectrum at small magnetic aspect angles. In the literature, the effects of Coulomb collisions on incoherent scattering have been studied before, but using simplified models. Based on a qualitative analysis, Farley [1964] recognized that ion Coulomb collisions would be responsible for the lack of O+ gyroresonance signatures on IS observations at F-region heights. This analysis was later verified by the collisional IS theory of Woodman [1967], theory that was based on the simplified Fokker-Planck collision model of Dougherty [1964]. These two studies concluded that electron Coulomb collisions would have very small effects 2 on IS observations at small magnetic aspect angles. However, Sulzer and González [1999] showed that this is not the case. Based on the more complex Fokker-Planck Coulomb collision model of Rosenbluth et al. [1957], Sulzer and González carried out intensive com- putations to accurately model the effect of electron collision on the IS spectra for a set of plausible ionospheric plasma parameters. Their results verified that the collisional spec- trum is indeed narrower than what the collisionless theory predicts, a result that is in better agreement with Jicamarca observations. Furthermore, they found that the effect of electron collisions extends up to relatively “large” magnetic aspect angles (as large as 3 or 4 degrees). This important result was the key to solve another puzzle of Jicamarca IS observations (using antenna beams pointed off-perpendicular to B), which was the sys- tematic measurement of electron temperatures lower than ion temperatures at nighttime when thermal equilibrium is expected [Aponte et al., 2001]. The collisional incoherent scatter spectrum calculations of Sulzer and González are based on the statistics of simulated electron trajectories randomized by Coulomb colli- sions. In order to simplify their computations, Sulzer and González assumed that the geomagnetic field is sufficiently strong to keep the electrons moving along the same ge- omagnetic field lines, such that diffusion across the field could be neglected. Under this assumption, the problem can be reduced to the simulation of the electron motion in a single dimension (i.e., in the direction parallel to B). This “strong” magnetic field approx- imation is valid for Jicamarca IS observations at aspect angles larger than 0.1◦ (limit of Sulzer and González simulations). However, for radar observations perpendicular to B, an IS spectrum model valid for zero magnetic aspect angle is required. Given the small range of magnetic aspect angles that remains to be modeled, one may think that extrapolation of Sulzer and González spectrum would be just good enough; however, it is in this range of aspect angles (between 0◦ and 0.1◦) that the spectrum becomes “super” sharp with its bandwidth being reduced by a factor larger than ten. This is the cause of the sharp spec- trum measured at Jicamarca with the antenna pointed perpendicular to B. Therefore, in order to quantify all the features of the measured IS spectrum, an accurate model for all magnetic aspect angles is needed. 3 More recently, Kudeki et al. [2003] developed a new incoherent scatter radar technique for the measurement of ionospheric plasma densities using antenna beams pointed perpen- dicular to B. An electromagnetic wave propagating through the ionosphere experiences changes in its polarization caused by the presence of the Earth’s magnetic field. These magneto-ionic propagation effects are significant for detection at 50 MHz [Farley , 1969]. The technique consists of transmitting radiowave pulses using a single linearly polarized antenna. The waves excite both the ordinary and extraordinary modes of propagation in a magneto-ionic medium. In reception, two orthogonal linearly polarized antennas are used to collect the backscattered signals. The differential phase of the signals detected by both orthogonal antennas depends on the difference between the refractive indices of the modes of propagation. The model of these measurements was developed based on the Appleton- Hartree equation for electromagnetic wave propagation in a magnetized plasma. Inversion of the radar measurements yields electron density estimates as a function of height. Addi- tional details of the experiment, the data model, and the estimation technique are given in Feng et al. [2003] and Feng et al. [2004]. Depending on the direction of propagation, radiowaves traveling through the iono- sphere experience different types of magneto-ionic effects [Yeh et al., 1999a, b]. If the magnetic field is perpendicular to the propagation direction, changes in the polarization are described by the Cotton-Mouton effect, and the characteristic modes of propagation are linearly polarized. At relatively “large” magnetic aspect angles, the Faraday effect takes place and the characteristic modes become circularly polarized. In the case of radiowave propagation at 50 MHz, the transition between these two regimes happens at a critical angle of approximately 0.5◦ measured with respect to the perpendicular-to-B direction. As the Jicamarca antenna beam has a finite angular width, different modes of propagation for small magnetic aspect angles are excited, modes that in general are elliptically polarized. The resultant magneto-ionic effects on the differential-phase measurements were modeled by Feng et al. [2003]. The model gets further complicated as the incoherent scatter radar cross section (RCS) is also dependent on magnetic aspect angle. When the plasma is not in thermal equilibrium, i.e., when the electron temperature is greater than the ion 4 temperature, the incoherent scatter RCS becomes larger at small magnetic aspect angles [Farley , 1966]. This feature of the IS signal was not fully modeled by Feng et al. [2003]. In their approach, they considered a heuristic expression for the RCS since the IS theories available at that time were not appropriate for this type of radar measurements. The differential phase technique can be improved if an accurate incoherent scatter model for modes propagating perpendicular to B is applied. 1.2 Studying and Modeling Coulomb Collisions and Magneto-Ionic Propagation Effects In this dissertation, Coulomb collisions and magneto-ionic propagation effects on the incoherent scatter radar measurements have been studied and analyzed in detail. This study aims at modeling radar observations of the equatorial ionosphere carried out at the Jicamarca Radio Observatory (Lima, Peru) using antenna beams pointed perpendicular to the Earth’s magnetic field B. The result of our studies is the development of a new incoherent scatter spectrum model valid for all magnetic aspect angles (including the direction perpendicular to B). This work is a continuation of the efforts of Bhattacharyya [1998] and Feng [2002] to improve the capabilities of the incoherent scatter Jicamarca radar. Our investigations were divided in two stages: 1. The study of the effects of Coulomb collisions in determining the shape of the inco- herent scatter spectrum. 2. The modeling of the incoherent scatter radar spectrum considering magneto-ionic propagation effects. The first stage of this research involved the development of a procedure to model the incoherent scatter spectrum including the effects of Coulomb collisions. The proposed method is based on the stochastic simulation of the trajectories of charged particles em- bedded in a plasma. This procedure is effectively an extension of the method of Sulzer and González [1999] into three dimensions. In the simulations, particle trajectories are 5 randomized by Coulomb interactions as described by the Fokker-Planck collision model of Rosenbluth et al. [1957]. A statistical analysis of the simulated trajectories shows that the electron statistics are poorly approximated by simplified collision models. The application of this procedure in routine computations of collisional IS spectra is not practical because of the large amount of computations that it requires. In order to circumvent this difficulty, a numerical library of the statistics of the electron trajectories has been developed. The library spans a set of densities, temperatures, and magnetic fields as needed for Jicamarca F-region applications. The second stage of this project corresponds to the modeling of the incoherent scatter beam-weighted radar spectrum measured with antennas pointed perpendicular to B. At 50 MHz, the polarizations of radiowaves propagating through the ionosphere are modified due to magneto-ionic effects. Therefore, the consideration of these effects in IS spectrum models is necessary. Further complications arise because, within the range of small magnetic aspect angles illuminated by the antenna beam, a variety of characteristic modes of propagation are excited. These effects result in a distortion of the shape of the propagating wavefronts and, thus, in a modification of the measured IS signals. On the other hand, the rapid variation of the IS spectrum at small magnetic aspect angles imposes an additional level of complexity to the model. Because of these features of the IS signal, the shape of the radar beam has to be considered in the modeling, as it effectively determines the shape of the measured spectrum. A model of the beam-weighted incoherent scatter radar spectrum that takes into account magneto-ionic propagation and collisional effects was developed. A comparison of measured and simulated incoherent scatter data validates to some extent the accuracy of our model. The ultimate application of this model will be the estimation of ionospheric physical parameters with Jicamarca antenna beams pointed perpendicular to the geomagnetic field. 1.3 Dissertation Outline The present dissertation is organized as follows. The general framework for incoherent scatter theories is introduced in Chapter 2. The framework formulates the spectrum of 6 incoherently scattered radiowaves in terms of the so-called Gordeyev integrals [Kudeki and Milla, 2010]. As discussed in Section 2.2, the Gordeyev integral is defined as the one-sided Fourier transform of the characteristic function of the displacements of an electron or an ion in the absence of collective interactions. As Coulomb collisions modify the particle trajectories, collisions do have an impact on the shapes of these characteristic functions, also termed single-particle ACF’s. It is through the modeling of the particle trajectories that the effects of Coulomb collisions are introduced in our incoherent scatter spectrum model. The motion of a charged particle embedded in a plasma is modeled by a generalized version of the Langevin equation. In this equation, the effect of Coulomb collisions is described by the action of a friction force and random diffusion forces. The expected values of these forces are given by the Fokker-Planck friction and diffusion coefficients of Rosenbluth et al. [1957]. The general relationship between the Fokker-Planck and the generalized Langevin equations is also discussed in Chapter 2. A Monte Carlo procedure was developed to simulate the particle trajectories in a magnetized plasma. The approach effectively extends the procedure of Sulzer and González [1999] into three dimensions by letting the charged particles diffuse across the magnetic field lines. Since the Langevin equation models the particle trajectories as Markov processes, discrete time-update equations for the velocity and the position of a test particle can be readily obtained. The computer program for the simulation of the particle trajectories was developed based on these equations. The program is outlined in Chapter 3. Both electron and ion trajectories were simulated. The use of these trajectories for the numerical estimation of the single-particle ACF’s and the calculation of the corresponding Gordeyev integrals are described in Section 3.2. The statistics of our computed trajectories deserve a careful analysis. In Section 3.3, velocity and displacement distributions resultant from the simulations are presented. We have verified that the probabilistic distributions of the ion displacements in the directions parallel and perpendicular to B can be very well approximated by independent Gaussian distributions. In consequence, the ion trajectory can be described by a Brownian motion 7 process with constant friction and diffusion coefficients. This is convenient because, in this particular case, an analytical expression for the Gordeyev integral can be derived [Wood- man, 1967; Holod et al., 2005, and others]. On the other hand, however, electron statistics show a different behavior. While in the direction perpendicular to B the displacement distribution can still be approximated by a Gaussian function, in the parallel direction the distribution is different. For this reason, the electron motion cannot be modeled as a simple Brownian motion process. Moreover, no analytical expression for the corresponding Gordeyev integral was found; therefore, a numerical library of electron Gordeyev integrals for an oxygen plasma had to be developed. The library spans a set of densities, tempera- tures, and magnetic fields as needed for Jicamarca F-region applications. In Chapter 4, the computed electron and ion Gordeyev integrals are utilized in the general framework of IS models to produce theoretical IS spectra for different magnetic aspect angles, including the direction perpendicular to B. The analysis of our results shows that the modeled spectrum saturates at very small magnetic aspect angles (< 0.01◦). Additionally, we have also verified that the width of the perpendicular-to-B spectrum is proportional to the electron Coulomb collision frequency. As the value of this collision frequency is inversely proportional to the electron temperature, we found that the modeled spectrum narrows at high temperatures and widens at low temperatures; this result is discussed further in Section 4.2. Moreover, the dependence on temperature of our simulated spectra for different magnetic aspect angles is analyzed at the end of this section. The second stage of our investigations corresponds to the study of the magneto-ionic propagation effects on the IS radar signals. As mentioned above, a radiowave propagat- ing through the ionosphere experiences changes in its polarization due to the presence of the Earth’s magnetic field. To account for these magneto-ionic propagation effects on perpendicular-to-B radar observations at Jicamarca, a model for incoherent scatter spectrum and cross-spectrum measurements is formulated in Chapter 5. In this model, the voltages detected by the radar antennas are represented as beam-weighted sums of the ionospheric backscattered signals arriving from the range of magnetic aspect angle directions illuminated by the antenna beams. The model is effectively an extension of 8 the Appleton-Hartree solution for a homogeneous magneto-plasma to the case of an inho- mogeneous ionosphere. Simulation studies based on our model show that magneto-ionic propagation effectively modifies the shapes of the radar beams and does have an impact on the IS radar measurements as described in further detail in Chapter 6. In Chapter 7, the incoherent scatter data model of Chapter 5 is utilized in the inversion of ionospheric physical parameters from radar data collected in three-beam radar exper- iments conducted at Jicamarca. These experiments interleave radar observations with perpendicular-to-B and off-perpendicular antenna beams. The data model matches very closely the different features of the measured data, e.g., it predicts the enhancement of the measured power in the direction perpendicular-to-B at ionospheric altitudes where the electron temperature is greater than the ion temperature. F-region electron density and temperature ratio (Te/Ti) estimates were obtained by applying a least-squares inversion algorithm. The estimated results show a good agreement with ionosonde data, validating our model for incoherent scatter radar measurements. This dissertation is concluded in Chapter 8 with a discussion of the main results and future extensions of our studies. The ultimate application of the incoherent scatter model developed in this dissertation will be the estimation of ionospheric temperatures, simulta- neously with electron densities and plasma drifts, using antenna beams pointed perpen- dicular to the geomagnetic field at Jicamarca. This will be the fulfillment of objectives set forth more than a decade ago, when spectral temperature estimations were first attempted [Bhattacharyya, 1998]. 9 CHAPTER 2 INCOHERENT SCATTER SPECTRUM, LANGEVIN EQUATION, AND FOKKER-PLANCK COLLISION MODEL The problem of modeling and calculating the IS spectrum was first addressed by Dougherty and Farley [1960], Fejer [1960], Salpeter [1960], Hagfors [1961], and others. Different methods of derivation were followed but identical results were obtained in equiv- alent physical settings. The fundamental principles that link these different approaches constitute a framework for incoherent scatter spectral theories. A detailed description of this general framework is presented by Kudeki and Milla [2010]. Aspects of the theory relevant to our work are discussed in the first two sections of the present chapter. The framework offers the alternative of formulating the IS spectrum problem in terms of the statistics of single-particle trajectories. Then, to include the effects of Coulomb collisions into the problem, a stochastic dynamical equation governing the motion of a charged particle in a collisional plasma is introduced in Section 2.3. This equation, which has the form of a generalized Langevin equation, is linked to the Fokker-Planck collision model of Rosenbluth et al. [1957] in Section 2.4. The components of this collision model, friction and diffusion coefficients, are then utilized to complete the formulation of the Langevin equation. These coefficients are specified in Section 2.5 for a Maxwellian plasma. 2.1 Soft-Target Radar Equation and Ambiguity Function In soft-target Bragg scattering, each infinitesimal volume dv = r2drdΩ of a radar field-of-view behaves like a hard-target with a radar cross section dv ´ dω 2pi σv(k, ω), where 10 σv(k, ω), by definition, is the soft-target RCS per unit volume per unit Doppler frequency ω/2pi, and k = −2korˆ denotes the relevant Bragg vector for a radar with a carrier wavenum- ber ko and associated wavelength λ. As a consequence, the hard-target radar equation generalizes for soft-target Bragg-scatter radars as [Milla and Kudeki , 2006] Pr(t) = ˆ dr dΩ dω 2pi G(rˆ)A(rˆ) (4pir)2 Pt(t− 2r c )σv(k, ω), (2.1) where Pt(t) ∝ |f(t)|2 is the transmitted power carried by a pulse waveform f(t), G(rˆ) and A(rˆ) are the gain and effective area of the radar antenna as a function of radial unit vector rˆ, r is the radar range, and c the speed of light. Furthermore, equation (2.1), which assumes an open radar bandwidth, can be recast for a match-filtered receiver output as Pr(t) EtKs = ˆ dr dΩ dω 2pi G2(rˆ)A(rˆ) 4pi |χ(t− 2rc , ω)|2 4pir2 σv(k, ω), (2.2) where Et is the total energy of the transmitted radar pulse, Ks denotes a system calibration constant including loss factors ignored in (2.1), and χ(τ, ω) ≡ 1 T ˆ dt ejωt f(t)f∗(t− τ) (2.3) has a magnitude known as the radar ambiguity function [e.g., Levanon, 1988]. In (2.3) the normalization constant T denotes the duration of the pulse waveform f(t) and (2.3) itself is effectively the normalized cross-correlation of the Doppler-shifted pulse f(t) ejωt with a delayed pulse echo proportional to f(t − τ) that would be expected from a point target located at a radar range R ≡ c τ/2. Hence, variable τ in (2.3) represents not only a time delay but also a corresponding radar range R. Let us denote by Sr(ω) the power spectrum of the signal detected by the radar receiver. Thus, applying Parseval’s theorem, we have that Pr = ´ dω 2pi Sr(ω). Assuming that σ(k, ω) varies slowly with the radar range r, and that the ambiguity function is almost flat within the bandwidth of the RCS spectrum (which is a valid approximation in the case of short- pulse radar applications), we can use equation (2.2) to obtain 11 Sr(ω) EtKs ≈ δR 4piR2 ˆ dΩ G(rˆ)A(rˆ) 4pi σv(k, ω), (2.4) where R ≡ ct/2 is the measured radar range and δR ≡ ˆ dr |χ(2 c (R− r), 0)|2 (2.5) is the effective range depth of the radar scattering volume. An important system parameter for soft-target radars is the backscatter aperture ABS , defined as ABS ≡ ˆ dΩ G(rˆ)A(rˆ) 4pi = λ2 16pi2 ˆ dΩG2(rˆ) (2.6) as a function of the two-way antenna gain G2(rˆ). The volumetric RCS spectrum of a quiescent ionosphere can be written as σv(k, ω) = 4pir 2 e 〈|ne(k, ω)|2〉, (2.7) where re is the classical electron radius and 〈|ne(k, ω)|2〉 denotes the space-time spectrum of electron density fluctuations ne to be formulated in the next section. 2.2 General Framework of IS Spectral Models In a plasma, electron and ion currents due to random thermal motions produce the impressed electron and ion density fluctuations nte and nti, which can be regarded as statistically independent random processes having the frequency spectra 〈|nte,i(k, ω)|2〉 = Ne,i ˆ dτe−jωτ 〈ejk·∆re,i〉 (2.8) expressed in terms of ambient densities Ne,i and random particle displacements ∆re,i as described in Kudeki and Milla [2010]. Electrostatic fields generated by the imbalance between nte and nti in turn drive a conduction current so that the total current, including the displacement current, vanishes in the electrostatic approximation. In the k-ω space, 12 the spectrum of total electron density fluctuations ne then exhibits a dependence on the spectra of nte and nti that is given by [e.g., Kudeki and Milla, 2006, 2010] 〈|ne(k, ω)|2〉 = |jωo + σi(k, ω)| 2〈|nte(k, ω)|2〉+ |σe(k, ω)|2〈|nti(k, ω)|2〉 |jωo + σe(k, ω) + σi(k, ω)|2 , (2.9) where σe(k, ω) and σi(k, ω) are the electron and ion conductivities, and o is the permit- tivity of free space. This expression, valid for stationary single-ion plasmas, is derived in Kudeki and Milla [2010] without making any assumption regarding the statistical distri- butions of electrons and ions. Expressions for the conductivity σs(k, ω) and the spectrum of thermal fluctuations 〈|nts(k, ω)|2〉 of each plasma species can be independently derived from plasma kinetic equations. Note that macroscopic electric interactions between different particle species have to be disregarded in the formulation of the kinetic equations as their effects have already been considered in the derivation of 〈|ne(k, ω)|2〉. However, if each species is in thermal equilibrium — i.e., if its velocity distribution is Maxwellian — it can be verified that σs(k, ω) and 〈|nts(k, ω)|2〉 are related by the fluctuation-dissipation or generalized Nyquist theorem [e.g., Callen and Welton, 1951; Kubo, 1966] ω2q2s k2 〈|nts(k, ω)|2〉 = 2KTsRe{σs(k, ω)}, (2.10) where qs and Ts are the charge and temperature of particles of type s, and K is the Boltzmann constant. Kudeki and Milla [2010] make use of this relation as follows: First, the spectrum of thermal fluctuations is derived from a microscopic model of particle dy- namics (where macroscopic “collective” interactions have been neglected), and then, the fluctuation-dissipation theorem is used to determine the real part of the conductivity of the corresponding particle species. Next, the imaginary part of σs(k, ω) is obtained by applying Kramers-Kronig relations [e.g., Clemmow and Dougherty , 1969; Chew , 1990]. As a result, it is found that 〈|nts(k, ω)|2〉 Ns = 2Re{Js(ω)} (2.11) 13 and σs(ω,k) jωo = 1− jωJs(ω) k2h2s , (2.12) where Ns is the mean species density, hs ≡ √ oKTs Nsq2s (2.13) is the corresponding Debye length, and Js(ω) ≡ ˆ ∞ 0 dτe−jωτ 〈ejk·∆rs〉 (2.14) is a Gordeyev integral, a one-sided Fourier transform of the characteristic function 〈ejk·∆rs〉 of random particle displacements ∆rs ≡ rs(t + τ) − rs(t) occurring over intervals τ in the absence of collective interactions [see Kudeki and Milla, 2010]. This characteristic function is also termed the single-particle ACF, since it can be interpreted as the normalized correlation of the signal scattered by a single particle exposed to a radar pulse. In this way, the problem of determining the incoherent scatter spectrum is reduced to the estimation of single-particle ACF’s and the evaluation of the corresponding Gordeyev integrals (in order to simplify the notation, subscript s is disregarded in the rest of this chapter). 2.3 Single-Particle ACF and Langevin Equation The single-particle ACF 〈ejk·∆r〉 for a species in a given plasma can be calculated directly from 〈ejk·∆r〉 = ˆ d(∆r) ejk·∆rf(∆r, τ), (2.15) where f(∆r, τ) is the pdf of the particle displacement ∆r at a time delay τ . This is a straightforward calculation as long as f(∆r, τ) is known either analytically or numerically. The standard procedure for determining f(∆r, τ) starts with solving the Boltzmann ki- netic equation for f(r, t) with a proper collision operator, e.g., the Fokker-Planck kinetic equation of Rosenbluth et al. [1957] for Coulomb collisions, subject to an initial condition 14 [e.g., Holod et al., 2005] f(r, 0) = δ(r). (2.16) Although, analytical solutions of simplified versions of the Fokker-Planck kinetic equation are available [e.g., Chandrasekhar , 1943; Dougherty , 1964], determining f(∆r, τ) would be a daunting task when the full equation is considered. An alternative approach to calculating 〈ejk·∆r〉 in the case of Coulomb collisions in- volves producing, as the result of some simulation procedure, suitable sets of particle dis- placement data ∆r for the same stochastic process described by the Fokker-Planck kinetic equation. Assuming the process to be ergodic, a sufficiently large set of simulated samples of particle trajectories r(t) can be used to compute 〈ejk·∆r〉 as well as any other statistical function of displacements ∆r over time delays τ . Random trajectory simulations of course require the availability of a stochastic equation describing how the particle velocities v(t) ≡ dr dt (2.17) may evolve under the influence of Coulomb collisions. Restricting v(t) to be a Markovian random process constrains the stochastic evolution equation of v(t) by very strict self-consistency conditions discussed by Gillespie [1996a, b] to acquire the form of a Langevin equation dv(t) dt = A(v, t) + C¯(v, t)W(t), (2.18) where vector A(v, t) and matrix C¯(v, t) consist of arbitrary smooth functions of arguments v and t. Above, W(t) is a random vector having statistically independent Gaussian white noise components. A more natural way of expressing the Langevin equation is to cast it in an update form, namely v(t+ ∆t) = v(t) + A(v, t) ∆t+ C¯(v, t) ∆t1/2U(t), (2.19) where ∆t is an infinitesimal update interval and U(t) is a vector composed of independent 15 zero-mean Gaussian random variables with unity variance, i.e., Ui(t) ≡ N (0, 1), (2.20) where N (µ, σ2) denotes the normal random variable with mean µ and variance σ2. The formal equivalency of (2.18) and (2.19) requires the i-th component of W(t) to be defined as Wi(t) = lim ∆t→0 (∆t)−1/2N (0, 1) = lim ∆t→0 N (0, 1/∆t), (2.21) compatible with the requirement that 〈Wi(t+ τ)Wi(t)〉 = δ(τ). Note that the Langevin equation (2.18) describing a Markovian process has the form of Newton’s second law of motion, with the terms on the right representing forces per unit mass that model the fluctuating Coulomb fields exerted on the particles in a plasma and that will be regarded as the collision forces. The first term is a deterministic friction force that causes particle deceleration, while the second term is a stochastic diffusion force that is responsible for the random walk of the simulated particles in space. Considering the Lorentz force on a charged particle in a magnetized plasma with a constant magnetic field B, and not violating the strict format of (2.18), we can modify the equation by adding a term qmv(t)×B to its right hand side. We would then have an additional term q mv(t)×B ∆t on the right-hand side of the update equation (2.19) as well. Another relevant fact is that a special type of Markov process characterized by a linear A(v, t) = −βv and a constant matrix C¯ = D1/2I¯, independent of v and t, is known as Brownian motion process [e.g., Uhlenbeck and Ornstein, 1930; Chandrasekhar , 1943], which is often invoked in simplified models of collisional plasmas [e.g., Dougherty , 1964; Woodman, 1967; Holod et al., 2005]. In these models, friction and diffusion coefficients, β and D, are constrained to be related by D = 2KT m β (2.22) for a plasma in thermal equilibrium. 16 2.4 The Link between Fokker-Planck and Langevin Equations In return for having restricted v(t) to the space of Markovian processes, we have gained a stochastic evolution equation (2.18) with a plausible Newtonian interpretation and with the potential of taking us beyond Brownian-motion-based collision models. Letting v(t) to be Markovian turns out to have a second consequence of importance in this work: namely, the evolution of probability density f(v, t) of a random variable v(t) is known to be governed, when v(t) is Markovian, by the Fokker-Planck kinetic equation having a “friction vector” and “diffusion tensor” 〈 ∆v ∆t 〉 c = A(v, t) (2.23) and 〈 ∆v∆vT ∆t 〉 c = C¯(v, t)C¯T(v, t), (2.24) respectively, specified in terms of the input functions of the Langevin equation. This intimate link between the Langevin equation and the Fokker-Planck kinetic equation — in describing Markovian processes from two different but mutually compatible perspectives — was first pointed out by Chandrasekhar [1943] and more recently by Gillespie [1996b]. This relationship is analogous to a similar relation between the BGK kinetic equation [Bhatnagar et al., 1954] for charged particles colliding with neutrals and a particle dynamics formalism in which this type of collisions is modeled as a discrete Poisson process, a model that is discussed further by Milla and Kudeki [2009] (see Appendix A). Since the Fokker-Planck friction vector and diffusion tensor for equilibrium plasmas with Coulomb interactions have already been worked out by Rosenbluth et al. [1957] as 〈 ∆v ∆t 〉 c = −β(v)v (2.25) and 〈 ∆v∆vT ∆t 〉 c = D⊥(v) 2 I¯ + ( D‖(v)− D⊥(v) 2 ) vvT v2 , (2.26) 17 in terms of scalar functions β(v), D‖(v), andD⊥(v) to be specified below, it follows that the Langevin update equation to be used in Coulomb collision particle trajectory simulations can be written as v(t+ ∆t) = v(t) + q m v(t)×B ∆t − β(v)∆tv(t) + √ D‖(v)∆t U1 vˆ‖ + √ D⊥(v) ∆t 2 (U2 vˆ⊥1 + U3 vˆ⊥2) , (2.27) including a DC magnetic field term, with vˆ‖(t), vˆ⊥1(t), and vˆ⊥2(t) denoting an orthogonal set of unit vectors parallel and perpendicular to the particle trajectory, and Ui denoting zero-mean Gaussian random variables with unit variance. 2.5 Spitzer Friction and Diffusion Coefficients Expressions for the Fokker-Planck friction vector and diffusion tensor for a particle moving in a magnetized plasma have been derived by different authors [e.g., Rostoker and Rosenbluth, 1960; Ichimaru and Rosenbluth, 1970]. However, when the plasma is “weakly” magnetized, i.e., if the particle Debye length is smaller than the mean gyroradius, magnetic field effects in the friction vector and diffusion tensor can be neglected [Rostoker and Rosenbluth, 1960] and the friction coefficient β(v) and velocity space diffusion coefficients D‖(v) andD⊥(v) introduced above can be specified in (2.27) as formulated by the standard collision model of Rosenbluth et al. [1957]. For a Maxwellian plasma, these coefficients take the forms given by Spitzer [Spitzer , 1962]. Specifically, if a particle of type s colliding against a background of particles of type s′ is considered, the Spitzer coefficients take the forms [Callen, 2006] β(v) = ∑ s′ (1 + ms ms′ )Ns′Γs/s′ 1 C2s′ ψ(zs′) v , (2.28) D||(v) = ∑ s′ 2Ns′Γs/s′ ψ(zs′) v , (2.29) D⊥(v) = ∑ s′ 2Ns′Γs/s′ φ(zs′)− ψ(zs′) v , (2.30) 18 where Cs = √ KTs/ms is the thermal speed of particles with mass ms, electric charge qs, and equilibrium temperature Ts. Additionally, Γs/s′ ≡ q2sq 2 s′ 4pi2om 2 s ln Λs/s′ (2.31) in terms of “Coulomb logarithm” ln Λs/s′ (see below), while zs ≡ v√ 2Cs (2.32) in terms of particle speed v, and ψ(z) ≡ φ(z)− zφ ′(z) 2z2 , (2.33) where φ(z) ≡ 2√ pi ˆ z 0 e−x 2 dx (2.34) is the well-known error function. Finally, defining a plasma Debye length hD ≡ (∑ s 1 h2s )−1/2 = (∑ s Nsq 2 s oKTs )−1/2 , (2.35) and a minimum impact parameter bmin ≡ qsqs′ 12piomss′C 2 ss′ , (2.36) where mss′ ≡ msmx′ms+ms′ is a reduced mass and C 2 ss′ ≡ C2s + C2s′ , we have Λs/s′ ≡ hD bmin , (2.37) which is known as the plasma parameter and is proportional to the number of particles of type s inside a Debye cube. See Appendix B for a discussion regarding the application of Spitzer friction and diffusion coefficients to the case of non-isothermal plasmas. 19 2.6 Summary The general framework of the incoherent scatter theory formulates the electron den- sity spectrum in terms of electron and ion Gordeyev integrals, functions that depend on the statistics of individual particle trajectories. Since collisions in the plasma modify these trajectories, they have an impact on the statistics of the particle displacements and also on the corresponding Gordeyev integrals. It is through the estimation of these functions that the effect of Coulomb collisions is considered in our incoherent scatter spectrum model. Although the statistics of particle trajectories, i.e., velocity and displacement density dis- tributions, can be obtained by numerically solving the Boltzmann kinetic equation with the Fokker-Planck collision operator, we chose to utilize a different approach based on the simulation of plasma particle trajectories using Langevin update equations. In the next chapter, we describe our computer simulations which provide us with sample electron and ion trajectories in magnetized and collisional F-region plasmas. 20 CHAPTER 3 MONTE CARLO COMPUTATION OF PARTICLE TRAJECTORIES AND SINGLE-PARTICLE ACF’S The Monte Carlo procedure for the simulation of charged particle trajectories in a plasma and the subsequent estimation of single-particle ACF’s is detailed in this chapter. In the simulations, ion and electron trajectories are computed using a stochastic update equation of the Langevin type where friction and diffusion forces account for the effects of Coulomb collisions on the particle motion. The computer program developed for the simulations is outlined in Section 3.1. The outcomes of the simulations, series of particle positions, are then utilized in the computation of single-particle ACF’s and associated Gordeyev integrals according to the algorithm described in Section 3.2. Velocity and displacement distributions obtained from the simulations are analyzed in detail in Section 3.3, where samples of ion and electron ACF’s are displayed and discussed. The following presentation follows closely the description of the procedure reported by Milla and Kudeki [2010]. 3.1 Computer Simulations of Plasma Particle Trajectories Consider the motion of a particle governed by the Langevin update equation (2.27). Denoting by v[n] the particle velocity vector at a discrete time n∆t, the updated velocity 21 Table 3.1 Typical plasma parameters for an equatorial F-region ionosphere. Parameter Value Plasma density Ne 1012 m−3 Magnetic field Bo 25µT Electron temperature Te 1000 K Ion temperature Ti 1000 K Ion composition O+ after a sufficiently small time interval ∆t can be determined by using v[n+ 1] = v[n] + q m v[n]×B ∆t − β(v)∆tv[n] + √ D‖(v)∆t U1[n] vˆ‖ + √ D⊥(v) ∆t 2 (U2[n] vˆ⊥1 + U3[n] vˆ⊥2) , (3.1) where m and q denote the mass and the electric charge of the simulated particle (which can be either an electron or an ion). The different terms in (3.1) on the right include the changes in the velocity produced by the magnetic and collisional forces discussed in Chapter 2. In addition, the particle position r[n] can be calculated from the velocities v[n] using r[n+ 1] = r[n] + v[n+ 1] + v[n] 2 ∆t, (3.2) another update equation that can be obtained by approximating the integral of the velocity vector over a time interval ∆t using the trapezoidal rule. We have developed a computer program for the simulation of F-region particle tra- jectories governed by (3.1) and (3.2). In the simulations, a homogeneous plasma in the presence of a uniform magnetic field B ≡ Bo zˆ is considered. The plasma has electron density Ne and electron and ion temperatures Te and Ti. The values of Bo, Ne, Te, and Ti, which are typical of the equatorial F-region ionosphere, are listed in Table 3.1. For completeness, a list of some representative plasma parameters is given in Table 3.2. The F-region simulation runs as follows. First, we randomly pick some initial velocity 22 T ab le 3. 2 R ep re se nt at iv e ph ys ic al pa ra m et er s of a pl as m a. D efi ni ti on s an d sa m pl e va lu es ar e pr ov id ed . T he va lu es w er e co m pu te d fo r an O + pl as m a w it h pa ra m et er s gi ve n in T ab le 3. 1. N ot e th at th e pl as m a an d gy ro fr eq ue nc y fo rm ul as ar e gi ve n fo r an gu la r fr eq ue nc ie s, bu t th e sa m pl e va lu es ar e gi ve n in H z. P ar am et er D efi ni ti on E le ct ro n/ Io n D efi ni ti on P la sm a T he rm al sp ee d C s = √ K T s /m s 12 3 .1 k m / s 71 8 .3 m /s — — P la sm a fr eq ue nc y ω p s = √ N s q2 s / ( o m s ) 8 .9 8 M H z 52 .3 8 k H z ω p = ( ∑ s ω 2 p s ) 1/2 8. 98 M H z D eb ye le ng th h s = C s /ω p s 2. 18 m m 2. 18 m m h D = ( ∑ s 1/ h 2 s ) −1/2 1 .5 4 m m G yr of re qu en cy Ω s = q s B o /m s 69 9 .8 k H z 23 .8 2 H z Ω p = ( ∑ s Ω 2 s ) 1/2 69 9. 8 k H z M ea n gy ro ra di us ρ s = √ 2 C s /Ω s 39 .6 m m 6 .7 9 m — — 23 v[0] and set the starting position of the particle to the origin (i.e., r[0] = 0). Next, the values of friction and diffusion coefficients for the given particle speed are calculated. The velocity increments resulting from the action of each of the simulated forces (magnetic, friction, and diffusion) in a time interval ∆t are then computed. For the calculation of the diffusive forces, we have to randomly pick three normal numbers. These values are computed using the Ziggurat method for random generation of Gaussian variables [Marsaglia and Tsang , 2000]. Once all the velocity increments are calculated, we proceed to compute the new value of the velocity vector using (3.1), and subsequently, the position of the particle is updated using (3.2). The algorithm is then repeated in a loop. The resulting series of velocities and positions are stored in sequences of length M (typically equal to 217 samples). About 104 of these sequences are generated, which provides more than 109 simulated velocities and positions for a single particle in a given plasma configuration. Note that the starting velocity and position of each new sequence is set equal to the last velocity and position of the previous sequence. Only the initial velocity of the first sequence is randomly chosen; the remaining samples of the series are calculated using the update velocity equation (3.1). A sample trajectory of an electron moving in an O+ plasma with densityNe = 1012 m−3 and temperatures Te = Ti = 1000 K is presented in Figure 3.1. The intensity of the magnetic field was taken as 25µT. On the top, the left panel shows the particle trajectory in 3D space, while the panel on the right is the projection of the trajectory on the xy- plane (i.e., the plane perpendicular to B). On the bottom, the displacement of the electron in the direction parallel to B is displayed as a function of time. As shown, the electron moves approximately in a spiral path that is randomized by the action of collisions. As the electron accelerates and decelerates at random rates, not only does the guiding center of its trajectory randomly drift in space, but the diameter of the particle orbits also changes as a function of time. Because of this, there are time intervals in which the electron orbits have smaller or larger radii than the mean gyroradius ρe. For instance, the trajectory presented in Figure 3.1 corresponds to a period when the orbit radius is smaller than ρe (which, in this example, is approximately 40 mm). We can also see that, in the plane perpendicular 24 −0.1 0 0.1 −0.1 0 0.1 −5 0 5 10 15 x-axis [m] Time Interval: 0-400 µs y -axis [m] z -a x is [m ] −0.1 −0.05 0 0.05 0.1 −0.1 −0.05 0 0.05 0.1 x-axis [m] y -a x is [m ] Displacement ⊥ to B 0 50 100 150 200 250 300 350 400 −5 0 5 10 15 Time [µs] z -a x is [m ] Displacement ‖ to B Figure 3.1 Sample trajectory of an electron moving in an O+ plasma with density Ne = 1012 m−3 and temperatures Te = Ti = 1000 K. The presence of a uniform magnetic field B parallel to the z-axis was considered in the simulation. The intensity of the field was assumed to be 25µT. On the top, the left panel shows the particle trajectory in 3D space, while the panel on the right is the projection of the trajectory on the xy-plane (i.e., the plane perpendicular to B). On the bottom, the displacement of the electron in the direction parallel to B is displayed as a function of time. The trajectory was sampled every ∆t = 0.1µs. A total of 4 × 103 samples are displayed corresponding to a time interval of 400µs. In all three panels, the red dots show the starting location of the particle (the origin). In addition, the red curve depicts the first 1.4µs of the simulated trajectory (about one gyroperiod). 25 to B, the particle never returns to the same position after a gyroperiod. Additionally, note that there is a large difference between the distances covered by the electron in the directions parallel and perpendicular to B. These characteristics of the simulated electron trajectories have implications for the shape of the computed single-particle ACF’s and their corresponding Gordeyev integrals presented later in this chapter. The simulation procedure just outlined was motivated by the earlier work of Sulzer and González [1999]. Both simulation studies make use of Spitzer friction and diffusion coefficients in order to include the effects of Coulomb collisions in IS spectral models. However, the equations of motion employed by Sulzer and González [1999] differed from our Langevin-based 3D update procedure. Sulzer and González considered the effect of Coulomb collisions on particle displacements only in the direction of the ambient field B, neglecting the diffusion and random walk effects across the field lines. That limits the applicability of their results to magnetic aspect angles larger than about 0.1◦, in relation to the 50 MHz Jicamarca radar observations. Our simulation results, by contrast, enable ACF and Gordeyev integral calculations for all aspect angles, including the direction of exact perpendicularity to the geomagnetic field B. Furthermore, our Langevin-based update procedure does not suffer numerical instability issues when used with sufficiently small update intervals ∆t. 3.2 Estimation of the Single-Particle ACF’s and Gordeyev Integrals In this section, we will describe our processing procedures of the particle trajectory data for estimating the single-particle ACF’s and the corresponding Gordeyev integrals. Consider a long sequence of N simulated particle positions uniformly sampled in time. For a given radar Bragg vector k, the unbiased estimator for the single-particle ACF at a discrete time delay m∆t is given by R̂[m] ≡ ̂〈ejk·∆r[m]〉 = 1 N −mcrr[m], 0 ≤ m ≤ N − 1, (3.3) 26 where crr[m] = N−m−1∑ n=0 exp (jk · (r[n+m]− r[n])) = N−1∑ n=m exp (jk · (r[n]− r[n−m])) (3.4) is the discrete correlation of a sequence of samples ejk·r[n]. In general, discrete correlations can be computed efficiently using the FFT technique. In our case, however, we cannot directly apply this procedure because the entire sequence of N particle positions is very large and it cannot be fully allocated in the RAM of a computer (for instance, 109 sample positions stored in double float format would require about 24 gigabytes of RAM). Since we are just interested in the first M samples of this correlation (about 105 samples), let us divide the full set of positions into L segments of length M , such that N = LM . We can then re-express crr[m] as crr[m] = L−1∑ l=0 2M−1∑ n=0 gl[n]h ∗ l [(n−m)2M ], 0 ≤ m ≤M − 1, (3.5) where (a)b denotes the modulo operation (amod b), and functions gl and hl are gl[n] ≡  exp(jk · r[n+ lM ]), 0 ≤ n ≤M − 1, 0, M ≤ n ≤ 2M − 1, (3.6) and hl[n] ≡  gl[n], 0 ≤ n ≤M − 1, gl−1[n−M ], M ≤ n ≤ 2M − 1, (3.7) respectively (note that h0 = g0). Using equation (3.5), we can compute the correlation crr[m] in an iterative way reducing significantly the amount of required computer mem- ory. The inner summation in (3.5) is the circular cross-correlation of functions gl and hl, calculation of which is performed using FFT’s. Estimates Rˆ[m] of the single-particle ACF 〈ejk·∆r〉 are then obtained by dividing crr[m] by LM −m. This means that the statistical variance of our estimates increases with m, since for larger m fewer samples are used to estimate the ACF. Assuming that at every 27 M samples the simulated positions are independent random variables, we find that in the worst-case scenario (i.e., when m = M − 1), the variance of our estimates will be at most inversely proportional to L − 1. Thus, in order to secure small estimation errors, large values of L have to be considered. In our calculations we used L = 104. The procedure outlined above is used to compute 〈ejk·∆r〉 for different values of Bragg vector k. In our calculations, we define k = 2pi λB (cos(α) xˆ+ sin(α) zˆ) , (3.8) where λB is the Bragg wavelength (e.g., λB = 3 m in the case of Jicamarca), and α is the magnetic aspect angle that is the complement of the angle between k and the magnetic field vector B = Bo zˆ. Notice that to estimate 〈ejk·∆r〉 for different aspect angles, we do not need to generate a new set of particle positions. We took advantage of this and performed many of these calculations in parallel on a computer. Our definition of vector k is somewhat arbitrary. Notice that, for the direction perpendicular to B, we could have chosen k to be 2piλB yˆ instead of 2pi λB xˆ. Therefore, for α = 0◦, there are two alternative ways of computing the single-particle ACF that provide two sets of almost statistically independent estimates of 〈ejk·∆r〉. Averaging these two sets, we could have reduced the statistical errors of our estimates. In the future, we will take advantage of this in our calculations. Our next task is the computation of the Gordeyev integral J(ω) from discrete estimates of 〈ejk·∆r〉. Since our samples are uniformly distributed in time, we could have simply taken the FFT of the ACF estimates Rˆ[m] to perform the Gordeyev integral calculations, but then samples of J(ω) would have been restricted to a discrete set of frequencies. We instead used a chirp-Z transform algorithm described by Li et al. [1991] in our Gordeyev integral computations, which allows the evaluation of the transform over any desirable range of frequencies ω. Given the set of frequencies ωk = ωo + k∆ω for 0 ≤ k ≤ K − 1, we can approximate 28 the Gordeyev integral (2.14) using the following summation J(ωk) ≈ ∆t M−1∑ m=0 w[m] Rˆ[m] e−jωkm∆t = ∆t M−1∑ m=0 w[m] Rˆ[m] e−jωom∆t e−jkm∆ω∆t, (3.9) where, in order to improve the accuracy of the approximation, ACF estimates Rˆ[m] are weighted by composite quadrature coefficients w[m]. In our calculations, we use composite Simpson’s rule coefficients which are given by w[m] =  1 3 , m = 0,M, 4 3 , m = 1, 3, . . . ,M − 1, 2 3 , m = 2, 4, . . . ,M − 2, (3.10) where M is assumed to be an even number. To derive the integration algorithm, let us use the identity km = 1 2 (k2 +m2 − (k −m)2) (3.11) to re-express equation (3.9) as J(ωk) ≈ ∆t M−1∑ m=0 w[m] Rˆ[m] e−jnωo∆tW k 2/2Wm 2/2W−(k−m) 2/2 (3.12) where W = e−j∆ω∆t. Note that the previous equation is a convolution sum; then, if we define x[m] ≡  w[m] Rˆ[m] e−jmωo∆tWm2/2, 0 ≤ m ≤M − 1, 0, M ≤ m ≤M +K − 1, (3.13) and y[m] ≡  W−m2/2, 0 ≤ m ≤ K − 1, W−(M+K−m)2/2, K ≤ m ≤M +K − 1, (3.14) 29 Table 3.3 Plasma parameters spanned by the library of single-electron ACF’s. Plasma parameter Initial value Final value Step log10Ne 11 13 0.5 Te 600 K 3000 K 200 K Ti 600 K 2000 K 200 K Bo 20µT 30µT 5µT we can rewrite equation (3.12) in the following form J(ωk) ≈ ∆tW k2/2 M+K−1∑ m=0 x[m] y[(k −m)M+K ], 0 ≤ k ≤ K − 1, where the summation is the circular convolution of sequences x[m] and y[m], calculations that can be efficiently performed using FFT’s. In addition, note that in order to avoid aliasing artifacts, we need ωK∆t < pi. In the above discussion, we made no distinction between electrons and ions since identi- cal procedures are used to compute the Gordeyev integrals for each species. The estimated electron and ion Gordeyev integrals are then used to compute theoretical incoherent scatter spectra. As a result of our simulation studies, we found that the ion Gordeyev integrals can be well represented by analytical means, as discussed in the next section of this chapter. However, closed-formed expressions could not be found for the electron Gordeyev inte- grals; therefore, a numerical electron ACF library had to be constructed for a wide range of plasma parameters and magnetic aspect angles, as needed by Jicamarca applications. The plasma parameters that this library spans are given in Table 3.3. The construction of the library was computationally demanding and was carried out using a cluster of computers. The Turing Cluster, maintained by the Computational Science and Engineering Program at the University of Illinois, was used for this purpose. The Turing system consists of 768 Apple Xserves, each with two 2 GHz G5 processors and 4 GB of RAM. In the cluster, hundreds of simulations run simultaneously, which allowed us to build the library in less than three weeks. The same task would have taken many 30 months (probably more than two years) using a single desktop computer. 3.3 Statistics of the Ion and Electron Trajectories If a plasma is in thermal equilibrium, i.e., if Te = Ti, one can show analytically that the steady-state solution of the Fokker-Planck equation is given by the Maxwell-Boltzmann velocity distribution [Montgomery and Tidman, 1964] f(v) = 1 (2piC2)3/2 e− 1 2C2 (v2x+v2y+v2z), (3.15) where vx, vy, and vz are the components of the particle velocity vector, and C = √ KT/m is the corresponding thermal speed. Notice that this pdf can be written as the product of three independent Gaussians, one for each of the velocity components. Since the Fokker- Planck and the Langevin equations are alternative representations of the same Markov process, we expect the distributions of our simulated velocities to be Gaussians. Velocity distributions resulting from independent ion and electron simulations in an O+ plasma are presented in Figure 3.2. In the simulations, the plasma was considered to be in thermal equilibrium with temperatures Te = Ti = 1000 K (the rest of the simulation plasma param- eters are given in Table 3.1). The distributions were built using more than 109 sampled velocities. In each panel, the expected Gaussian pdf is plotted on top of our simulation results. As we can see, the computed distributions match exactly the Gaussian curves in these two examples. We repeated the same test for different plasma parameters and found that the agreement was excellent in all cases. Because of these results, we have the confidence that our simulation procedure is working properly. As we mentioned in Chapter 2, the single-particle ACF is a characteristic function; therefore, its shape is determined by the time evolution of the pdf f(∆r, τ) of particle displacements ∆r. If the plasma is magnetized, particles are forced to diffuse slowly in the plane perpendicular to B, making the pdf f(∆r, τ) narrower in that plane than in the parallel direction (for a given τ). In order to analyze the behavior of f(∆r, τ), we have computed from the simulated trajectories the distributions of ion and electron 31 −3 −2 −1 0 1 2 3 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Normalized velocity (v/Ci) Io n v el o ci ty d is tr ib u ti o n vx pdf vy pdf vz pdf Gaussian −3 −2 −1 0 1 2 3 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Normalized velocity (v/Ce) E le ct ro n v el o ci ty d is tr ib u ti o n vx pdf vy pdf vz pdf Gaussian (a) (b) Figure 3.2 Probability distributions of (a) ion and (b) electron velocities resulting from the simulations. In each plot, the pdf’s of the velocity components (displayed in colors) are compared to a Gaussian distribution (black line). Notice that the velocity axes are normal- ized by the thermal speeds of the corresponding particle species. The plasma parameters considered in the simulations are given in Table 3.1. displacements in the directions parallel and perpendicular to B. Since B = Bo zˆ, the distributions of the parallel displacement correspond to f(∆z, τ). On the other hand, the distributions of the perpendicular displacement were computed by averaging f(∆x, τ) and f(∆y, τ). In addition, the variance and covariance of the components of the displacement vector ∆r were computed from the simulated trajectories. These mean-square quantities are defined as 〈∆ri ∆rj〉 = 〈(ri(t+ τ)− ri(t)) (rj(t+ τ)− rj(t))〉, (3.16) where i and j denote the Cartesian coordinates. These quantities were estimated following a procedure similar to the one used to compute 〈ejk·∆r〉. Before presenting our statistics of the simulated ion and electron displacements, notice that expression (2.15) for the ACF 〈ejk·∆r〉 is simply the Fourier transform of f(∆r, τ). In principle, we could have computed 〈ejk·∆r〉 using the Fourier transforms of the dis- placement distributions obtained in the simulations. This would have required, however, a significant increase in the number of simulated particle positions in order to reduce the statistical estimation errors to an acceptably small level. Since we are interested in eval- uating 〈ejk·∆r〉 only for some discrete values of k, we consider the procedure outlined in 32 0 2 4 6 8 10−2 0 2 0 0.2 0.4 Time delay [ms] Ion displacement distributions ⊥ to B ∆r⊥(τ )/σ⊥(τ ) ∆ r ⊥ p d f 0 2 4 6 8 10−2 0 2 0 0.2 0.4 Time delay [ms] Ion displacement distributions ‖ to B ∆r‖(τ )/σ‖(τ ) ∆ r ‖ p d f −3 −2 −1 0 1 2 3 0 0.1 0.2 0.3 0.4 ∆r⊥(τ )/σ⊥(τ ) ∆ r ⊥ p d f Ion displacement distributions ⊥ to B τ = 0.0ms τ = 2.0ms τ = 4.0ms τ = 6.0ms τ = 8.0ms τ = 10.0ms Gaussian −3 −2 −1 0 1 2 3 0 0.1 0.2 0.3 0.4 ∆r‖(τ )/σ‖(τ ) ∆ r ‖ p d f Ion displacement distributions ‖ to B τ = 0.0ms τ = 2.0ms τ = 4.0ms τ = 6.0ms τ = 8.0ms τ = 10.0ms Gaussian Figure 3.3 Probability distributions of the displacements of a simulated ion in the direc- tions perpendicular (top panels) and parallel (bottom panels) to the magnetic field. On the left, the displacement pdf’s are displayed as functions of time delay τ . On the right, sample cuts of the pdf’s are compared to a Gaussian distribution. Note that all distribu- tions at all time delays are normalized to unit variance. The displacement axis of each distribution at every delay τ is scaled with the corresponding standard deviation of the simulated displacements. the previous two sections to be more accurate, involving fewer computations. 3.3.1 Statistics of the Ion Displacements In Figure 3.3, we show the distributions of ion displacements in the directions perpen- dicular and parallel to the magnetic field. Note that, at every delay τ , the distributions have been normalized to unit variance by scaling the displacement axis of each distribu- tion with the corresponding standard deviation of the simulated displacements.1 On the 1Note that in the limit τ → 0, we can write lim τ→0 ∆ri(τ) σi(τ) = lim τ→0 ri(t+ τ)− ri(t) Cτ = vi(t)/C. Therefore, the distribution of ∆ri(τ)/σi(τ) at τ = 0 is equal to the distribution of the corresponding velocity component. 33 left panels, the distributions are displayed as functions of τ , while, on the right panels, sample cuts of these distributions are compared to a Gaussian pdf. For the time interval considered in these plots, we can see that the shapes of the distributions do not change with time delay and, also, that they match almost perfectly to the Gaussian curves. We have observed that for τ up to and exceeding 10 ms, the distribution of the displacement in the direction perpendicular to B remains Gaussian in shape. We also found that at time delays of the order of hundreds of milliseconds, the parallel distribution eventually becomes spikier than a Gaussian. However, by that time the single-ion ACF is negligibly small — typical correlation times of the ion ACF for λB = 3 m are of the order of 1 ms — and for that reason, ion displacements can be regarded as Gaussian random variables for all practical purposes. Additionally, we verified that the components of the vector displacement (i.e., ∆rx, ∆ry, and ∆rz) are mutually uncorrelated. The simulation results presented here were computed for the plasma configuration of Table 3.1. The analysis was repeated for other possible ionospheric plasmas and similar results were observed. For the case of uncorrelated and jointly Gaussian ∆r components, it is known that [e.g., Kudeki and Milla, 2010] the single-particle ACF takes the following form 〈ejk·∆r〉 = e− 12k2 sin2 α〈∆r2‖〉 × e− 12k2 cos2 α〈∆r2⊥〉, (3.17) where, assuming a Brownian-motion process with distinct friction coefficients ν‖ and ν⊥ in the directions parallel and perpendicular to B, the mean square displacements will vary as 〈∆r2‖〉 = 2C2 ν2‖ ( ν‖τ − 1 + e−ν‖τ ) (3.18) and 〈∆r2⊥〉 = 2C2 ν2⊥ + Ω2 ( cos(2γ) + ν⊥τ − e−ν⊥τ cos(Ωτ − 2γ) ) (3.19) in which γ ≡ tan−1(ν⊥/Ω), and C ≡ √ KT/m and Ω ≡ qB/m are, respectively, the thermal speed and gyrofrequency of the simulated particles. In the literature, friction coefficients ν‖ and ν⊥ are also referred to as “collision frequencies.” However, in the context of Coulomb collisions, they should be interpreted not as the rate of collisions between 34 0 1 2 3 4 50 5 10 〈∆ r 2 ⊥ 〉[ m 2 ] Variance of ion displacement ⊥ to B 〈|∆x|2〉 〈|∆y|2〉 Best fit ν⊥ = 5.88 Hz 0 1 2 3 4 50 5 10 Time delay [ms] 〈∆ r 2 ‖〉 [m 2 ] Variance of ion displacement ‖ to B 〈|∆z |2〉 Best fit ν‖ = 5.13 Hz 0 100 200 300 400 5000 50 100 150 〈∆ r 2 ⊥ 〉[ m 2 ] Variance of ion displacement ⊥ to B 〈|∆x|2〉 〈|∆y|2〉 Best fit ν⊥ = 5.88 Hz 0 100 200 300 400 5000 2 4 6 x 104 Time delay [ms] 〈∆ r 2 ‖〉 [m 2 ] Variance of ion displacement ‖ to B 〈|∆z |2〉 Best fit ν‖ = 5.13 Hz Figure 3.4 Variances of the displacements of a simulated ion in the directions perpen- dicular (top panels) and parallel (bottom panels) to the magnetic field. The simulation results (color lines) are displayed for two time intervals: 5 ms (left panels) and 500 ms (right panels). The dashed lines correspond to the fits of the results to the theoretical expressions for 〈∆r2⊥〉 and 〈∆r2‖〉 of the Brownian-motion model. particles but instead as the rate of loss of momentum (deceleration) experienced by a particle due to its interaction with the Coulomb fields generated by the swarm of electrons and ions surrounding it. To test the viability of a Brownian-motion model to represent the simulated ion data, we fitted the expressions (3.18) and (3.19) to the variances of ion displacements obtained in our simulations and compared the best-fit parameters ν‖ and ν⊥ to the Spitzer collision frequency2 for ion-ion interactions, which for the case of a single-ion plasma is given by [Callen, 2006] νi/i = Ne e 4 ln Λi 12pi3/2 2om 2 i C 3 i , (3.20) where Λi ≡ 12piNeh2ihD is the ion plasma parameter, and hi and hD are the ion and plasma Debye lengths defined in Table 3.2. An example of our fitting results is presented in Figure 3.4. In general, we found that ν⊥ ≈ νi/i. For instance, for the case shown in Figure 3.4, the best-fit collision frequency is ν⊥ = 5.88 Hz, while the Spitzer collision frequency is νi/i = 5.94 Hz. We also found that 2Spitzer collision frequency νs/s′ is the Maxwellian-averaged momentum relaxation rate of a particle of type s due to its interaction with a background of particles of type s′. The inverse of νs/s′ can be interpreted as the time interval over which the particle velocity vector rotates by about 90◦, an accumulated effect due to many small Coulomb deflections known as Spitzer collisions. 35 0 0.5 1 1.5 2 2.5 30 0.2 0.4 0.6 0.8 1 Time [ms] (a) Ion ACF (λB = 3 m) α = 0.0◦ α = 0.5◦ α = 1.0◦ α = 90.0◦ Brownian 0 0.05 0.1 0.15 0.2 0.25 0.30 0.2 0.4 0.6 0.8 1 Time [ms] (b) Ion ACF (λB = 0.3 m) α = 0.0◦ α = 0.5◦ α = 1.0◦ α = 90.0◦ Brownian Figure 3.5 Simulated single-ion ACF’s at different magnetic aspect angles α for two radar Bragg wavelengths: (a) λB = 3 m and (b) λB = 0.3 m. The simulation results (color lines) are compared to theoretical ACF’s computed using expression (3.22) of the Brownian- motion approximation (dashed lines). Note that there is effectively no dependence on aspect angle α. the best-fit parallel collision frequency ν‖ is smaller than the Spitzer frequency by a factor of ∼ 1.15, i.e., ν‖ < νi/i. Finally, we note in the same figure that the variances of ∆r‖ and ∆r⊥ are very similar for τ < 5 ms. This is expected because for τ in that interval ν‖τ  1 and ν⊥τ < Ωiτ  1, in which case (3.18) and (3.19) indicate that 〈∆r2‖〉 ≈ 〈∆r2⊥〉 ≈ C2i τ2. (3.21) Using this approximation in expression (3.17), we find that the single-ion ACF simplifies to 〈ejk·∆ri〉 ≈ e− 12k2C2i τ2 , (3.22) which is a well-known result for collisionless and non-magnetized plasmas that fits well the simulated single-ion ACF’s for different magnetic aspect angles and Bragg wavelengths as shown in Figure 3.5. Evidently, in a magnetized F-region plasma with Coulomb collisions, the oxygen ions diffuse along and across the magnetic field lines like Brownian-motion particles having 36 isotropic friction coefficients ν⊥ ≈ ν‖ ≈ νi/i. (3.23) But more significantly, the single-ion F-region ACF’s at 50 MHz are essentially the same as in collisionless and non-magnetized plasmas because (a) the ions move by many Bragg wavelengths λB = 2pi/k between successive Spitzer collisions, and (b) the ions are unable to return to within λB/2pi of their starting positions after a gyroperiod as a consequence of the ion-ion interactions (giving rise to Spitzer collisions). As a result, ion diffusion is sufficient to eliminate the gyroresonance peaks from the single-ion ACF’s, as we have shown above. As an upshot, we will be able to handle the ion terms analytically in the general framework equations. 3.3.2 Statistics of the Electron Displacements Next, we study the effects of Coulomb collisions on electron trajectories using pro- cedures similar to those applied to ions. In Figure 3.6, the displacement distributions resulting from the simulation of an electron moving in an O+ plasma are presented. The plasma parameters considered in this simulation are also given in Table 3.1. The top and bottom panels in Figure 3.6 correspond, respectively, to the distributions in the directions perpendicular and parallel to B. On the left, the distributions are displayed as functions of τ , while on the the right, sample cuts of the distributions are compared to a Gaus- sian pdf. As in the ion case, we note that the normalized distributions for the direction perpendicular to B do not vary with τ and match almost perfectly to Gaussian curves. However, for displacements parallel to B, the normalized distributions do vary with τ , and the shapes are distinctly non-Gaussian for intermediate values of τ . More specifically, at very small time delays (τ . 1µs), the distributions are Gaussian, but then, in less than a millisecond, the distribution curves become more “spiky” (positive kurtosis) than a Gaus- sian. Although, at even longer delays τ the distributions once again relax to a Gaussian shape, it is clear that the electron displacement in the direction parallel to B is not a Gaussian random variable at all time delays τ . Will a Brownian-motion model still prove useful to fit the simulation data for the elec- 37 0 2 4 6 8 10−2 0 2 0 0.2 0.4 Time delay [ms] Electron displacement distributions ⊥ to B ∆r⊥(τ )/σ⊥(τ ) ∆ r ⊥ p d f 0 2 4 6 8 10−2 0 2 0 0.2 0.4 Time delay [ms] Electron displacement distributions ‖ to B ∆r‖(τ )/σ‖(τ ) ∆ r ‖ p d f −3 −2 −1 0 1 2 3 0 0.1 0.2 0.3 0.4 0.5 ∆r⊥(τ )/σ⊥(τ ) ∆ r ⊥ p d f Electron displacement distributions ⊥ to B τ = 0.0ms τ = 2.0ms τ = 4.0ms τ = 6.0ms τ = 8.0ms τ = 10.0ms Gaussian −3 −2 −1 0 1 2 3 0 0.1 0.2 0.3 0.4 0.5 ∆r‖(τ )/σ‖(τ ) ∆ r ‖ p d f Electron displacement distributions ‖ to B τ = 0.0ms τ = 2.0ms τ = 4.0ms τ = 6.0ms τ = 8.0ms τ = 10.0ms Gaussian Figure 3.6 Same as Figure 3.3 but for the case of a simulated electron. All distributions at all time delays are normalized to unit variance. Note that the distributions of the displacements parallel to B become narrower than a Gaussian distribution in less than a millisecond. trons even though the displacement process seems non-Gaussian? To answer this question we first attempted to fit (3.18) and (3.19) to the simulated variances of the electron dis- placements (as we did earlier for the ion displacements). An example of our fitting results is presented in Figure 3.7. In general, we found that the simulated variance data are well fitted by the Brownian-motion expressions with the best-fit results of ν‖ ≈ νe/i (3.24) and ν⊥ ≈ νe/i + νe/e, (3.25) where νe/e = Ne e 4 ln Λe 12pi3/2 2om 2 e C 3 e (3.26) 38 0 0.005 0.01 0.015 0.020 1 2 3 x 10−3 〈∆ r 2 ⊥ 〉[ m 2 ] Variance of electron displacement ⊥ to B 〈|∆x|2〉 〈|∆y|2〉 Best fit ν⊥ = 2440.90 Hz 0 0.005 0.01 0.015 0.020 2 4 6 Time delay [ms] 〈∆ r 2 ‖〉 [m 2 ] Variance of electron displacement ‖ to B 〈|∆z |2〉 Best fit ν‖ = 1468.74 Hz 0 2 4 6 8 100 0.02 0.04 〈∆ r 2 ⊥ 〉[ m 2 ] Variance of electron displacement ⊥ to B 〈|∆x|2〉 〈|∆y|2〉 Best fit ν⊥ = 2440.90 Hz 0 2 4 6 8 100 1 2 x 105 Time delay [ms] 〈∆ r 2 ‖〉 [m 2 ] Variance of electron displacement ‖ to B 〈|∆z |2〉 Best fit ν‖ = 1468.74 Hz Figure 3.7 Same as Figure 3.4 but for the case of a simulated electron. In this case, the results are displayed for time intervals of 0.02 ms (left panels) and 10 ms (right panels). Note the different scales for the displacement variances. and νe/i = √ 2νe/e = √ 2Ne e 4 ln Λe 12pi3/2 2om 2 e C 3 e (3.27) are the Spitzer electron-electron and electron-ion collision frequencies [Callen, 2006] with Λe ≡ 12piNeh2ehD and electron Debye length he defined in Table 3.2. For instance, in Figure 3.7 the best-fit frequencies were ν‖ = 1.469 kHz and ν⊥ = 2.441 kHz, while νe/i = 1.439 kHz and νe/i + νe/e = 2.457 kHz. Next we examined whether the Brownian single-particle ACF model (3.17) can be used to fit the simulated electron ACF’s using the best-fit friction coefficients identified above. In Figure 3.8 we present the single-electron ACF’s that were computed for different magnetic aspect angles and λB = 3 m. The blue curves correspond to the ACF’s obtained from our simulations, while the green curves are the electron ACF’s calculated using expression (3.17) together with our approximations for ν‖ and ν⊥. Additionally, the electron ACF’s for a collisionless magnetized plasma are also plotted (red curves). We note that the simulated and the Brownian ACF’s matched almost perfectly at α = 0◦, and also that the agreement is still good at very small magnetic aspect angles (see panels a, b, and c). However, substantial differences between the Brownian and simulated ACF’s become evident as the magnetic aspect angle increases (see panels d, e, and f). In summary, our study revealed that the single-electron ACF’s needed for ISR spectral 39 0 50 10 0 15 0 20 0 25 0 30 0 0 0. 2 0. 4 0. 6 0. 81 T im e [m s] (a ) E le c tr o n A C F (α = 0 ◦ ) S im u la ti o n B ro w n ia n 0 50 10 0 15 0 20 0 25 0 30 0 0 0. 2 0. 4 0. 6 0. 81 T im e [m s] (b ) E le c tr o n A C F (α = 0 .0 1 ◦ ) S im u la ti o n B ro w n ia n N o -c o ll is io n s 0 10 20 30 40 50 60 0 0. 2 0. 4 0. 6 0. 81 T im e [m s] (c ) E le c tr o n A C F (α = 0 .0 5 ◦ ) S im u la ti o n B ro w n ia n N o -c o ll is io n s 0 5 10 15 20 25 0 0. 2 0. 4 0. 6 0. 81 T im e [m s] (d ) E le c tr o n A C F (α = 0 .1 ◦ ) S im u la ti o n B ro w n ia n N o -c o ll is io n s 0 0. 5 1 1. 5 2 2. 5 0 0. 2 0. 4 0. 6 0. 81 T im e [m s] (e ) E le c tr o n A C F (α = 0 .5 ◦ ) S im u la ti o n B ro w n ia n N o -c o ll is io n s 0 0. 2 0. 4 0. 6 0. 8 1 0 0. 2 0. 4 0. 6 0. 81 T im e [m s] (f ) E le c tr o n A C F (α = 1 ◦ ) S im u la ti o n B ro w n ia n N o -c o ll is io n s F ig ur e 3. 8 Si m ul at ed el ec tr on A C F ’s fo r λ B = 3 m at di ffe re nt m ag ne ti c as pe ct an gl es : (a ) α = 0◦ , (b ) α = 0. 01 ◦ , (c ) α = 0. 05 ◦ , (d ) α = 0. 1 ◦ , (e ) α = 0. 5 ◦ , an d (f ) α = 1◦ . N ot e th e di ffe re nt ti m e sc al es us ed in ea ch pl ot . 40 models cannot be accurately modeled as a Brownian-motion process over all aspect angles. The fundamental reason for this is the deviation of the electron displacements parallel to B from Gaussian statistics, despite the fact that displacement variances are well fitted by the Brownian model. Certainly, a non-Gaussian process cannot be fully characterized by a model that specifies its first and second moments only; this is particularly true for the estimation of the characteristic function of the process 〈ejk·∆re〉 that depends on all the moments of the process distribution. What is the physical reason for the electron displacements to be non-Gaussian in the direction parallel to B at intermediate time delays? We believe the answer is related to the fact that at low and high electron speeds (in comparison with Ce), the Fokker-Planck collision coefficients are dominated by different types of physical processes corresponding to electron-ion and electron-electron interactions, respectively. While the electron-electron collisions dominant at high speeds give rise to random changes in the electron speed, the same quantity is conserved in electron-ion collisions (due to the large mass difference of electrons and ions) dominating at low speeds. This asymmetry leads to the forma- tion of a low-speed population of electrons (roughly after one electron-electron collision time) staying close to their starting locations for longer periods of time (during which many electron-ion collisions may take place) than would have been expected in a (speed independent) Brownian-motion model for Coulomb collisions. As a result, the electron dis- placements at intermediate delays τ have distributions that are sharper than a Gaussian, a shape which is consistent with an enhanced population at small electron displacements. This population eventually disperses, and the distributions relax back to a Gaussian form at large time delays corresponding to many electron-electron collision times, as would be expected in view of the central limit theorem. This description of what we suspect is going on is also consistent with ν‖ fitting best νe/i (as opposed to νe/i + νe/e) at the interme- diate time scales of importance to our fitting procedure. Note that despite our best-fit result ν‖ ≈ νe/i, electron-electron Coulomb collisions still play a crucial role in shaping the single-electron ACF by determining the time scale over which the displacements parallel to B return back to normal (Gaussian) distribution. 41 As for the absence of similar effects in the direction perpendicular to B, this can be explained by the fact that the main role the collisions play in the dynamics of the electrons in the perpendicular plane is to knock them off their gyrocenters, a process that is not sensitive to the distinctions between electron-electron and electron-ion collisions (and hence ν⊥ ≈ νe/i + νe/e as we have found out). Returning back to Figure 3.8, the fact that the correlation time of the simulated electron signal is longer than what is predicted by the other models is a signature effect of Coulomb collisions encountered at small (but non-zero) aspect angles. As discussed in the next chapter, this feature of the electron ACF determines the width of the electron Gordeyev integral; therefore, it has an impact on the shape of the incoherent scatter spectrum averaged over small magnetic aspect angles. Since the Brownian-motion model cannot be used for the calculation of electron ACF’s and no other simplified model was found, a numerical library of electron ACF’s had to be built, as already described in the previous section. The collisional IS spectra to be presented in the next chapter were produced using the numerical library. 42 CHAPTER 4 COLLISIONAL INCOHERENT SCATTER SPECTRUM: RESULTS AND COMPARISONS In the present chapter, ion and electron Gordeyev integrals are utilized in the gen- eral framework formulation of Chapter 2 to compute incoherent scatter spectra that take into account Coulomb collision effects. The Gordeyev integrals are computed from their corresponding single-particle ACF’s using the chirp-Z transform algorithm detailed in the previous chapter. The electron Gordeyev integral, which resembles the shape of the IS spectrum in the direction perpendicular to B, is analyzed first in Section 4.1. The char- acteristics of the collisional spectra are discussed in Section 4.2, where our results are also compared to other IS spectrum models. 4.1 Simulated Electron Gordeyev Integrals Figure 4.1 shows the plots of Re{Je(ω)}, the real part of the electron Gordeyev integral proportional to 〈|nte(k, ω)|2〉, computed using our new numerical library as a function of Doppler frequency and magnetic aspect angle for two different Bragg wavelengths, λB = 3 m and λB = 0.3 m. Notice the difference between the frequency axes used in these plots. In the first one, we are considering a range of ω/2pi from −1 kHz to 1 kHz, while in the second plot, the frequency range is ten times wider. A few degrees away from perpendicular to B, where the effect of collisions can be neglected, the bandwidth of the electron Gordeyev integral Je(ω) is proportional to the product of Bragg wavenumber k and electron thermal speed Ce. Thus, if the radar wavelength is reduced by a factor of ten, the bandwidth of Je(ω) will increase by the same factor. If the same were true for viewing directions close 43 Frequency [kHz] A sp ec t a n g le [d eg ] Electron Gordeyev integral (λB = 3m) −1 −0.5 0 0.5 10 0.1 0.2 0.3 0.4 0.5 R e {J e }[ d B ] −45 −40 −35 −30 −25 −20 −15 −10 −5 0 Frequency [kHz] A sp ec t a n g le [d eg ] Electron Gordeyev integral (λB = 0.3m) −10 −5 0 5 100 0.1 0.2 0.3 0.4 0.5 R e {J e }[ d B ] −45 −40 −35 −30 −25 −20 −15 −10 −5 0 (a) (b) Figure 4.1 Simulated electron Gordeyev integrals as functions of Doppler frequency and magnetic aspect angle for two radar Bragg wavelengths: (a) λB = 3 m and (b) λB = 0.3 m. An O+ plasma with electron density Ne = 1012 m−3, temperatures Te = Ti = 1000 K, and magnetic field Bo = 25µT is considered. to perpendicular to B, the plots displayed above would have looked identical. However, that is not the case, and our results illustrate that the dependence of Je(ω) on k and Ce deviates from kCe at small aspect angles. We observe that both Gordeyev integrals become narrower as the magnetic aspect angle decreases. However, in the case of λB = 0.3 m, Re{Je(ω)} stops shrinking at an angle around 0.02◦, while in the case of λB = 3 m, the same happens at a smaller angle (∼ 0.01◦). This is another effect of Coulomb collisions, namely the saturation of electron Gordeyev integral width at very small α. Without having considered collisions, Re{Je(ω)} for α = 0◦ would have approached a Dirac delta function at ω = 0 (accompanied by a series of smaller deltas located at multiples of the electron gyrofrequency Ωe). Trying out the Brownian-motion expression for the electron ACF as a guide, we noticed that the bandwidth of Re{Je(ω)} in the limit α = 0◦ varies according to k2C2e ν⊥ ν2⊥ + Ω2e . (4.1) This expression is obtained by placing (3.19) into (3.17) and considering ν⊥τ  1, a limit in which the electron ACF becomes an exponential function. Notice that the bandwidth 44 Figure 4.2 Simulated incoherent scatter spectra as a function of magnetic aspect angle and Doppler frequency for λB = 3 m (e.g., Jicamarca radar Bragg wavelength). An O+ plasma with physical parameters given in Table 3.1 is considered. dependence on kCe has changed from linear (at relatively large α) to quadratic (at very small α), a change that is related to the transition between the collisionless and collisional regimes of the electron diffusion. Furthermore, using ν⊥ ≈ νe/i+ νe/e from the last section and taking Ωe  ν⊥, we can verify that bandwidth dependence (4.1) is proportional to Ne√ Te . (4.2) Since at very small aspect angles the electron Gordeyev integral dominates the shape of the overall incoherent scatter spectrum, we expect the IS spectral width to exhibit the same dependence (4.2) on density and temperature. We will demonstrate this to be the case by specific examples shown in the following section. 4.2 Simulated Collisional IS Spectra Figure 4.2 shows a surface plot constructed from full IS spectrum calculations for λB = 3 m (e.g., Jicamarca radar) as a function of aspect angle α and Doppler frequency 45 ω/2pi. In this self-normalized surface plot, we observe how the IS spectrum sharpens significantly at small aspect angles. For instance, just in the range between 0.1◦ and 0◦, the amplitude of the spectrum becomes ten times larger while its bandwidth is reduced by the same factor. Sample cuts of the surface plot taken at a number of magnetic aspect angles are dis- played in Figure 4.3. In these plots, the simulated spectra (blue curves) are compared to other spectral models. The green curves correspond to spectra computed using the Brownian-motion model discussed in the previous section, while the red curves correspond to the collisionless electron case (i.e., the case when Coulomb collisions for the electrons are neglected but collisions for the ions are still considered). We note that the agreement with the Brownian-motion spectra is excellent at very small aspect angles (α ≤ 0.01◦). However, as α increases, our collisional spectra become narrower than the Brownian-motion spectra. We also observe that the collisionless spectra are in general wider than our simulation results. Additionally, note that at very small aspect angles, the collisionless spectrum has a humped shape (similar to the one expected at large aspect angles but with a much nar- rower spectral width). In this regime, electrons and ions exchange their roles in defining the shape of the collisionless spectrum so that the spectral width becomes proportional to the electron thermal speed Ce and sinα, as previously discussed by Kudeki et al. [1999]. In the limit of α → 0◦, the collisionless spectrum develops a delta function ω = 0 implying an infinite signal correlation time in the collisionless limit. The reason for this is that in the absence of collisions, the magnetic field restricts the motion of the electrons in the plane perpendicular to B, forcing them to gyrate always around the same magnetic field lines. With collisions, the electrons manage to diffuse across the field lines, and consequently the correlation time of the IS signal becomes finite and its spectrum broadens out of its delta function limit. However, collisions do not cause spectral broadening at all aspect angles. Notice that in Figure 4.3 the collisional spectrum for α > 0.01◦ is narrower than the collisionless spectrum, implying that collisions in this range of aspect angles have caused an increase of the IS signal correlation time. The explanation for this is that at aspect angles larger than 0.01◦ the shape of the IS spectrum is dominated by electron 46 0 5 10 15 20 0246 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e (a ) IS -S P C (α = 0 ◦ ) × 10 − 2 S im u la ti o n B ro w n ia n W o o d m a n [2 0 0 4 ] 0 5 10 15 20 0123456 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e (b ) IS -S P C (α = 0 .0 1 ◦ ) × 10 − 2 S im u la ti o n B ro w n ia n N o -c o ll is io n s W o o d m a n [2 0 0 4 ] 0 25 50 75 10 0 0 0. 51 1. 5 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e (c ) IS -S P C (α = 0 .0 5 ◦ ) × 10 − 2 S im u la ti o n B ro w n ia n N o -c o ll is io n s W o o d m a n [2 0 0 4 ] 0 50 10 0 15 0 20 0 25 0 0123456 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e (d ) IS -S P C (α = 0 .1 ◦ ) × 10 − 3 S im u la ti o n B ro w n ia n N o -c o ll is io n s W o o d m a n [2 0 0 4 ] 0 25 0 50 0 75 0 10 00 12 50 0 0. 2 0. 4 0. 6 0. 81 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e (e ) IS -S P C (α = 0 .5 ◦ ) × 10 − 3 S im u la ti o n B ro w n ia n N o -c o ll is io n s S & G [1 9 9 9 ] 0 30 0 60 0 90 0 12 00 15 00 0123456 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e (f ) IS -S P C (α = 1 ◦ ) × 10 − 4 S im u la ti o n B ro w n ia n N o -c o ll is io n s S & G [1 9 9 9 ] F ig ur e 4. 3 Sa m pl e cu ts of th e si m ul at ed sp ec tr a of F ig ur e 4. 2 ta ke n at di ffe re nt m ag ne ti c as pe ct an gl es : (a ) α = 0◦ , (b ) α = 0. 01 ◦ , (c ) α = 0. 05 ◦ , (d ) α = 0. 1 ◦ , (e ) α = 0. 5◦ , an d (f ) α = 1◦ . O ur si m ul at io n re su lt s (b lu e lin es ) ar e co m pa re d to IS sp ec tr a co m pu te d us in g th e B ro w ni an -m ot io n m od el (g re en lin es ), th e co lli si on le ss el ec tr on m od el (r ed lin es ), th e si m ul at io n re su lt s of Su lz er an d G on zá le z [1 99 9] (l ig ht -b lu e lin es ), an d th e as pe ct an gl e de pe nd en t co lli si on fr eq ue nc y m od el of W oo dm an [2 00 4] (m ag en ta lin es ). N ot e th e di ffe re nt fr eq ue nc y sc al es us ed in ea ch pl ot . 47 diffusion in the direction parallel to B. As collisions impede the free motion of particles, electrons diffuse slower in a collisional plasma than in a collisionless one, which implies that the electrons stay closer to their original locations for longer periods of time. As a result, the correlation time of the signal scattered by the electrons also becomes longer, causing the broadening of the IS signal ACF and the associated narrowing of the signal spectrum in this aspect angle regime, as first explained by Sulzer and González [1999]. Also in Figure 4.3, we compare our simulated spectra to the results of Sulzer and González [1999] for α ≥ 0.5◦ (see panels e and f). The comparison is good except for a minor offset which was traced to a minor coding error in the simulation program used by Sulzer and González [1999]. At smaller aspect angles, comparisons are made with the spectrum model of Woodman [2004], which is effectively an ad hoc extrapolation of the Sulzer and González [1999] model to very small magnetic aspect angles (including the α→ 0◦ limit). Woodman [2004] fitted the Brownian-motion spectral model to the results of Sulzer and González [1999] and developed an empirical formula for the electron collision frequency that depends on the magnetic aspect angle. As can be observed in panels a, b, c, and d, our simulated spectra are not quite as sharp as Woodman’s spectra (though the differences are small). These differences can be expected since Woodman [2004] did not really have the spectra to fit at very small aspect angles; therefore, his collision frequency formula may require a little adjustment for α < 0.1◦. Figure 4.4 shows the dependence of the simulated IS spectra on the mean electron density Ne that controls the Coulomb collision rates. As discussed in the previous section, we expect the spectrum at α = 0◦ to become wider as Ne grows (see relation (4.2)). This behavior is clearly illustrated in the left panel of Figure 4.4. As Ne increases, collisions become more frequent and consequently the electrons diffuse out to longer distances across the magnetic field lines. Hence, the spectrum bandwidth increases with Ne, so long as the collision frequency remains smaller than the gyrofrequency, otherwise the collisions start impeding electron motion rather than just perturbing the gyromotions to facilitate diffusion. In a highly collisional regime, characterized by collisional forces much larger than magnetic forces (or by ν⊥  Ωe), electrons will find it very difficult to move in any 48 0 1 2 3 4 50 1 2 3 4 5 Frequency (ω/2pi) 〈|n e (k ,ω )|2 〉/ N e IS-SPC (α = 0◦) × 10−1 log10Ne = 11.0 log10Ne = 11.5 log10Ne = 12.0 log10Ne = 12.5 0 50 100 150 2000 2 4 6 Frequency (ω/2pi) 〈|n e (k ,ω )|2 〉/ N e IS-SPC (α = 0.1◦) × 10−3 No-collisions log10Ne = 11.0 log10Ne = 11.5 log10Ne = 12.0 log10Ne = 12.5 0 500 1000 15000 2 4 6 Frequency (ω/2pi) 〈|n e (k ,ω )|2 〉/ N e IS-SPC (α = 1◦) × 10−4 No-collisions log10Ne = 11.0 log10Ne = 11.5 log10Ne = 12.0 log10Ne = 12.5 Figure 4.4 Electron density dependence of the simulated IS spectra for λB = 3 m at aspect angles α = 0◦ (left panel), α = 0.1◦ (central panel), and α = 1◦ (right panel). An O+ plasma in thermal equilibrium with Te = Ti = 1000 K and Bo = 25µT is considered. Note the different frequency scales used in each plot. direction. As a result, the bandwidth of the IS spectrum at α = 0◦ will stop increasing with Ne and eventually it will decrease. A similar effect can be observed as the aspect angle increases. A few hundredths of a degree away from perpendicular to B, the shape of the IS spectrum starts being controlled by diffusion along the magnetic field lines. In this regime, the spectrum becomes narrower as Ne grows because the increasing number of collisions impede the movement of the electrons (as already discussed above). This behavior is observed in the central panel of Figure 4.4. Also note that the collisional spectrum approaches the collisionless spectrum as Ne is reduced, which is consistent with the dependence of collision frequencies on Ne. Finally, a comparison of the central and rightmost panels of Figure 4.4 shows that collision effects become less significant at even larger aspect angles where the spectrum is shaped by ion dynamics. In that regime, the spectral shapes become independent of Ne as long as khe  1. The temperature dependence of the simulated IS spectra at different aspect angles will be examined with the help of Figure 4.5. Considering an O+ plasma with Ne = 1012 m−3 and Ti = 1000 K, we have computed the IS spectra for a Bragg wavelength of λB = 3 m and for a set of electron temperatures Te (1000 K, 1500 K, 2000 K, 2500 K, and 3000 K). Figure 4.5 displays the results for three different aspect angles, α = 0◦ (left panels), α = 0.02◦ (center panels), and α = 0.1◦ (right panels). Notice that the same spectra are plotted in the top and bottom panels, but in the bottom panels we have normalized each spectrum by its peak value in order to emphasize the changes in spectral widths. 49 0 2 4 6 8 10 0 0. 51 1. 52 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e IS -S P C (α = 0 ◦ ) × 10 − 1 T e / T i = 1 .0 T e / T i = 1 .5 T e / T i = 2 .0 T e / T i = 2 .5 T e / T i = 3 .0 0 10 20 30 40 50 0123 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e IS -S P C (α = 0 .0 2 ◦ ) × 10 − 2 T e / T i = 1 .0 T e / T i = 1 .5 T e / T i = 2 .0 T e / T i = 2 .5 T e / T i = 3 .0 0 50 10 0 15 0 20 0 012345 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/N e IS -S P C (α = 0 .1 ◦ ) × 10 − 3 T e / T i = 1 .0 T e / T i = 1 .5 T e / T i = 2 .0 T e / T i = 2 .5 T e / T i = 3 .0 0 2 4 6 8 10 0 0. 2 0. 4 0. 6 0. 81 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/〈|n e (k,0)| 2 〉 IS -S P C (α = 0 ◦ ) T e / T i = 1 .0 T e / T i = 1 .5 T e / T i = 2 .0 T e / T i = 2 .5 T e / T i = 3 .0 0 10 20 30 40 50 0 0. 2 0. 4 0. 6 0. 81 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/〈|n e (k,0)| 2 〉 IS -S P C (α = 0 .0 2 ◦ ) T e / T i = 1 .0 T e / T i = 1 .5 T e / T i = 2 .0 T e / T i = 2 .5 T e / T i = 3 .0 0 50 10 0 15 0 20 0 0 0. 2 0. 4 0. 6 0. 81 F re q u e n c y (ω / 2 pi ) 〈|n e (k,ω)| 2 〉/〈|n e (k,0)| 2 〉 IS -S P C (α = 0 .1 ◦ ) T e / T i = 1 .0 T e / T i = 1 .5 T e / T i = 2 .0 T e / T i = 2 .5 T e / T i = 3 .0 F ig ur e 4. 5 T e / T i de pe nd en ce of th e si m ul at ed IS sp ec tr a fo r λ = 3 m at as pe ct an gl es α = 0◦ (l ef t pa ne ls ), α = 0. 02 ◦ (c en tr al pa ne ls ), an d α = 0 .1 ◦ (r ig ht pa ne ls ). T he sa m e sp ec tr a ar e pl ot te d on th e to p an d bo tt om pa ne ls ; ho w ev er , on th e bo tt om , ev er y sp ec tr um ha s be en no rm al iz ed to it s pe ak va lu e. A n O + pl as m a w it h N e = 10 1 2 m −3 , T i = 10 00 K , an d B o = 25 µ T is co ns id er ed . E ac h cu rv e co rr es po nd s to di ffe re nt va lu es of el ec tr on te m pe ra tu re s T e (1 00 0 K ,1 50 0 K ,2 00 0 K ,2 50 0 K ,a nd 30 00 K ). N ot e th e di ffe re nt fr eq ue nc y sc al es us ed in ea ch co lu m n. 50 Figure 4.5 includes a wealth of information pertinent to the ultimate application of our theory in estimating the ionospheric temperatures from Jicamarca radar data collected with beams perpendicular to B. First, we note that the left panels of Figure 4.5 show how the spectrum bandwidth at α = 0◦ is inversely proportional to the square root of Te as discussed earlier. Second, we note that this dependence of spectral width on Te changes at larger aspect angles as shown in the remaining panels of the figure — for aspect angles larger than about α = 0.01◦, spectral width increases with increasing Te. Third, the spectral amplitudes change rapidly with aspect angle. For instance, in the case of Te = Ti = 1000 K, the peak of the spectrum at α = 0◦ is about ten times larger than the one at α = 0.1◦. Finally, the amplitude differences for different α are enhanced at higher values of Te (in comparison with Ti), which is related to the aspect angle dependence of the volumetric radar cross section (RCS) [e.g., Farley , 1966; Milla and Kudeki , 2006] σv ≡ 4pir2eNeη(k) (4.3) of incoherent scatter radar echoes. In (4.3), re is the classical electron radius and η(k) is an electron scattering efficiency factor [see Milla and Kudeki , 2006] defined as η(k) ≡ ˆ dω 2pi 〈|ne(k, ω)|2〉 Ne , (4.4) which is primarily dependent on the temperature ratio Te/Ti and magnetic aspect angle α. A plot of this factor obtained from our collisional simulation results is shown in Figure 4.6. As we can observe, if the plasma is in thermal equilibrium (i.e., if Te = Ti), this factor is 1/2 at all angles α. We can also appreciate that η(k) at α = 0◦ is clearly increasing as Te/Ti grows. On the other hand, at very large magnetic aspect angles, we can see that the efficiency factor decreases as Te/Ti increases. In particular, note that our calculations for α = 90◦ match the well-known formula (1 + Te/Ti)−1, as expected for moderate values of Te/Ti and negligible Debye length [e.g., Farley , 1966]. Note that for α ≈ 1◦ the factor is approximately independent of Te/Ti, but otherwise it increases and decreases with the 51 1 1.5 2 2.5 30.2 0.3 0.4 0.5 0.6 0.7 0.8 Te/Ti η (k ) α = 0. 00 ◦ 0.10 ◦ 0.25 ◦ 0.50 ◦ 1.00◦ 2.00◦ 5.00◦ 90.00◦ Figure 4.6 Electron scattering efficiency factor η(k) resulting from the frequency integra- tion of the simulated incoherent scatter spectra as a function of electron-to-ion temperature ratio Te/Ti and magnetic aspect angle α. An O+ plasma withNe = 1012 m−3, Ti = 1000 K, and Bo = 25µT is considered. temperature ratio at small and large aspect angles, respectively. The efficiency factor η(k) plays an important role in Ne estimation from ISR power data, and its accurate determi- nation in the collisional regime at small aspect angles is one of the main contributions of this work. The next stage of our studies, which is described in the following chapters, is the modeling of incoherent scatter spectrum measurements carried out with Jicamarca antenna beams pointed perpendicular to B. For this purpose, the measured spectrum has to be carefully modeled because within the range of small magnetic aspect angles illuminated by the antenna beams, the theoretical IS spectra vary quite rapidly. In addition, magneto- ionic propagation effects are also considered into the model as these effects are significant at 50 MHz. This further complicates the description of Jicamarca measurements as different modes of propagation are excited within the width of the antenna beams; these modes vary from linearly polarized (at α = 0◦) to circularly polarized (at α > 1◦) with the transition between these two regimes happening at an angle of approximately 0.5◦. The modeling of the beam-weighted incoherent scatter radar spectrum including magneto-ionic propagation and collisional effects is the subject of the next chapter of this dissertation. 52 CHAPTER 5 MAGNETO-IONIC PROPAGATION EFFECTS ON THE INCOHERENT SCATTER RADAR MEASUREMENTS A radiowave propagating through the ionosphere experiences changes in its polarization caused by the presence of the Earth’s magnetic field. In this chapter, a model for incoherent scatter spectrum and cross-spectrum measurements that takes into account magneto-ionic propagation effects is developed. A mathematical description of radiowave propagation in a homogeneous magneto- plasma is derived in the first section of this chapter. This basic formulation is extended to the case of an inhomogeneous ionosphere in Section 5.2. The resultant wave propagation model is used in Section 5.3 to reformulate the soft-target radar equation of Chapter 2 in order to account for the magneto-ionic propagation effects on incoherent scatter spectrum and cross-spectrum measurements. 5.1 Propagation of Electromagnetic Waves in a Homogeneous Magnetoplasma In this section, the characteristic modes of electromagnetic wave propagation in a ho- mogeneous magnetoplasma are derived. We start by formulating Maxwell’s curl equations 53 in differential form ∇×E = −∂B ∂t , (5.1) ∇×H = J + ∂D ∂t . (5.2) Considering wave field solutions proportional to ejωt−jk·r, where ω is the wave frequency and k is the wave propagation vector, we can reformulate the previous equations in plane wave form as k×E = ωµoH, (5.3) −jk×H = J + jωoE, (5.4) where the differential operators ∇ and ∂/∂t have been replaced by −jk and jω, respec- tively. Above, we have considered B = µoH and D = oE, which is a valid assumption as long as current J includes components that are due to free and/or bound charged carri- ers represented in terms of conductivity and/or susceptibility tensors. Eliminating H by taking k× k×E, we have that k× k×E = jωµo [J + jωoE] . (5.5) Above, the propagation vector can be written as k = kon kˆ, where ko ≡ ω√µoo is the wavenumber in free-space, n is the refractive index of the medium, and kˆ is the unit vector parallel to the propagation direction. Now, we can rewrite (5.5) in the following form n2 kˆ× kˆ×E︸ ︷︷ ︸ −E⊥ = − 1 jωo J−E, (5.6) which simplifies to E− n2E⊥ = − 1 jωo J, (5.7) 54 where E⊥ is the component of E that is perpendicular to kˆ. Without loss of generality, we consider kˆ = zˆ; thus, breaking the previous equation into its components, we have that (1− n2)Ex = − Jx jωo , (5.8) (1− n2)Ey = − Jy jωo , (5.9) Ez = − Jz jωo . (5.10) Based on these equations, we can define the following polarization relations R = Ex Ey = Jx Jy (5.11) S = Ez Ex = Jz Jx (1− n2). (5.12) In addition, the refractive index can be written as n2 = 1 + Jx/Ex jωo = 1 + Jy/Ey jωo . (5.13) In order to obtain the refractive index of the medium, we need to define the relations Jx/Ex or Jy/Ey from an appropriate conductivity model. For instance, in an isotropic conductor, the current and electric field vectors are related by J ≡ σE, where σ is the conductivity. Thus, Jx/Ex = Jy/Ey = σ, which implies that n2 = 1 + σjωo , as is well known. Note that J = σE also implies that there is no coupling between Ex and Ey. Therefore, R could take any value, but S is equal to zero because there is no current along the propagation direction. In order to determine the conductivity of a magnetized plasma with a DC magnetic field Bo, let us consider an electron moving in this medium under the influence of the Lorentz force. Then, the equation of motion is given by me dv dt = −e(E + v ×Bo), (5.14) where me and −e are the mass and charge of the electron. In phasor form, the motion 55 Bo = By yˆ + Bz zˆ k θ Ek Eφ Eθ x y z θˆ = −yˆ φˆ = xˆ Figure 5.1 Coordinate system used for analyzing wave propagation in a magnetized plasma. The magnetic field Bo is on the yz-plane at an angle θ from the propagation vector k which is parallel to the zˆ-direction. The wave field E has three mutually orthog- onal components Ek, Eθ, and Eφ in directions kˆ = zˆ, θˆ = −yˆ, and φˆ = xˆ, respectively. equation becomes jωmev = −e(E + v ×Bo), (5.15) which can be written as − e me E = jωv + e me v ×Bo. (5.16) Next, multiplying both sides by −eNe and defining J ≡ −eNev, we have that Nee 2 jωme E = J− e jωme Bo × J. (5.17) In the coordinate system of Figure 5.1, we have considered Bo = Byyˆ +Bzzˆ with compo- nents By = Bo sin θ and Bz = Bo cos θ, where Bo is the magnitude of the field and θ is the angle with respect to the z-axis. Now, we can recast the previous equation as − jωoXE = J + jYT yˆ × J + jYLzˆ× J (5.18) in which X ≡ ω 2 p ω2 , YL ≡ Ωe ω cos θ, and YT ≡ Ωe ω sin θ, (5.19) where ωp ≡ √ Neq2e/ome is the plasma frequency and Ωe = eBo/me is the electron 56 gyrofrequency. Expanding equation (5.18) into its vector components, we have that −jωoXEx = Jx + jYTJz − jYLJy, (5.20) −jωoXEy = Jy + jYLJx, (5.21) −jωoXEz = Jz − jYTJx. (5.22) Early in this section, we found that Jz + jωoEz = 0. Then, using this result together with (5.22), we can find that Jz Jx = j YT 1−X . (5.23) Next, using the polarization relation R = Ex/Ey = Jx/Jy and the previous expression, we can rewrite equations (5.20) and (5.21) in the following forms −jωoXEx Jx = 1− Y 2 T 1−X − jYL 1 R , (5.24) −jωoXEy Jy = 1 + jYLR. (5.25) Noting that Ex/Jx = Ey/Jy, we combine the equations above and find that jYLR+ Y 2T 1−X + jYL 1 R = 0, (5.26) an equation that can be written as jYL(1−X)R2 + Y 2TR+ jYL(1−X) = 0. (5.27) Defining F ≡ −jYLR, the previous equation becomes (1−X)F 2 − Y 2T F − (1−X)Y 2L = 0. (5.28) Solutions of (5.28) are given by F± = Y 2T ± √ Y 4T + 4Y 2 L (1−X)2 2 (1−X) , (5.29) 57 which are valid for any given ω and direction θ. In order to find the refractive indices, note that Ey Jy = −1 + jYLR jωoX = − 1− F jωoX . (5.30) Thus, inserting this relation into (5.13), we obtain the following refractive indices n2± = 1− X 1− F± , (5.31) which characterize the two possible modes of propagation in a magnetoplasma. This result is known as the Appleton-Hartree equation [e.g., Budden, 1961]. Each of these refractive indices, n+ and n−, is accompanied by an associated polarization vector. The relations between the components of these vectors are given by the polarization relations R and S defined before. In order to formulate these relations, we considered a coordinate system in which kˆ = zˆ, θˆ = −yˆ, and φˆ = xˆ. In this system, the electric field components are given by Ek, Eθ, and Eφ. Thus, after some math, we can find the following relations Ek Eθ = −SR = F 2± − Y 2L − F±Y 2T YLYT (1− F±) = F±XYT (1− F±) (1−X)YL , (5.32) Ek Eφ = S = j F 2± − Y 2L − F±Y 2T YTF± (1− F±) = j XYT (1− F±) (1−X) , (5.33) Eθ Eφ = − 1 R = j YL F± , (5.34) where Y = √ Y 2L + Y 2 T = Ωe/ω. Let us discuss now some particular cases. In the limit of θ = 0◦ corresponding to propagation parallel to the ambient magnetic field Bo (i.e., longitudinal propagation), we have that YL = Y and YT = 0. Thus, F+ = Y and F− = −Y, which implies that refractive indices (5.31) simplify to n+ = √ 1− X 1− Y and n− = √ 1− X 1 + Y . (5.35) The polarization vectors associated with these modes are given by θˆ−jφˆ√ 2 and θˆ+jφˆ√ 2 , which correspond to right- and left-circularly polarized waves. Note that the polarization vectors 58 are transverse to the propagation direction. Additionally, in the limit of θ = 90◦ corre- sponding to propagation perpendicular to Bo (i.e., transverse propagation), we have that YL = 0 and YT = Y. This implies that F+ = Y 2 1−X and F− = 0, which leads to the following expressions for the refractive indices: n+ = √ 1− X 1− F+ and n− = √ 1−X. (5.36) While the polarization of the mode associated with n− is linear and parallel to the zˆ- direction (i.e., parallel to Bo), the polarization vector of n+ has components in kˆ- and φˆ-directions. If X  1 and Y  1, the kˆ-component can be considered negligible and the polarization becomes linear and parallel to φˆ (a direction that is perpendicular to Bo). Note that n− is equal to the refractive index of a plasma without magnetic field; because of this, the modes associated with n− and n+ are known as the ordinary and extraordinary modes (O- and X-modes), respectively. In general, for propagation at intermediate angles θ, the propagation modes are ellip- tically polarized. However, if X  1 and Y  1, e.g., the case of VHF propagation in the ionosphere, the propagation modes can be considered quasi-longitudinal for a wide range of angles θ. The approximation is valid as long as Y 4T  4Y 2L (1−X)2 ⇒ Y 2T  2YL (1−X) , (5.37) which implies the following condition: tan θ sin θ  21−X Y ≈ 2 Y . (5.38) Note that, since sin θ ≤ 1, a more stringent condition is tan θ  2/Y ; thus, we can define a critical angle θc ≡ tan−1 ( 2 Y ) (5.39) that corresponds to the transition between longitudinal and transverse modes of propa- gation. For instance in the case of Bo = 25µT and f = ω/2pi = 50 MHz, the electron 59 gyrofrequency is Ωe/2pi = 699.8 kHz and the critical angle is about 89.6◦. At Jicamarca, antenna beams pointed perpendicular to Bo have a beam-width of the order of a degree; thus, different modes of propagation ranging from longitudinal to transverse are excited within the range of aspect angles illuminated by the beam. As discussed above, in the presence of an ambient magnetic field, there are two possible modes of electromagnetic wave propagation. Labeling the modes as ordinary and extraor- dinary and denoting the refractive indices given in (5.31) as nO = n− and nX = n+, we can express the transverse component of a propagating electric field as Et = AO ( θˆ − jφˆFO YL ) e−jkonOr +AX ( θˆ − jφˆFX YL ) e−jkonXr, (5.40) where AO and AX are the amplitudes of the O- and X-waves. Above, r is the distance from the origin along the propagation direction kˆ and ko ≡ ω√µoo. In matrix notation, we can express (5.40) as  Eθ Eφ  =  e−jkonOr e−jkonXr −j FOYL e−jkonOr −j FX YL e−jkonXr   AO AX  , (5.41) where Eθ and Eφ are the transverse field components in θˆ and φˆ directions. Defining Eθ,o and Eφ,o as the field components at the origin, we have that  Eθ,o Eφ,o  =  1 1 −j FOYL −j FX YL   AO AX  . (5.42) Then, noting that FOFX Y 2L = −1, the propagating electric field (5.41) can be recast as  Eθ Eφ  =  e−jkonOr + a2e−jkonXr 1 + a2 ja ( e−jkonOr − e−jkonXr) 1 + a2 −ja (e−jkonOr − e−jkonXr) 1 + a2 a2e−jkonOr + e−jkonXr 1 + a2   Eθ,o Eφ,o  , (5.43) where the parameter a ≡ FOYL = − YL FX determines the type of wave polarization and can take values within the range 0 ≤ |a| ≤ 1. Note that the limits 0 and 1 correspond to the 60 cases of linearly and circularly polarized propagation modes. Defining n¯ ≡ nO+nX2 and ∆n ≡ nO−nX2 , we can rewrite (5.43) as Eθ Eφ  = e−jkon¯r 1 + a2  e−jko∆nr + a2ejko∆nr 2a sin(ko∆nr) −2a sin(ko∆nr) a2e−jko∆nr + ejko∆nr  ︸ ︷︷ ︸ T¯  Eθ,o Eφ,o  , (5.44) where T¯ is a propagator matrix that maps the fields at the origin into the fields at a distance r. Note that in the case of waves traveling in −kˆ direction, the same matrix T¯ can be used to propagate the fields from a distance r to the origin. The propagator matrix can be eigenvalue-decomposed as T¯ = P¯ S¯ P¯H where S¯ = e−jkon¯r  e−jko∆nr 0 0 ejko∆nr  (5.45) is a diagonal matrix composed of the eigenvalues of T¯, and P¯ = 1√ 1 + a2  1 −ja −ja 1  (5.46) is an unitary matrix with columns equal to the eigenvectors of T¯ (corresponding to the O- and X-modes, respectively). Returning to the xyz-coordinate system, we can express the transverse electric field as Et = Eθ θˆ + Eφ φˆ. Then, using (5.44) and the eigenvalue decomposition of T¯, we can show that Et = e−jkon¯r [ e−jko∆nrpOpHO + e jko∆nrpXp H X ] Eto, (5.47) where Eto is the wave field at the origin, and pO = θˆ − ja φˆ√ 1 + a2 and pX = −ja θˆ + φˆ√ 1 + a2 (5.48) are the polarization vectors of the O- and X-waves expressed in cartesian coordinates. 61 This form of the transverse field is suitable for computational evaluation because it only involves vector operations and there is no need to explicitly build the propagator matrix. 5.2 Model for Radiowave Propagation in an Inhomogeneous Ionosphere A radiowave propagating through the ionosphere will experience refraction and polar- ization effects. Wave refraction is caused by the altitudinal variation of the refractive index of the medium, which in this case is dependent on the amount of ionospheric electron den- sity. Additionally, wave polarization variations are caused by the presence of the Earth’s magnetic field, which is the controlling agent of the anisotropic nature of the ionospheric plasma. At VHF frequencies and for oblique propagation, refraction effects can be considered negligible for most propagation directions because the plasma frequency ωp at ionospheric altitudes is smaller than the wave frequency ω (i.e., X  1). It can be verified that this assumption is valid as long as the propagation angle ψ (measured with respect to the vertical axis) satisfies the condition tanψ  2 X . (5.49) Thus, refraction might have an impact only at very low elevation angles (an angular regime that is not of our interest). But for the same set of frequencies, polarization changes are still significant despite the fact that the electron gyrofrequency Ωe is much smaller than the wave frequency ω (i.e., Y  1). The reason for this is that the distances traveled by the propagating fields are long enough (hundreds of kilometers) so that phase differences between wave components propagating in distinct modes accumulate to detectable levels. Taking these elements into consideration and noting that, at VHF frequencies, the longi- tudinal components of the wave fields are negligibly small (as X  1 and Y  1), waves propagating through the ionosphere can be represented as TEM (transverse electromag- netic) waves. 62 zBi ∆r θi ri+1 ri ri−1 r3 r2 r1 k yx kˆ ψ Figure 5.2 Geometry of wave propagation in an inhomogeneous magnetized ionosphere. Consider plane wave propagation in an inhomogeneous magnetized ionosphere in an arbitrary direction kˆ. To model the electric field of the propagating wave, we can divide the ionosphere in slabs of equal width perpendicular to the propagation direction such that within each slab the physical plasma parameters (as electron density, electron and ion temperatures, and magnetic field) can be considered constants (see Figure 5.2).1 The transverse component of the wave electric field propagates from the bottom to the top of the i-th slab according to (5.47), Ei = e −jkon¯i∆r [ e−jko∆ni∆rpOpHO + e jko∆ni∆rpXp H X ] ︸ ︷︷ ︸ T¯i Ei−1, (5.50) which is the superposition of the O- and X-modes of magneto-ionic propagation detailed in the previous section. Above, T¯i denotes the i-th propagator matrix (expressed in cartesian coordinates), where ko ≡ 2pi/λo is the free-space wavenumber and ∆r is the width of the slab, and where n¯i ≡ nO,i+nX,i2 and ∆ni ≡ nO,i−nX,i 2 are the mean and half difference between the refractive indices of the propagation modes in the i-th layer. The polarization 1In the ionosphere, electron density and plasma temperatures can be considered to be functions of altitude f(z). Thus, the values of these physical parameters at any position r from a radar placed at the origin are given by f(r cosψ) where r is the radar range and ψ is the zenith angle. 63 vectors of the O- and X-modes are pO = θˆi − jai φˆi√ 1 + a2i and pX = −jai θˆi + φˆi√ 1 + a2i , (5.51) where ai ≡ FO,iYL,i = − YL,i FX,i is the polarization parameter, and θˆi and φˆi are a pair of mu- tually orthogonal unit vectors perpendicular to kˆ whose directions depend on the relative orientation of the propagation vector k and the magnetic field Bi (see Figure 5.1). At the interfaces between the slabs, the propagating fields are totally transmitted (i.e., wave reflection is neglected), an assumption that is valid because ionospheric plasma parameters change very slowly with altitude (having scale heights much larger than the wavelength of the propagating fields). Therefore, the field components of an upgoing plane wave propa- gating in the +kˆ direction (at a distance ri = i∆r from the origin) can be computed by the successive application of the propagator matrices; that is, Eui = T¯i · · · T¯2T¯1Euo , (5.52) where Euo is the wave field at the origin (perpendicular to kˆ), and T¯1 . . . T¯i are the prop- agator matrices from the bottom layer to the i-th layer. As the propagator matrices are normal, it can be shown that |Eui |2 = |Euo |2 , (5.53) which is in agreement with the assumption of total transmission. Similarly, taking advan- tage of the bidirectionality of the propagator matrices, the field components of a downgoing plane wave propagating in the −kˆ direction (from the i-th layer to the ground) can be written as Edo = T¯1T¯2 · · · T¯iEdi , (5.54) where Edi is the field at the top of the i-th layer. In radar experiments, the transverse field component of the signal backscattered from 64 a radar range ri = i∆r can be modeled as Ero ∝ κi T¯1T¯2 · · · T¯iT¯i · · · T¯2T¯1︸ ︷︷ ︸ Π¯i Eto, (5.55) where Eto and Ero are the fields transmitted and received by the radar antenna in the kˆ direction. Above, Π¯i denotes a two-way propagator matrix that accounts for the polariza- tion effects on the waves incident on and backscattered from the radar range ri (upgoing and downgoing waves, respectively). In addition, κi is a scattering factor that depends on the radar cross section (RCS) of the scatterers at the range ri (scatterers that in the case of the ionosphere are randomly moving electrons). As an example, consider the following incoherent scatter radar configuration. Two or- thogonal linearly polarized antennas with very narrow beams (pencil beams) are placed at Jicamarca. The antennas are steerable in the north-south plane such that different mag- netic aspect angle directions can be probed by the radar beams. Considering a coordinate system in which the x-, y-, and z-axes are parallel to the east, north, and zenith directions, we can express the radar pointing direction as kˆ = cosψ zˆ + sinψ yˆ, where ψ is the zenith angle. Furthermore, the radar antennas are oriented such that their polarization vectors are given by pˆ1 = xˆ and pˆ2 = − sinψ zˆ + cosψ yˆ, (5.56) respectively. Note that pˆ1, pˆ2, and kˆ constitute a set of orthogonal unit vectors. Operating at 50 MHz, the antenna beams scan the ionosphere from north to south to build power maps of the backscattered signals. In every pointing direction, narrow pulses are transmitted so that range filtering effects (due to the convolution of the pulse shape with the response of the ionosphere) can be ignored. In transmission, only the first antenna is excited, while, in reception, both antennas are used to collect the backscattered signals. As a result, the co- and cross-polarized antenna voltages sampled at each range ri can be expressed as v1(kˆ) ∝ κi pˆT1 Π¯i pˆ1 and v2(kˆ) ∝ κi pˆT2 Π¯i pˆ1, (5.57) 65 0 2 4 6 0 100 200 300 400 500 600 700 800 900 N e [1011 m−3] H ei gh t [k m] 1 1.5 2 0 100 200 300 400 500 600 700 800 900 T e /Ti H ei gh t [k m] Figure 5.3 Electron density and Te/Ti profiles as functions of height. where the two-way propagator matrix Π¯i (defined above) is dependent on the electron density and magnetic field values along kˆ up to the range ri. As the scattering factor κi is a random variable, statistics of the collected voltages have to be estimated in order to analyze the characteristics of the backscattered signals. Thus, we can model the mean square values of v1 and v2 as 〈|v1|2〉 ∝ σv Γ1 and 〈|v2|2〉 ∝ σv Γ2, (5.58) where σv = 〈|κi|2〉 is the volumetric RCS of the medium, which is dependent on the electron density, temperature ratio, and magnetic aspect angle at any given range. In addition, Γ1 and Γ2 are polarization coefficients defined as Γ1 = ∣∣∣pˆT1 Π¯i pˆ1∣∣∣2 and Γ2 = ∣∣∣pˆT2 Π¯i pˆ1∣∣∣2 . (5.59) Note that, given the orthogonality properties of the polarization vectors pˆ1 and pˆ2, it can be verified that the polarization coefficients in (5.59) satisfy the relation Γ1 + Γ2 = 1. To simulate the radar voltages, an ionosphere with the electron density and Te/Ti profiles displayed in Figure 5.3 was considered. The density profile was obtained from 66 ionosonde measurements conducted at Jicamarca on June 8, 2004 (around local noon). The Te/Ti profile is modeled as a layer of Gaussian shape centered at an altitude of 250 km with a width of 120 km, the bottom and peak values of Te/Ti are 1 and 2, respectively. In addition, the magnetic field above Jicamarca was computed using the International Geomagnetic Reference Field (IGRF) model [e.g., Olsen et al., 2000] for the same date of the density profile. Let us first analyze magneto-ionic propagation effects on the radar voltages, disregard- ing scattering effects. For this purpose, polarization coefficients Γ1 and Γ2 are displayed in Figure 5.4 as functions of distance and altitude from Jicamarca (in the plots, the positive horizontal axis is directed north). We can also think of these plots as the signal power detected by the radar antennas in the case that the RCS of the medium was constant at all heights. Note that, at low altitudes, where there is no ionosphere, signal returns will be detected only by the co-polarized antenna (i.e., by the same antenna used on transmis- sion). However, as the signal propagates farther through the ionosphere, magneto-ionic effects start taking place. We can appreciate that, for most of the propagation directions, the polarization vector of the detected field rotates such that signal from one polarization goes to the other as the radar range increases (Faraday rotation effect). Note, however, that there is a direction in which the wave polarization does not rotate much. In this direction, the antenna beams are pointed perpendicular to the Earth’s magnetic field, and it can be observed that the polarization of the detected fields varies progressively from linear to circular as a function of height (Cotton–Mouton effect). Finally, note that at higher altitudes, where the ionosphere vanishes, no more magneto-ionic effects take place, and the polarization of the detected signal approaches a final state. Next, scattering and propagation effects are considered in the simulation of the backscat- tered power collected by the pair of orthogonal antennas described above. The incoherent scatter volumetric RCS formulated in the previous chapter is used in the calculations. In Figure 5.5, the simulated co-polarized (left panel) and cross-polarized (middle panel) power data are displayed as functions of distance and altitude from Jicamarca. In addi- tion, the right panel depicts the total power detected by both antennas. Note that power 67 Figure 5.4 Polarization coefficients for the mean square voltages detected by a pair of orthogonal linearly polarized antennas placed at Jicamarca. The antennas have very nar- row beams and scan the ionosphere from north to south probing different magnetic aspect angle directions. Note that, for most pointing directions, the polarization of the detected fields rotates (Faraday rotation effect), except in the direction where the beam is pointed perpendicular to B, in which case, the type of polarization changes from linear to circular (Cotton-Mouton effect). Figure 5.5 Co-polarized (left panel), cross-polarized (middle panel), and total (right panel) backscattered power detected by a pair of orthogonal linearly polarized antennas (see caption of Figure 5.4). Power levels are displayed in units of electron density. In each plot, the dashed white lines indicate the directions half a degree away from perpendicular to B, while the continuous lines correspond to the directions one degree off. 68 levels are displayed as volumetric radar cross sections divided by 4pir2e (i.e., power levels are in units of electron density). In each plot, the dashed white lines indicate the directions half a degree away from perpendicular to B, while the continuous lines correspond to the directions one degree off. In the plots, we can observe that there is negligible backscattered power at low altitudes. In real measurements, however, coherent echoes stronger than the incoherent scatter signal can be detected at these altitudes. For instance, echoes from the electrojet at 100 km altitude or from the so-called 150-km region can be measured daily at Jicamarca (note that these echoes are only typical of the equatorial ionosphere). In our model, the presence of this type of echo is ignored because (to our knowledge) there is no accurate physical model for their radar cross sections. In the data analysis presented in the next chapter, we treat these echoes as clutter. At higher altitudes between approximately 200 and 700 km (where polarization effects are significant), co- and cross-polarized power maps exhibit features that are similar to the ones observed in Figure 5.4. Note, however, that there is an enhancement of the detected power in the direction where the antenna beams are pointed perpendicular to B; this can be observed more clearly in the plot of the total power (right panel of Figure 5.5). This feature is characteristic of the incoherent scatter process for probing directions perpendicular to B and for heights where electron temperature exceeds the ion temperature (i.e., Te > Ti) as described in the previous chapter. At even higher altitudes, scattered signals become weaker and weaker as the ionospheric electron density vanishes. Power measurements similar to the simulation results presented above were conducted for the first time a few years ago. Operating at a different frequency and located at a different latitude, the UHF (422 MHz) ALTAIR radar was used to scan the low-latitude ionosphere from north to south [e.g., Milla and Kudeki , 2006; Kudeki et al., 2006]. The radar is located on Roi-Namur (an island of the Kwajalein Atoll in the Pacific sector), and the measurements were carried out as part of the EQUIS 2 NASA rocket campaign in September 2004 [e.g., Lehmacher et al., 2006; Friedrich et al., 2006]. The antenna system at ALTAIR is dual-polarized; however, power signal was detected only by the 69 Figure 5.6 An example of the ALTAIR scan measurements collected during the 2004 EQUIS 2 campaign. In colors, levels of backscattered power calibrated to match electron densities in the probed region are displayed. co-polarized antenna because magneto-ionic propagation effects are negligible at UHF fre- quencies (for the range of ionospheric altitudes probed in these experiments). An example of the ALTAIR scan measurements is presented in Figure 5.6; note the similarity with our simulation results (right panel of Figure 5.5). At Jicamarca, however, the antenna array has limited steering capabilities, which makes scan measurements of the ionosphere not possible in practice. But instead, the modularity of the array permits the simultaneous observation of the ionosphere with different sections of the antenna phased to point into different directions. Multi-beam radar configurations are currently under test and evalu- ation at Jicamarca. The goal of these configurations is to make use of the aspect angle dependence of the incoherent scatter signal in the data inversion of ionospheric densities, temperatures, and drifts, all at the same time. The results obtained using one of these configurations will be presented in the next chapter of this dissertation. The data model is based on the propagation model presented in this section; however, an extra level of complexity has to be considered because, within the range of aspect angles illuminated by 70 the antenna beams, propagation and scattering effects vary quite rapidly. For this rea- son, the measured backscattered signals were carefully modeled taking into account the shapes of the antenna beams. A model for the beam-weighted incoherent scatter spectrum that considers magneto-ionic propagation and collisional effects is formulated in the next section. 5.3 Soft-Target Radar Equation and Magneto-Ionic Propagation In this section, the soft-target radar equation derived in Chapter 2 is reformulated using the wave propagation model described above. Consider a radar system composed of a set of antenna arrays (located in the same area) with matched filter receivers connected to the antennas used in reception. The mean square voltage at the output of the i-th receiver can be expressed as 〈|vi(t)|2〉 = EtKi ˆ dr dΩ dω 2pi |T (rˆ)|2 |Ri(rˆ)|2 k2 Γi(r) |χ(t− 2rc , ω)|2 4pir2 σv(k, ω), (5.60) where t is the radar delay, Et is the total energy of the transmitted radar pulse, and Ki is the i-th calibration constant (a proportionality factor that accounts for the gains and losses along the i-th signal path). Integrals are taken over range r, solid angle Ω, and Doppler frequency ω/2pi. In addition, k = −2korˆ denotes the relevant Bragg vector for a radar with a carrier wavenumber ko and associated wavelength λ. Above, T (rˆ) and Ri(rˆ) are the antenna factors of the arrays used in transmission and reception. Note that |T (rˆ)|2 and |Ri(rˆ)|2 are antenna gain patterns and the product |T (rˆ)|2 |Ri(rˆ)|2 is the corresponding two-way radiation pattern. The polarization coefficient Γi(r) is defined as Γi(r) = ∣∣∣pˆTi Π¯(r) pˆt∣∣∣2 , (5.61) where pˆt and pˆi are the polarization unit vectors of the transmitting and receiving anten- nas, and Π¯(r) is the two-way propagator matrix for the wave field components propagating 71 along rˆ (incident on and backscattered from the range r). Note that pˆt and pˆi are normal to rˆ, but they are not necessarily perpendicular to each other. In addition, χ(t, ω) is the radar ambiguity function and σv(k, ω) is the volumetric RCS spectrum, functions that have been defined in previous chapters. Similarly, the cross-correlation of the voltages at the outputs of the i-th and j-th receivers can be expressed as 〈vi(t)v∗j (t)〉 = EtKi,j ˆ dr dΩ dω 2pi |T (rˆ)|2Ri(rˆ)R∗j (rˆ) k2 Γi,j(r) |χ(t− 2rc , ω)|2 4pir2 σv(k, ω), (5.62) where Ki,j is a cross-calibration constant (dependent on gains and losses along the i-th and j-th signal paths), and Γi,j(r) is a cross-polarization coefficient defined as Γi,j(r) = ( pˆTi Π¯(r) pˆt )( pˆTj Π¯(r) pˆt )∗ . (5.63) Note that dispersion of the pulse shape due to wave propagation effects has been neglected in our model. Denoting by Si(ω) the self-spectrum of the signal at the output of the i-th receiver and applying Parseval’s theorem, we have that 〈|vi(t)|2〉 = ˆ dω 2pi Si(ω). (5.64) Likewise, the cross-spectrum Si,j(ω) and the cross-correlation of the signals at the outputs of the i-th and j-th receivers are related by 〈vi(t)v∗j (t)〉 = ˆ dω 2pi Si,j(ω). (5.65) Assuming that the ambiguity function is almost flat within the bandwidth of the RCS spectrum σv(k, ω) (which is a valid approximation in the case of short-pulse radar appli- cations), we can use equations (5.60) and (5.62) to obtain the following beam-weighted spectrum and cross-spectrum models: Si(ω) = EtKi ˆ dr |χ(t− 2rc )|2 4pir2 ˆ dΩ |T (rˆ)|2 |Ri(rˆ)|2 k2 Γi(r)σv(k, ω) (5.66) 72 and Si,j(ω) = EtKi,j ˆ dr |χ(t− 2rc )|2 4pir2 ˆ dΩ |T (rˆ)|2Ri(rˆ)R∗j (rˆ) k2 Γi,j(r)σv(k, ω), (5.67) where χ(t) = 1 T f∗(−t) ∗ f(t) (5.68) is the normalized auto-correlation of the pulse waveform f(t). In the radar equations (5.66) and (5.67), the polarization coefficients Γi(r) and Γi,j(r) effectively modify the radiation patterns; thus, the spectrum shapes are dependent not only on the scattering process but also on the modes of propagation. This dependence further complicates the spectrum analysis of radar data and the inversion of physical parameters. 73 CHAPTER 6 MAGNETO-IONIC EFFECTS ON SIMULATED AND MEASURED ISR DATA In this chapter, signatures of magneto-ionic propagation effects on incoherent scatter radar measurements are analyzed and discussed in detail. This analysis is performed us- ing simulated and measured data of the differential-phase technique used at Jicamarca to measure F-region electron densities with antenna beams pointed perpendicular to the geomagnetic field. First, the radar configuration of the differential-phase technique is de- scribed. Next, magneto-ionic effects on the differential-phase antenna patterns are modeled and simulated using the radar equations of Section 5.3. In the simulation, the ionospheric plasma configuration specified in Section 5.2 is considered. The patterns are then utilized to compute IS signal power and cross-correlation profiles. The different features of the computed profiles are discussed and linked to the altitude variations of the simulated ra- diation patterns. The computed profiles are then compared to real measurements carried out at Jicamarca. The results of the comparisons show good agreement between mod- eled and measured data. Finally, magneto-ionic effects on spectrum and cross-spectrum measurements are identified and analyzed based on our simulation results. 6.1 The Differential-Phase Radar Configuration The differential-phase technique developed for the Jicamarca incoherent scatter radar system is used to measure high-precision drifts and ionospheric densities simultaneously at F-region heights [e.g., Kudeki et al., 2003]. In this mode of operation, the Jicamarca antenna is phased to point perpendicular to the Earth’s magnetic field. The Jicamarca 74 antenna consists of two superimposed square arrays of half-wave dipoles orthogonal to each other with a diagonal oriented ∼ 6◦ east of the geographic north [Ochs, 1965]. The arrays, labeled as “up” and “down,” are northeast and southeast polarized. As shown in the antenna diagram of Figure 6.1, the up and down arrays are subdivided into quarters. In transmission, all quadrants of the down (southeast polarized) array are excited. In reception, backscattered fields detected by the up and down polarizations of the north and south antenna quarters are combined using RF hybrids to synthesize meridional and zonal polarized field components. After matched filter detection, the radar signals are digitalized and stored in a computer. Further details of the differential-phase technique are given by Feng et al. [2003, 2004]. The four antenna channels of the differential-phase radar configuration are: • Channel 1: Meridional polarization of the north antenna quadrant. • Channel 2: Zonal polarization of the north antenna quadrant. • Channel 3: Meridional polarization of the south antenna quadrant. • Channel 4: Zonal polarization of the south antenna quadrant. Statistics of the signals detected by these channels, spectrum and cross-spectrum measure- ments, can be modeled using the soft-target radar equations of Section 5.3. Assuming that very short pulses are transmitted, the spectra measured by the i-th channel can take the following form Si(ω) = EtKiδr 4pir2 ˆ dΩWi(r)σv(k, ω), (6.1) where δr is the effective radar pulse width defined as δr ≡ ´ dr |χ(t− 2rc )|2. Similarly, the cross-spectra of the i-th and j-th channels can be modeled as Si,j(ω) = EtKi,jδr 4pir2 ˆ dΩWi,j(r)σv(k, ω). (6.2) In these equations, Wi(r) and Wi,j(r) denote effective two-way radiation patterns that take into account the magneto-ionic propagation effects on the radar signals. Due to these effects, the effective patterns vary smoothly as functions of the radar range. Note that 75 ND H1 NU SU TX RX DN -p ol UP -p ol M ag ne tic No rth SD H2 Ch 1 Ch 2 Ch 3 Ch 4 I+ I- O- O+ I+ I- O- O+ F ig ur e 6. 1 A nt en na di ag ra m of th e di ffe re nt ia l-p ha se te ch ni qu e de ve lo pe d fo r th e Ji ca m ar ca in co he re nt sc at te r ra da r sy st em . In th is m od e of op er at io n, th e Ji ca m ar ca an te nn a is ph as ed to po in t pe rp en di cu la r to B . In tr an sm is si on ,t he so ut he as t po la ri ze d do w n an te nn a qu ar te rs ar e ex ci te d. In re ce pt io n, ba ck sc at te re d fie ld s de te ct ed by th e up an d do w n po la ri za ti on s of th e no rt h an d so ut h qu ar te rs ar e co m bi ne d to sy nt he si ze m er id io na la nd zo na lp ol ar iz ed fie ld co m po ne nt s. 76 the solid angle integrations in (6.1) and (6.2) need to be carefully evaluated because the incoherent scatter RCS spectrum σv(k, ω) varies quite rapidly at small magnetic aspect angles, the angular regime that is illuminated by the radar beams. The effective radiation patterns to be considered in the self-spectrum model (6.1) are given next. For the meridional polarized antenna configurations, channels 1 and 3, the patterns are defined as W1(r) = W3(r) = 1 k2 Gt(rˆ)Gr(rˆ) Γm(r), (6.3) while, for the zonal polarized antennas, channels 2 and 4, we have that W2(r) = W4(r) = 1 k2 Gt(rˆ)Gr(rˆ) Γz(r). (6.4) In these expressions, k is the corresponding Bragg wavenumber (k = 4pi/λ) that normal- izes the solid angle integration of the two-way patterns to units of backscatter aperture (i.e., units of area). Above, Gt(rˆ) and Gr(rˆ) are the free-space radiation patterns of the antenna arrays used in transmission and reception (full and quarter arrays, respectively). In addition, the polarization coefficients of the meridional and zonal patterns, Γm(r) and Γz(r), are defined as Γm(r) = ∣∣∣pTm Π¯(r) pˆt∣∣∣2 and Γz(r) = ∣∣∣pTz Π¯(r) pˆt∣∣∣2 , (6.5) where pˆt is the polarization unit vector of the field transmitted in the direction rˆ, and pm and pz are the polarization vectors of the meridional and zonal components of the scattered fields detected at the antenna level. In addition, Π¯(r) denotes the two-way propagator matrix defined in Section 5.2. In the Jicamarca frame of reference — a coordinate system in which the antenna is parallel to the xy-plane and the dipoles of the down- and up-arrays are aligned with the x- and y-axes — the polarization unit vectors of the down and up antennas can be expressed 77 as pˆd = rˆ× rˆ× xˆ |ˆr× rˆ× xˆ| and pˆu = rˆ× rˆ× yˆ |ˆr× rˆ× yˆ| (6.6) for fields propagating along rˆ. Since, in transmission, the down-array is excited, we have that pˆt = pˆd, (6.7) while, in reception, the meridional and zonal polarization vectors are given by pm = 1√ 2 (pˆu − pˆd) and pz = 1√ 2 (pˆu + pˆd) . (6.8) To model the cross-spectrum measurements, expressions for the effective two-way ra- diation patterns corresponding to each of the six pairs of antenna configurations are given next. For the cross-spectra between cross-polarized channels of the same antenna quad- rants, the effective patterns for the north pair (channels 1 and 2) and the south pair (channels 3 and 4) are equal to W1,2(r) = W3,4(r) = 1 k2 Gt(rˆ)Gr(rˆ) Γm,z(r), (6.9) where the polarization coefficient Γm,z is defined as Γm,z(r) = ( pˆTm Π¯(r) pˆt )( pˆTz Π¯(r) pˆt )∗ . (6.10) In the case of the cross-spectra between co-polarized channels of different antenna quad- rants, the pattern for the pair of meridional antenna channels takes the form W1,3(r) = ejkrˆ·D k2 Gt(rˆ)Gr(rˆ) Γm(r), (6.11) while, for the pair of zonal channels, the pattern is equal to W2,4(rˆ) = ejkrˆ·D k2 Gt(rˆ)Gr(rˆ) Γz(r), (6.12) 78 where D is a distance vector from the center of the south quadrant to the center of the north quadrant. Finally, we have the case of the cross-spectra between cross-polarized channels of different antenna quadrants. The pattern for the pair composed of the meridional channel of the north quadrant and the zonal channel of the south quadrant is given by W1,4(r) = ejkrˆ·D k2 Gt(rˆ)Gr(rˆ) Γm,z(r), (6.13) and the pattern for the pair composed of the zonal channel of the north quadrant and the meridional channel of the south quadrant is equal to W2,3(r) = ejkrˆ·D k2 Gt(rˆ)Gr(rˆ) Γ ∗ m,z(r). (6.14) In Figure 6.2, we present the two-way radiation pattern Gt(rˆ)Gr(rˆ) of the differential- phase radar configuration in free space.1 In the left panel, the pattern (in dB units relative to its peak value) is displayed as a function of θx and θy (direction cosines with respect to the x- and y-axes of the Jicamarca frame of reference). In the right panel, the pattern is displayed as a contour plot in a rotated coordinate system in which the x- and y-axes are now aligned with the east and north directions (the Jicamarca antenna frame of reference is rotated 51◦ positive azimuth). In the plots, the continuous black lines trace the positions where the propagation directions are perpendicular to the geomagnetic field at radar ranges r = 100, 500, 1000 km. 6.2 Magneto-Ionic Effects on the Radiation Patterns Simulation results showing the altitudinal variation of the radiation patterns of the differential-phase radar configuration caused by the magneto-ionic propagation of the radar signals are presented next. In the simulation, the electron density and Te/Ti profiles dis- played in Figure 5.3 are considered. The Earth’s magnetic field components for June 2004 are calculated using the IGRF model. The simulated effective radiation patterns are dis- 1The Jicamarca antenna radiation patterns can be efficiently computed using a 2D FFT-based algorithm that takes advantage of the spatial uniformity of the antenna array. 79 θ x θ y D 2W = 7 9. 4 dB − A BS = 5 97 1. 2 − 0. 05 0 0. 05 − 0. 08 − 0. 06 − 0. 04 − 0. 020 0. 02 0. 04 0. 06 0. 08 − 40 − 35 − 30 − 25 − 20 − 15 − 10 − 50 10 0 km 50 0 km 10 00 k m θ x θ y D 2W = 7 9. 4 dB − A BS = 5 97 1. 2 − 0. 05 0 0. 05 − 0. 08 − 0. 06 − 0. 04 − 0. 020 0. 02 0. 04 0. 06 0. 08 − 40 − 35 − 30 − 25 − 20 − 15 − 10 − 50 10 0 km 50 0 km 10 00 k m F ig ur e 6. 2 T w o- w ay an te nn a ra di at io n pa tt er n of th e di ffe re nt ia l-p ha se ra da r co nfi gu ra ti on (i n fr ee -s pa ce ). In th e le ft pa ne l, th e ra di at io n pa tt er n (i n dB un it s re la ti ve to it s pe ak va lu e) is di sp la ye d as a fu nc ti on of di re ct io n co si ne s w it h re sp ec t to th e x - an d y -a xe s of th e Ji ca m ar ca fr am e of re fe re nc e. In th e ri gh t pa ne l, th e pa tt er n is di sp la ye d as a co nt ou r pl ot in a ro ta te d co or di na te sy st em in w hi ch th e x - an d y -a xe s ar e no w al ig ne d w it h th e ea st an d no rt h di re ct io ns (t he Ji ca m ar ca fr am e of re fe re nc e is ro ta te d 51 ◦ po si ti ve az im ut h) . 80 played in Figures 6.3, 6.4, and 6.5 in the rotated coordinate system specified above. Each row of plots corresponds to different radar ranges. In the left and right panels, the merid- ional and zonal polarized radiation patterns are displayed as functions of direction cosines. In the middle panels, the polarization states of the wave field components corresponding to different propagation directions are shown. In Figure 6.3, we can observe that the initial polarization vector (at r = 0 km) is pointed southeast; therefore, both meridional and zonal field components are equally excited and the associated patterns look identical to each other. Note that the initial patterns are actually the radiation patterns in free-space (disregarding near field effects). In the first 100 km, the radiation patterns vary very little, but as soon as the ionosphere appears, magneto-ionic effects start taking place. At 200 km, Faraday rotation effects are noticeable at aspect angle directions around one degree (∼ 0.02 dc) to the north and south of the loci of perpendicularity (depicted as black lines). In the north, the polarization of the wave fields rotates counter-clockwise, while, in the south, the contrary happens. Because of this, we can observe that power from the northern part of the meridional pattern goes to the zonal pattern, while power from the southern part of the zonal pattern goes to the meridional one. Note, however, that in the close vicinity of the perpendicular-to-B direction, rotation does not take place. In Figure 6.4, beam shape differences between the meridional and zonal patterns be- come more evident. Note that at 250 km and 300 km, the sidelobes of the patterns have clearly changed their shapes. At higher altitudes around 350 km and 400 km (close to the peak of the density profile), the polarization of the backscattered fields at one degree from the north and south of the loci of perpendicularity have rotated about 90 degrees. Also note that around the perpendicular direction, the polarization states have changed from linear to elliptical. Moreover, we can see that the shapes of the mainlobes have started to look different. Finally, in Figure 6.5, we can observe that not only have the mainlobe shapes become different but their centers have also been displaced. Note that the zonal beam pattern is now pointing to the north of the locus of perpendicularity, while the meridional beam 81 θ x θ y Meridional Beam Height = 0.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 0.0 km θ x Zonal Beam Height = 0.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 100.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 100.0 km θ x Zonal Beam Height = 100.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 150.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 150.0 km θ x Zonal Beam Height = 150.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 200.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 200.0 km θ x Zonal Beam Height = 200.0 km −0.02 0 0.02 −40 −30 −20 −10 0 Figure 6.3 Meridional and zonal effective radiation patterns for r = 0, 100, 150, 200 km. 82 θ x θ y Meridional Beam Height = 250.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 250.0 km θ x Zonal Beam Height = 250.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 300.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 300.0 km θ x Zonal Beam Height = 300.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 350.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 350.0 km θ x Zonal Beam Height = 350.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 400.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 400.0 km θ x Zonal Beam Height = 400.0 km −0.02 0 0.02 −40 −30 −20 −10 0 Figure 6.4 Same as Figure 6.3 but for r = 250, 300, 350, 400 km. 83 θ x θ y Meridional Beam Height = 450.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 450.0 km θ x Zonal Beam Height = 450.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 500.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 500.0 km θ x Zonal Beam Height = 500.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 550.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 550.0 km θ x Zonal Beam Height = 550.0 km −0.02 0 0.02 −40 −30 −20 −10 0 θ x θ y Meridional Beam Height = 600.0 km −0.02 0 0.02 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.02 0 0.02 θ x Polarization Vector (RX) Height = 600.0 km θ x Zonal Beam Height = 600.0 km −0.02 0 0.02 −40 −30 −20 −10 0 Figure 6.5 Same as Figure 6.3 but for r = 450, 500, 550, 600 km. 84 pattern is still pointing, approximately, perpendicular to B. As the altitude increases and the electron density vanishes, we can appreciate that magneto-ionic propagation ef- fects are less severe, and the radiation patterns eventually relax toward their final shapes. Note that the final polarization states have changed a lot compared to their initial linear state. At 600 km, for instance, the polarization of the wave field components propagating perpendicular to B have become more or less circular. We can also observe at the same altitude, but one degree to the north or south of the perpendicular direction, that the polarization states have remained linear but rotated many degrees. The behavior of the radiation patterns described above has an impact on the measured signals, as we discuss next. 6.3 Simulated and Measured Power and Cross-Correlation Profiles The effective radiation patterns simulated in the previous section are used to compute signal power and cross-correlation profiles. Based on the soft-target radar equations (6.1) and (6.2), we can model the power measured by the i-th channel as follows, 〈|vi|2〉 = ˆ dωSi(ω) = EtKiδr 4pir2 ˆ dΩWi(r)σv(k), (6.15) while the cross-correlation between the i-th and j-th channels can be expressed as 〈viv∗j 〉 = ˆ dωSi,j(ω) = EtKi,jδr 4pir2 ˆ dΩWi,j(r)σv(k), (6.16) where Wi(r) and Wi,j(r) are the effective two-way radiation patterns, and σv(k) is the volumetric IS radar cross section (see Chapter 4). As discussed in previous chapters, the incoherent scatter RCS is a function of the angle between the wave propagation direction k and the ambient magnetic field B. At the magnetic equator, it can be verified that σv(k) varies progressively in the north-south plane (which is nearly parallel to the geomagnetic field lines). In fact, in the vicinity of the 85 direction perpendicular to B, the magnitude of σv(k) varies even more rapidly. However, in the east-west plane, σv(k) changes very slowly. To take advantage of this behavior in the numerical evaluation of the solid angle integrals in (6.15) and (6.16), we consider a coordinate system in which the x- and y-axes are aligned with the east and north directions of the Jicamarca location (see, for instance, the right panel of Figure 6.2). After recasting the integrals in terms of θx and θy (direction cosines with respect to the x- and y-axes of the coordinate system considered), we compute (6.15) and (6.16) applying a finite-element integration method. At every radar range, values of σv(k) are finely sampled in the θy direction but very coarsely in the θx direction. Then, the integrals are computed by adding the samples of σv(k) weighted by numerical coefficients that account for the shapes of the effective radiation patterns. Signal power and cross-correlation profiles are simulated for the differential-phase radar configuration of June 8, 2004, considering the set of ionospheric plasma parameters given in the previous chapter. As mentioned before, four reception channels are used in the differential-phase experiments. In the left panel of Figure 6.6, the simulated power profiles of the four antenna channels are plotted in arbitrary units. Channels 1 and 2 (blue and green curves) correspond to the meridional and zonal polarized signals detected by the north quadrant of the Jicamarca array, while channels 3 and 4 (red and light-blue curves) correspond to the signals detected by the south quadrant. Note that the meridional power profiles (blue and red curves) overlap each other. In the case of the zonal profiles (green and light-blue curves), the same happens. However, both the meridional and zonal profiles are different from each other. This dissimilarity is an effect caused by the magneto- ionic propagation of the radar signals in the ionosphere. According to the simulation results presented in the previous section, the shapes of the effective radiation patterns and, therefore, their corresponding radar volumes vary as functions of the radar range. These variations, however, are different for different antenna polarizations. Therefore, non-co-polarized antennas, such as the ones used by the differential-phase technique, detect distinct signal power profiles. The simulated cross-correlations of the antenna signals are used to compute coherence profiles for different pairs of antenna channels. Coherence is 86 0 2 4 0 100 200 300 400 500 600 700 800 Signal R an ge [k m] [1] North M−pol [2] North Z−pol [3] South M−pol [4] South Z−pol 0 0.25 0.5 0.75 1 0 100 200 300 400 500 600 700 800 Coherence [1,2] [3,4] [1,3] [2,4] −180 −90 0 90 180 0 100 200 300 400 500 600 700 800 Phase [deg] [1,2] [3,4] [1,3] [2,4] Figure 6.6 Simulated incoherent scatter signal power and coherence (magnitude and phase) profiles for the differential-phase radar configuration of June 8, 2004. defined as ρij = 〈viv∗j 〉√〈|vi|2〉〈|vj |2〉 . (6.17) The magnitude and phase of the computed coherence profiles are displayed in the middle and right panels of Figure 6.6. The coherences between the cross-polarized signals of the north (channels 1 and 2) and of the south (channels 3 and 4) antenna quadrants are plotted as blue and green curves. As expected, these coherence profiles overlap each other. At low altitudes, where magneto-ionic effects are negligible, the coherence is effectively one, since the beam shapes of the meridional and zonal polarized antenna arrays are still the same. However, it can be observed that at 200 km the coherence magnitudes start to decrease. This happens because, as the ionospheric density increases, the meridional and zonal radi- ation patterns change their shapes due to magneto-ionic propagation effects. Eventually, the patterns become somewhat different from each other, which leads to a reduction of the magnitude of the coherence. The coherences between the co-polarized signals of different antenna quadrants, meridional signals (channels 1 and 3) and zonal signals (channels 2 87 and 4), are plotted as red and light-blue curves. Note that the magnitude and phase of the coherence profiles vary as functions of height differently for distinct polarizations. We can observe, once more, that at low altitudes the coherence profiles look similar. As alti- tude increases, note in both cases that the magnitude of the coherence first rises, which implies that the radiation patterns have become more compact. At higher altitudes, the coherence of the meridional polarized signals decreases, while that of the zonal signals increases. The reason for this is also related to the altitudinal variations of the effective radiation patterns. While in the case of the meridional beams, the patterns become wider (power is spread over aspect angle), in the case of the zonal beams, the patterns become more compact (power is concentrated in some region). The phase of the coherence profiles brings further information regarding the effective centers of the beam shapes. Note that, when there is no ionosphere, the meridional and zonal beams are pointing to the same direction (which is slightly to the north of the loci of perpendicularity and corresponds to a coherence phase of about 90◦). In the case of the meridional coherence, we can observe that the phase (red curve) is decreasing. This is in agreement with the fact that, as the radar range increases, the centers of the meridional patterns move to the south of the loci of perpendicularity. This effect can be observed in the simulated patterns displayed in the previous section. In the case of the zonal coherence, the contrary happens, i.e., the phase (light-blue curve) increases. This is also expected because, at ionospheric altitudes, the zonal polarized patterns are pointing to the north of their initial pointing direction. In addition to the simulation results, we present in Figure 6.7, measured power and coherence profiles obtained after a five-minute integration of incoherent radar signals col- lected during the differential-phase experiment conducted on June 8, 2004. We can observe that the simulated profiles resemble very closely the actual measurements. This gives us confidence that our incoherent scatter model is capable of reproducing the different features of real measurements. In the following section, the characteristics of spectrum and cross- spectrum measurements are analyzed based on the simulation results presented above. 88 0 100 200 300 0 100 200 300 400 500 600 700 800 Signal R an ge [k m] 0 0.25 0.5 0.75 1 0 100 200 300 400 500 600 700 800 Coherence 2004−06−08 14:30:00 −180 −90 0 90 180 0 100 200 300 400 500 600 700 800 Phase [deg] [1] North M−pol [2] North Z−pol [3] South M−pol [4] South Z−pol [1,2] [3,4] [1,3] [2,4] [1,2] [3,4] [1,3] [2,4] Figure 6.7 Measured incoherent scatter signal power and coherence (magnitude and phase) profiles that correspond to the differential-phase experiment conducted on June 8, 2004. 6.4 Magneto-Ionic Effects on Spectrum and Cross-Spectrum Measurements In this section, some of the features of incoherent scatter spectrum and cross-spectrum measurements are analyzed. The sample spectra considered in this analysis were acquired during the differential-phase experiment conducted at Jicamarca on June 8, 2004. In the left panel of Figure 6.8, meridional and zonal polarized self-spectrum measurements are dis- played as functions of Doppler velocity (or frequency) for different radar ranges. In general, the incoherent scatter spectrum measured with perpendicular-to-B antennas is composed of a narrow spectral component with a sharp peak and a wide spectral component that is typically aliased in pulse-to-pulse spectrum observations. We can appreciate these fea- tures of the measured spectra in the plots of Figure 6.8. Note that at low altitudes, where the ionosphere and the magneto-ionic propagation effects are negligible, the four spectra have the same shape. However, as altitude increases, we can observe that the meridional 89 IS R D op pl er Sp ec tr a − 20 04 − 06 − 08 14 :3 0: 00 IS R C oh er en ce Sp ec tr a − 20 04 − 06 − 08 14 :3 0: 00 N or th − So ut h ba se lin e in te rf er om et er 0 0. 01 0. 02 R an ge = 5 28 k m M ag ni tu de 0 0. 01 0. 02 R an ge = 4 52 k m 0 0. 01 0. 02 R an ge = 3 76 k m 0 0. 01 0. 02 R an ge = 3 00 k m − 15 0 − 10 0 − 50 0 50 10 0 15 0 0 0. 01 0. 02 R an ge = 2 28 k m Ve lo ci ty [m /s] M −p ol (N ort h) Z− po l (N ort h) M −p ol (S ou th) Z− po l (S ou th) 0 0. 2 0. 4 0. 6 0. 81 R an ge = 5 28 k m M ag ni tu de − 18 0 − 9009018 0 Ph as e [de g] 0 0. 2 0. 4 0. 6 0. 81 R an ge = 4 52 k m − 18 0 − 90090 0 0. 2 0. 4 0. 6 0. 81 R an ge = 3 76 k m − 18 0 − 90090 0 0. 2 0. 4 0. 6 0. 81 R an ge = 3 00 k m − 18 0 − 90090 − 10 0 0 10 0 0 0. 2 0. 4 0. 6 0. 81 R an ge = 2 28 k m Ve lo ci ty [m /s] − 10 0 0 10 0 − 18 0 − 90090 Ve lo ci ty [m /s] M −p ol (N × S* ) Z− po l (N × S* ) F ig ur e 6. 8 In co he re nt sc at te r D op pl er sp ec tr a (l ef t) an d no rt h- so ut h ba se lin e co he re nc e sp ec tr a (r ig ht ) m ea su re d w it h th e di ffe re nt ia l- ph as e ra da r te ch ni qu e at Ji ca m ar ca on Ju ne 8, 20 04 . 90 spectra become sharper than the zonal spectra. This behavior is related to the altitudi- nal variation of the effective radiation patterns presented earlier in this chapter. Due to magneto-ionic propagation effects, the centers of the radiation patterns are displaced at F-region heights. While the meridional pattern points approximately perpendicular to the magnetic field, the zonal pattern points to the north of the locus of perpendicularity. This implies that incoherent scatter signals coming from off-perpendicular directions contribute more to shaping the zonal spectra than the meridional spectra. As off-perpendicular signal contributions have wide spectral bandwidths, the zonal polarized signals have wider spec- tra than the meridional signals. Although the differences between the meridional and zonal spectra are noticeable, they are not very large and can be masked by the statistical uncer- tainties of the measurements. Because of this, a large number of collected spectra need to be averaged in order to reduce the uncertainties to a level that allows the identification of the spectrum differences. In the right panel of Figure 6.8, the magnitude and phase of the cross-spectra between co-polarized signals are displayed for different radar ranges. The blue curves correspond to the cross-spectra of the meridional signals (channels 1 and 3), while the green curves correspond to the cross-spectra of the zonal signals (channels 2 and 4). As expected, the meridional and zonal cross-spectra look the same at low altitudes; however, as altitude increases, the magnitudes and phases of both cross-spectra become different. Note, for instance, that at altitudes above ∼ 400 km, the meridional cross-spectra are still peaky, but the zonal spectra have become almost flat. This is related to the fact that, at these altitudes, the zonal radiation patterns point away from the perpendicular-to-B direction. Therefore, zonal signals measured by the radar are mostly composed of wide spectrum (off-perpendicular) components. Also note that the phases of the measured cross-spectra vary in altitude. In general, the zonal phase is greater than the meridional phase, which is related to the fact that the effective zonal beam points to the north of the meridional beam. In this chapter, we have developed an incoherent scatter model that takes into account magneto-ionic propagation. Simulated power and cross-correlation data have been com- 91 pared with real measurements. The results of the comparisons show excellent agreement between simulated and measured data. The comparison of spectrum measurements with simulated spectra obtained with our model will be the subject of future studies. Also part of these investigations will be the analysis of the sensitivity of the spectrum and cross-spectrum models to electron and ion temperatures. In the following chapter, our incoherent scatter model is applied to the inversion of ionospheric densities and Te/Ti pro- files using power and cross-correlation measurements obtained with a multi-beam radar configuration. 92 CHAPTER 7 THREE-BEAM EXPERIMENTS AT JICAMARCA In the present chapter, the incoherent scatter data model described in this dissertation is applied to the inversion of ionospheric plasma parameters using radar data collected in recent multi-beam experiments at Jicamarca. In these experiments, different sections of the Jicamarca antenna array are phased to point west, east, or south of the local zenith. While the west and east beams are aimed perpendicular to the magnetic field, the south beams are pointed in a direction ∼ 4◦ to the south of the loci of perpendicularity. The details of the radar experiments and the antenna configurations are given in Section 7.1. Statistics of the backscattered signals collected in these radar experiments are processed to estimate state parameters of the equatorial ionosphere. For instance, vertical and zonal plasma drifts are determined from the Doppler shifts of the frequency spectra measured with the antenna beams pointed perpendicular to B (the west and east beams). The spectral fitting technique developed to estimate the drifts was the subject of previous research and was reported a decade ago by Kudeki et al. [1999]. In this chapter, we describe a different inversion procedure developed for the estimation of electron densities and electron-to-ion temperature ratios. Using the signal power and cross-correlation data collected in these experiments (with all the antenna beams), we infer electron density profiles corrected for the Te/Ti variations throughout the bottom-side of the F-region. The estimation procedure exploits the magnetic aspect angle dependence of the incoherent scatter RCS at heights where the electron temperature exceeds the ion temperature (i.e., when Te > Ti). The development of this technique was motivated by the work of Milla and 93 Kudeki [2006]. Samples of spectrum and cross-spectrum measurements, as well as power and cross-correlation data are displayed in Section 7.2. The forward model of the radar data and the parameter estimation technique are outlined in Section 7.3. Electron density and Te/Ti maps resultant from the inversion of the incoherent scatter radar data collected in these multi-beam experiments are presented in the last section of this chapter. The goodness-of-fit criterion is used to analyze the quality of the fitting results. In addition, our density estimates are compared to ionosonde measurements carried out simultaneously at Jicamarca. 7.1 Experiment Description and Antenna Configurations Two types of multi-beam incoherent scatter radar experiments with antenna beams pointed into three different directions were carried out at Jicamarca in recent years. The experiments were referred to as the DEWD 3Ba and 3Bb modes (names that stand for Densities and East-West Drifts measured with three beams). In a first data acquisition campaign, the radar was operated in the 3Ba mode from June 16 to 19, 2008. Radar experiments in the 3Bb mode were conducted in two subsequent campaigns, from June 23 to 24, 2008, and also from January 9 to 12, 2009. As mentioned before, the modularity of the Jicamarca array was exploited in order to configure different sections of the antenna to point west, east, or south of the Jicamarca zenith. A total of six antenna beams were configured, two for each pointing direction. While the west and east beams are aimed perpendicular to the magnetic field, the south beams are pointed in a direction ∼ 4◦ to the south of the loci of perpendicularity. The antenna configurations used in these experiments are displayed in Figure 7.1. In the 3Ba mode, the Jicamarca antenna is configured as follows. • West Beam 1: West quadrant of the up-polarization. • West Beam 2: East quadrant of the up-polarization. • East Beam 1: West quadrant of the down-polarization. 94 TX RX DN -p ol UP -p ol M ag ne tic No rth 2 15 4 36 5 6 RX 1 W U: W es t 1 (Q ua rte r) RX 2 EU : W es t 2 (Q ua rte r) RX 5 SU -N U: S ou th (C o- Po l) RX 3 W D: E as t 1 (Q ua rte r) RX 4 ED : E as t 2 (Q ua rte r) RX 6 SD -N D: S ou th (X -P ol) TX RX DN -p ol UP -p ol M ag ne tic No rth 1 12 3 34 5 6 RX 1 EU -W U: W es t ( Bo wt ie) RX 2 NU : W es t ( Qu ar te r) RX 5 SU : S ou th (C o- Po l) RX 3 ED -W D: E as t ( Bo wt ie) RX 4 ND : E as t ( Qu ar te r) RX 6 SD : S ou th (X -P ol) F ig ur e 7. 1 A nt en na co nfi gu ra ti on s of th e 3B a (l ef t pa ne l) an d 3B b (r ig ht pa ne l) ex pe ri m en ts . 95 • East Beam 2: East quadrant of the down-polarization. • South Beam 1: South and north quadrants of the up-polarization. • South Beam 2: South and north quadrants of the down-polarization. In the 3Bb mode, the configuration is somewhat different. • West Beam 1: East and west quadrants of the up-polarization (bowtie antenna). • West Beam 2: North quadrant of the up-polarization. • East Beam 1: East and west quadrants of the down-polarization (bowtie antenna). • East Beam 2: North quadrant of the down-polarization. • South Beam 1: South quadrant of the up-polarization. • South Beam 2: South quadrant of the down-polarization. Note that the west beams are co-polarized in the northeast direction (up-polarization). The east beams are also co-polarized but in the southeast direction (down-polarization). However, the south beams are cross-polarized: the first beam is polarized northeast, while the second one is polarized southeast. The two-way radiation patterns of the radar beams used in the 3Ba and 3Bb ex- periments are presented in Figures 7.2 and 7.3. The patterns are displayed in colors as functions of direction cosines, θx and θy, with respect to the x- and y-axes of a coordinate system in which the antenna is located at the origin. In this coordinate system, the z-axis is perpendicular to the plane of the Jicamarca antenna, and the x- and y- axes are aligned with the local east and north directions (i.e., the coordinate system is rotated 51◦ positive azimuth with respect to the direction of alignment of the “down” antenna dipoles). In the plots, the black lines trace the loci of magnetic perpendicularity for June 2008 at different altitudes (100, 500, and 1000 km). As indicated earlier in this section, the west-beam antennas in the 3Ba mode are quarters of the Jicamarca array aligned in the east-west direction. On the other hand, in the 3Bb mode, the west-beam antennas have different shapes and are aligned in the north-south direction. One of them is still a quarter array, 96 3B a tw o- w ay an te nn a pa tt er ns 3B b tw o- w ay an te nn a pa tt er ns θ y W es t B ea m 1 (U p) D 2W = 7 3. 91 d B AB S = 28 33 .6 − C O H = 0. 52 − 0. 1 − 0. 050 0. 050. 1 W es t B ea m 2 (U p) D 2W = 7 3. 91 d B AB S = 28 33 .6 − C O H = 0. 52 − 40 − 30 − 20 − 100 θ x θ y Ea st B ea m 1 (D n) D 2W = 7 3. 84 d B AB S = 28 30 .4 − C O H = 0. 52 − 0. 1 − 0. 05 0 0. 05 0. 1 − 0. 1 − 0. 050 0. 050. 1 θ x Ea st B ea m 2 (D n) D 2W = 7 3. 84 d B AB S = 28 30 .4 − C O H = 0. 52 − 0. 1 − 0. 05 0 0. 05 0. 1 − 40 − 30 − 20 − 100 θ y W es t B ea m 1 (U p) D 2W = 7 8. 51 d B AB S = 46 90 .5 − C O H = 0. 56 − 0. 1 − 0. 050 0. 050. 1 W es t B ea m 2 (U p) D 2W = 7 5. 10 d B AB S = 30 94 .3 − C O H = 0. 56 − 40 − 30 − 20 − 100 θ x θ y Ea st B ea m 1 (D n) D 2W = 7 8. 09 d B AB S = 43 84 .4 − C O H = 0. 54 − 0. 1 − 0. 05 0 0. 05 0. 1 − 0. 1 − 0. 050 0. 050. 1 θ x Ea st B ea m 2 (D n) D 2W = 7 4. 59 d B AB S = 29 02 .4 − C O H = 0. 54 − 0. 1 − 0. 05 0 0. 05 0. 1 − 40 − 30 − 20 − 100 F ig ur e 7. 2 T w o- w ay ra di at io n pa tt er ns of th e w es t (t op ro w ) an d ea st (b ot to m ro w ) be am s us ed in th e 3B a (l ef t) an d 3B b (r ig ht ) ex pe ri m en ts . T he pa tt er ns ar e di sp la ye d as fu nc ti on s of di re ct io n co si ne s in a co or di na te sy st em in w hi ch x -a nd y -a xe s ar e al ig ne d w it h th e ea st an d no rt h di re ct io ns at Ji ca m ar ca . N ot e th at in th e 3B a m od e, th e pa tt er ns of th e w es t be am s ar e id en ti ca l, bu t in th e 3B b m od e, th e w es t be am s ha ve di ffe re nt sh ap es . T he sa m e ca n be ob se rv ed in th e ca se of th e ea st be am s. 97 3Ba mode 3Bb mode θ x θ y South Beam D2W = 79.7 dB − ABS = 6942 −0.1 −0.05 0 0.05 0.1 −0.1 −0.05 0 0.05 0.1 −40 −35 −30 −25 −20 −15 −10 −5 0 θ x θ y South Beam D2W = 73.8 dB − ABS = 4353 −0.1 −0.05 0 0.05 0.1 −0.1 −0.05 0 0.05 0.1 −40 −35 −30 −25 −20 −15 −10 −5 0 Figure 7.3 Two-way radiation patterns of the south beams used in the 3Ba (left) and 3Bb (right) experiments. Note that, in the north-south direction, the 3Ba south beam is narrower than the 3Bb beam. but the other is a combination of two quarters forming a “bowtie” antenna configuration (see Figure 7.1). Because of this, we can observe in the right panels of Figure 7.2 that the 3Bb west-beam radiation patterns look different from each other. The same can be observed in the case of the east-beam radiation patterns. As we can expect, the bowtie antennas have narrower beams in the east-west direction. Moreover, their radiation pat- terns have greater directivities and ABS factors, which implies that the data collected by the bowtie antennas will have better SNR values. In each plot, the values of the expected coherences between the co-polarized signals of the west-beam and east-beam antenna pairs are given. In general, the expected coherence is ∼ 0.5 (a reference value that is calculated disregarding magneto-ionic propagation effects). In Figure 7.3, the two-way radiation patterns of the south beams used in the 3Ba and 3Bb experiments are displayed. Note that, in the north-south direction, the 3Ba south beam is narrower than the 3Bb beam. This is because the south-beam antennas in the 3Ba mode are the combination of two quadrants of the Jicamarca array (the north and south quadrants), but in the 3Bb mode, the antennas are only single quadrants. In addition, note that the ABS factors of the south beams are in general greater than those of the west and east beams. The reason for this is that the south beams have lower sidelobes than the 98 Table 7.1 Radar operating parameters considered in the 3Ba and 3Bb experiments. Radar parameter West-East South Inter-pulse period 1000 km 1000 km Pulse width 45 km 45 km Pulse code 3-baud Barker 3-baud Barker Number of pulses 128 32 Sampling rate 5 km 5 km Number of samples 198 198 other beams. Despite the different antenna configurations, the radar system was operated in the 3Ba and 3Bb experiments with the same transmission and reception parameters (see Table 7.1). In transmission, electromagnetic wave pulses were emitted in sequences that alternate between perpendicular-to-B and off-perpendicular antenna beams. Each sequence started with 128 pulses transmitted with the west and east beams (the four beams were excited at the same time). These pulses were then followed by 32 pulses transmitted with the up-polarization of the south beams. The separation between pulses (inter-pulse period or IPP) was 6.667 ms, which is equivalent to a radar range of 1000 km. The length of each transmitted pulse was 45 km. To improve the range resolution, the radar pulses were coded with a 3-baud Barker sequence + +− that was flipped every other pulse (the baud length was 15 km). In reception, the backscattered signals coming from the ionosphere were sampled at the antenna outputs using digital receivers. The demodulation and the matched filter detection of the signals are performed numerically in these receivers. Samples of the in-phase and quadrature radar voltages were obtained every 5 km (i.e., the radar returns were oversampled by a factor of three with respect to the baud length). These samples were recorded and stored in a computer for later analysis using the data processing and inversion techniques described in the following sections. 99 7.2 Radar Measurements As mentioned before, observations of the upper atmosphere were conducted alternating between perpendicular-to-B and off-perpendicular radar beams. The perpendicular-to-B observations were carried out with the west and east beams, while the off-perpendicular observations were performed with the south beams. The antenna voltages recorded during each sequence of radar pulses constitute a block of data. In each block, 128 samples per range gate per channel were collected with the west and east beams (4 channels). In addition, 32 samples (also per range gate per channel) were acquired with the south beams (2 channels). In both cases, a total of 198 range gates were sampled. Using the fast Fourier transform (FFT), time series of the radar voltages (for every radar range) are converted to frequency domain. As the time separation between these samples is equal to the interpulse period, the maximum frequency of the spectrum is equal to 75 Hz (i.e., half the sampling frequency). This is equivalent to a maximum Doppler velocity of 225 m/s.1 The resultant self- and cross-spectrum estimates are then averaged with the results of the subsequent data blocks. Approximately 280 spectra are averaged in 5 minutes of data integration. This procedure reduces the statistical uncertainties of the estimated spectra. In Figure 7.4, we present sample plots of the self- and cross-spectrum measurements carried out at Jicamarca while operating in the 3Ba and 3Bb modes. The top panels correspond to the data obtained in the 3Ba experiment of June 19, 2008, and the bottom panels correspond to the 3Bb experiment of January 9, 2009. The data presented in these plots were obtained after five minutes of incoherent integration. In each row, the first and second plots are the self-spectra of the signals collected with one pair of antenna beams (i.e., with either the west, east, or south beams). In addition, the third and fourth plots are the magnitude and phase of the cross-spectrum between the same pair of signals. Note that the magnitude of the cross-spectrum is normalized by the square root of the product of the associated self-spectra (i.e., the plot is the magnitude of the coherence spectrum). 1The Doppler velocity vd is related to the frequency shift f by vd = λ 2 f, where λ is the radar wavelength, e.g., λ = 6 m at Jicamarca. 100 R an ge [k m] W1−Up [dB] −200 0 200 300 400 500 600 700 W2−Up [dB] −200 0 200 W1 W2* [Mag] −200 0 200 W1 W2* [Phase] −200 0 200 R an ge [k m] E1−Dn [dB] −200 0 200 300 400 500 600 700 E2−Dn [dB] −200 0 200 E1 E2* [Mag] −200 0 200 E1 E2* [Phase] −200 0 200 Velocity [m/s] R an ge [k m] S−Up [dB] −200 0 200 300 400 500 600 700 Velocity [m/s] S−Dn [dB] −200 0 200 Velocity [m/s] SU SD* [Mag] −200 0 200 Velocity [m/s] SU SD* [Phase] −200 0 200 DEWD 3Ba − Date: 19−Jun−2008 13:00:00 50 52 54 56 50 52 54 56 0 0.2 0.4 0.6 0.8 1 −150 −100 −50 0 50 100 150 49 50 51 52 53 54 55 49 50 51 52 53 54 55 0 0.2 0.4 0.6 0.8 1 −150 −100 −50 0 50 100 150 44 45 46 47 48 49 50 44 45 46 47 48 49 50 0 0.2 0.4 0.6 0.8 1 −150 −100 −50 0 50 100 150 R an ge [k m] W1−Up [dB] −200 0 200 300 400 500 600 700 W2−Up [dB] −200 0 200 W1 W2* [Mag] −200 0 200 W1 W2* [Phase] −200 0 200 R an ge [k m] E1−Dn [dB] −200 0 200 300 400 500 600 700 E2−Dn [dB] −200 0 200 E1 E2* [Mag] −200 0 200 E1 E2* [Phase] −200 0 200 Velocity [m/s] R an ge [k m] S−Up [dB] −200 0 200 300 400 500 600 700 Velocity [m/s] S−Dn [dB] −200 0 200 Velocity [m/s] SU SD* [Mag] −200 0 200 Velocity [m/s] SU SD* [Phase] −200 0 200 DEWD 3Bb − Date: 09−Jan−2009 13:20:00 50 52 54 56 50 52 54 56 0 0.2 0.4 0.6 0.8 1 −150 −100 −50 0 50 100 150 49 50 51 52 53 54 55 49 50 51 52 53 54 55 0 0.2 0.4 0.6 0.8 1 −150 −100 −50 0 50 100 150 49 50 51 52 53 54 49 50 51 52 53 54 0 0.2 0.4 0.6 0.8 1 −150 −100 −50 0 50 100 150 Figure 7.4 Sample self- and cross-spectra measured in the 3Ba experiment of June 19, 2008 (top rows), and in the 3Bb experiment of January 9, 2009 (bottom rows). The sample data were obtained after five minutes of integration. The self- and cross-spectra are plotted as functions of Doppler velocity and radar range. 101 In the plots, we can appreciate that the spectra measured with the west and east beams have narrow frequency bandwidths, as is expected for radar observations with antenna beams pointed perpendicular to B. Note that the magnitude of the associated cross- spectrum measurements also have narrow shapes. As mentioned before, the Doppler shifts of the spectra measured with the west and east beams can be used to determine vertical and zonal plasma drifts. The fact that the perpendicular-to-B spectra have very sharp peaks makes it possible to estimate with high accuracy the Doppler shifts and the plasma drifts. The signal processing technique developed for the estimation of these parameters was the subject of previous research [Kudeki et al., 1999] and is not discussed in this dissertation. We can also see in Figure 7.4 that the spectra measured with the south beams, which are pointed off-perpendicular to B, are very flat. The reason for this is that the bandwidth of the off-perpendicular radar signal is greater than the pulse-to-pulse sampling frequency. Therefore, the measured spectra are aliased. In addition, note that the signals collected with the south beams experience Faraday rotation effects. As these beams are orthogonally polarized, we can see how the polarization of the backscattered signal rotates, i.e., how the signal goes from one polarization to the other as the radar range increases. These effects will be discussed in more detail later in this section. In view of Parseval’s theorem, the frequency components of the self- and cross-spectrum measurements are summed to obtain signal power and cross-correlation profiles. In Figure 7.5, we present range vs. time plots of the power and correlation data measured in the 3Ba experiment of June 19, 2008. In the left column, the power data collected with all the radar beams are displayed in linear scale. Additionally, the cross-correlation data of the three pairs of beams are presented in the right column. The magnitudes of the correlation are displayed in linear scale, and the phases are plotted in degrees. First, note that there are marked differences between the power and correlation data of the perpendicular-to-B (west and east) beams and the data of the off-perpendicular (south) beams. As studied in the previous chapter, radiowaves propagating through the ionosphere experience changes in their polarizations due to the presence of the Earth’s 102 R an ge [k m] DEWD 3Ba − Signal Power Map (West beams) − Date: 19−Jun−2008 Beam: W1−Up 250 300 350 400 450 500 550 600 Si gn al 0 100 200 300 400 Time [Hours] R an ge [k m] Beam: W2−Up 6 8 10 12 14 16 18 20 22 24 250 300 350 400 450 500 550 Si gn al 0 100 200 300 400 R an ge [k m] DEWD 3Ba − Cross−Correlation Map (West beams) − Date: 19−Jun−2008 Pair: W1 W2* 250 300 350 400 450 500 550 600 M ag ni tu de 0 50 100 150 200 Time [Hours] R an ge [k m] Pair: W1 W2* 6 8 10 12 14 16 18 20 22 24 250 300 350 400 450 500 550 Ph as e [de g] −150 −100 −50 0 50 100 150 R an ge [k m] DEWD 3Ba − Signal Power Map (East beams) − Date: 19−Jun−2008 Beam: E1−Dn 250 300 350 400 450 500 550 600 Si gn al 0 50 100 150 200 250 300 Time [Hours] R an ge [k m] Beam: E2−Dn 6 8 10 12 14 16 18 20 22 24 250 300 350 400 450 500 550 Si gn al 0 50 100 150 200 250 300 R an ge [k m] DEWD 3Ba − Cross−Correlation Map (East beams) − Date: 19−Jun−2008 Pair: E1 E2* 250 300 350 400 450 500 550 600 M ag ni tu de 20 40 60 80 100 120 140 Time [Hours] R an ge [k m] Pair: E1 E2* 6 8 10 12 14 16 18 20 22 24 250 300 350 400 450 500 550 Ph as e [de g] −150 −100 −50 0 50 100 150 R an ge [k m] DEWD 3Ba − Signal Power Map (South beams) − Date: 19−Jun−2008 Beam: S−Up 250 300 350 400 450 500 550 600 Si gn al 0 50 100 150 Time [Hours] R an ge [k m] Beam: S−Dn 6 8 10 12 14 16 18 20 22 24 250 300 350 400 450 500 550 Si gn al 0 50 100 150 R an ge [k m] DEWD 3Ba − Cross−Correlation Map (South beams) − Date: 19−Jun−2008 Pair: SU SD* 250 300 350 400 450 500 550 600 M ag ni tu de 0 10 20 30 40 50 Time [Hours] R an ge [k m] Pair: SU SD* 6 8 10 12 14 16 18 20 22 24 250 300 350 400 450 500 550 Ph as e [de g] −150 −100 −50 0 50 100 150 Figure 7.5 Range vs. time plots of the signal power and cross-correlation data measured in the 3Ba experiment of June 19, 2008. On the left, the power data collected by each of the radar beams are displayed in linear scale. On the right, the magnitudes of the cross-correlation data are also plotted in linear scale, while the phase data are plotted in degrees. 103 magnetic field. These polarization variations are dependent on the modes of propagation excited by the radar beams. As the properties of these modes vary with magnetic aspect angle, power and correlation profiles measured with antenna beams pointing in different directions will have, to some extent, distinct shapes. When the beam is pointed (at least one degree) away from perpendicular to B, the modes excited by the beam are in the Faraday rotation regime, i.e., the propagating fields rotate. This behavior can be observed in the power plots of the south beams. As these beams are orthogonally polarized, we can see how power goes from one polarization to the other as the radar range increases. A minimum of power in one polarization corresponds to a maximum in the other polarization. Note that the transition between consecutive minimum and maximum values in a power profile corresponds to a 90◦ rotation of the polarization of the detected field. Faraday rotation effects can also be observed in the cross-correlation plots. The transition between consecutive minimum values in a profile of cross-correlation magnitude, minima that coincide with sudden changes in the correlation phase, is related to a 180◦ rotation of the polarization vector. As studied in the previous chapter, the modes of magneto-ionic propagation vary quite rapidly around the perpendicular-to-B direction. The modes range from the Cotton- Mouton regime (at exact perpendicularity) to the Faraday rotation regime (at aspect angles around and greater than a degree). As the perpendicular-to-B (west and east) beams in the 3Ba (and also in the 3Bb) experiments have widths of the order of a degree, the signals collected with these beams are the result of the combination of different modes of propagation. This complicates the description of the measurements collected with these beams. Although there are no apparent signatures of Faraday rotation in the power plots of the west and east beams, some amount of power actually “rotates.” However, note that the west beams as well as the east beams are co-polarized; therefore, the fraction of power that goes to the orthogonal polarization is not measured in these experiments. Additionally, note that these pairs of antennas are in practice radar interferometers aligned in the east-west direction. As the characteristics of the propagating fields vary in the north- 104 south direction due to magneto-ionic propagation effects, the correlation profiles measured with these interferometers are not dependent on the aspect angle variations of the radar signals. Using the correlation and power information, we have computed coherence profiles and noticed that the magnitudes and phases of the coherences are constant in height — a result that is expected, as we discuss next. The coherence measured with a radar interferometer (assuming a beam-filling scattering process) depends on the shapes of the antenna beams along the direction of alignment of the antennas. The magnitude of the coherence is determined by the widths of the radar beams, while the phase is dependent on the pointing direction of the antennas. As mentioned before, magneto-ionic propagation effects do not modify the shapes of beams in the east-west direction; therefore, as we found in our measurements, the coherence values should be equal at all radar ranges. In addition, we noted that the magnitudes of the coherences are approximately equal to the theoretical values presented in the previous section (see Figure 7.2). This result gives us confidence that the actual shapes of the radar beams are very close to the theoretical patterns. Similar plots, but for the data measured in the 3Bb experiment of January 9, 2009, are presented in Figure 7.6. In general, the features of the power and correlation measurements presented in these plots are similar to the features of the 3Ba measurements displayed in Figure 7.5. However, there are some differences that need to be pointed out. First, note that the power measured with one of the west beams is stronger than the power of the other beam. A similar difference can be observed between the signals collected with the east beams. The reason for these differences is that the antennas of each pair are of different sizes, one of them is twice as large as the other (see the description of the antenna configurations in the previous section). In addition, note that the west- and east-beam pairs are also radar interferometers, but the antennas are now aligned in the north-south direction. Therefore, we expect that magneto-ionic propagation effects will have an impact on the cross-correlation profiles measured with these interferometers. However, we can observe that the phase of the correlation is almost constant as function of height. Although altitudinal variations of the phase were expected, there are probably two reasons to explain the lack of variations. First, note that, due to solar minimum 105 R an ge [k m] DEWD 3Bb − Signal Power Map (West beams) − Date: 09−Jan−2009 Beam: W1−Up 250 300 350 400 450 500 550 600 Si gn al 0 200 400 600 800 Time [Hours] R an ge [k m] Beam: W2−Up 6 8 10 12 14 16 18 20 22 250 300 350 400 450 500 550 Si gn al 0 200 400 600 800 R an ge [k m] DEWD 3Bb − Cross−Correlation Map (West beams) − Date: 09−Jan−2009 Pair: W1 W2* 250 300 350 400 450 500 550 600 M ag ni tu de 0 50 100 150 200 250 300 Time [Hours] R an ge [k m] Pair: W1 W2* 6 8 10 12 14 16 18 20 22 250 300 350 400 450 500 550 Ph as e [de g] −150 −100 −50 0 50 100 150 R an ge [k m] DEWD 3Bb − Signal Power Map (East beams) − Date: 09−Jan−2009 Beam: E1−Dn 250 300 350 400 450 500 550 600 Si gn al 100 200 300 400 500 600 Time [Hours] R an ge [k m] Beam: E2−Dn 6 8 10 12 14 16 18 20 22 250 300 350 400 450 500 550 Si gn al 100 200 300 400 500 600 R an ge [k m] DEWD 3Bb − Cross−Correlation Map (East beams) − Date: 09−Jan−2009 Pair: E1 E2* 250 300 350 400 450 500 550 600 M ag ni tu de 0 50 100 150 200 250 Time [Hours] R an ge [k m] Pair: E1 E2* 6 8 10 12 14 16 18 20 22 250 300 350 400 450 500 550 Ph as e [de g] −150 −100 −50 0 50 100 150 R an ge [k m] DEWD 3Bb − Signal Power Map (South beams) − Date: 09−Jan−2009 Beam: S−Up 250 300 350 400 450 500 550 600 Si gn al 0 200 400 600 Time [Hours] R an ge [k m] Beam: S−Dn 6 8 10 12 14 16 18 20 22 250 300 350 400 450 500 550 Si gn al 0 200 400 600 R an ge [k m] DEWD 3Bb − Cross−Correlation Map (South beams) − Date: 09−Jan−2009 Pair: SU SD* 250 300 350 400 450 500 550 600 M ag ni tu de 0 50 100 150 200 250 Time [Hours] R an ge [k m] Pair: SU SD* 6 8 10 12 14 16 18 20 22 250 300 350 400 450 500 550 Ph as e [de g] −150 −100 −50 0 50 100 150 Figure 7.6 Same plots as described in Figure 7.5, but for the signal power and cross- correlation data measured in the 3Bb experiment of January 9, 2009. 106 conditions, the ionospheric densities in the last few years are smaller than typical values; thus, magneto-ionic effects were relatively weak during these experiments. On top of that, note that the separation between the interferometric antennas is not as large as in the 3Ba configuration (or in the DVD configuration of the previous chapter). Because the phase variations are amplified by the separation between antennas, the 3Bb north-south interferometers are not very sensitive to magneto-ionic effects. In the following section, the different features of the power and cross-correlation data presented above are modeled based on the incoherent scatter and magneto-ionic propaga- tion theories developed in the previous chapters. The model is then used to invert physical parameters of the equatorial ionosphere. 7.3 Data Model and Inversion Technique In this section, we describe the forward model of the power and cross-correlation mea- surements gathered in the 3Ba and 3Bb experiments. The data inversion technique used for the estimation of electron densities and Te/Ti profiles is also described. Additionally, some sample plots of the outcomes of the fitting results are presented. In Chapter 5, a model for the propagation of radiowaves in a magnetized ionosphere was developed. Based on this model, soft-target radar equations (5.60) and (5.62) were formulated to represent the power and the cross-correlation of the incoherent scatter sig- nals collected by the radar receivers. These equations take into account the aspect angle variations of the radar signals due to incoherent scatter and magneto-ionic propagation effects. In order to simplify these equations, it can be assumed that the bandwidth of the incoherent scatter RCS spectrum σv(k, ω) is narrower than the frequency response of the matched-filter receivers used for the detection of the antenna signals (a valid approxima- tion in the case of short-pulse radar applications). After simplifying equations (5.60) and (5.62), the power of the signal collected by the i-th antenna receiver can be expressed as Pi = Ki ˆ dr |χ(t− 2rc )|2 4pir2 ˆ dΩWi(r)σv(k) +Ni, (7.1) 107 while the complex cross-correlation of the signals collected by the i-th and j-th antenna receivers can be modeled as Ci,j = Ki,j ˆ dr |χ(t− 2rc )|2 4pir2 ˆ dΩWi,j(r)σv(k) +Ni,j . (7.2) Above, Ki and Ki,j are calibration constants that account for the gains and losses along the paths of the radar signals. These constants can take the forms Ki = |gi|2 and Ki,j = gig∗j , (7.3) where gi and gj are complex proportionality factors for the voltages detected by the i- and j-th radar receivers. In addition, χ(t) is the normalized auto-correlation of the pulse waveform f(t) that can be expressed as χ(t) = 1 T f∗(−t) ∗ f(t). (7.4) In equations (7.1) and (7.2), the volumetric incoherent scatter RCS σv(k) is weighted by the effective two-way radiation patterns Wi(r) and Wi,j(r), which are smoothly varying functions of the radar range that account for the magneto-ionic propagation effects on the radar signals. In addition, Ni and Ni,j denote the power and cross-correlation of the noise signals detected by the antenna receivers. As mentioned before, in the 3Ba and 3Bb experiments, electromagnetic waves are transmitted simultaneously with the antennas pointed in the west and east directions. The west antennas are up-polarized, while the east ones are down-polarized. Despite the orthogonality of the polarizations of these antennas, some amount of power may go from the west beam to the east beam, or vice versa, due to magneto-ionic propagation effects. In order to take into account this possible “cross-talk,” the expressions for the effective radiation patterns in (7.1) and (7.2) become somewhat complicated. Denoting the pair of west beams by w1 and w2 and the pair of east beams by e1 and e2, the effective two-way 108 radiation patterns for the power model (7.1) can take the following forms, Ww1(r) = |Rw1(rˆ)|2 k2 ( |Tw(rˆ)|2 Γuu(r) + |Te(rˆ)|2 Γdu(r) + 2Re { Tw(rˆ) T ∗e (rˆ) Γu,du (r) }) , (7.5) Ww2(r) = |Rw2(rˆ)|2 k2 ( |Tw(rˆ)|2 Γuu(r) + |Te(rˆ)|2 Γdu(r) + 2Re { Tw(rˆ) T ∗e (rˆ) Γu,du (r) }) , (7.6) We1(r) = |Re1(rˆ)|2 k2 ( |Te(rˆ)|2 Γdd(r) + |Tw(rˆ)|2 Γud(r) + 2Re { Te(rˆ) T ∗w (rˆ) Γd,ud (r) }) , (7.7) We2(r) = |Re2(rˆ)|2 k2 ( |Te(rˆ)|2 Γdd(r) + |Tw(rˆ)|2 Γud(r) + 2Re { Te(rˆ) T ∗w (rˆ) Γd,ud (r) }) , (7.8) while the patterns for the cross-correlation model (7.2) can be expressed as Ww1,w2(r) = Rw1(rˆ)R∗w2(rˆ) k2 ( |Tw(rˆ)|2 Γuu(r) + |Te(rˆ)|2 Γdu(r) + 2Re { Tw(rˆ) T ∗e (rˆ) Γu,du (r) }) , (7.9) We1,e2(r) = Re1(rˆ)R∗e2(rˆ) k2 ( |Te(rˆ)|2 Γdd(r) + |Tw(rˆ)|2 Γud(r) + 2Re { Te(rˆ) T ∗w (rˆ) Γd,ud (r) }) , (7.10) where k = 2pi/λB = 4pi/λ is the Bragg wavenumber. Above, Tw(rˆ) and Te(rˆ) are the antenna-array factors of the west and east transmitting antennas. In addition, Rw1(rˆ), Rw2(rˆ), Re1(rˆ), and Re2(rˆ) are the antenna-array factors of the receiving antennas. These factors are properly normalized such that their squared magnitudes are antenna gain pat- terns. Note that the products of the forms |Tw(rˆ)|2|Rw(rˆ)|2 and |Te(rˆ)|2|Re(rˆ)|2 correspond to the two-way radiation patterns displayed in Figure 7.2. In the previous equations, the coefficients multiplying the antenna-array factors can be expressed in the following form, Γk,li,j (r) = ( pˆTi Π¯(r) pˆk )( pˆTj Π¯(r) pˆl )∗ , (7.11) and represent polarization coefficients that account for the magneto-ionic propagation effects on the polarizations of the returned fields. In this expression, Π¯(r) denotes the two-way propagator matrix defined in Section 5.2. In addition, pˆi, pˆj , pˆk, and pˆl are 109 polarization unit vectors. The indices i and j denote the receiving antennas, while k and l denote the transmitting antennas. In the equations for the effective radiation patterns, these indices are either u for the up-polarization, or d for the down-polarization of the Jicamarca antenna. Expressions for the polarization unit vectors pˆu and pˆd are given in the previous chapter. Furthermore, in order to simplify the notation, we consider the following equivalences Γki (r) ≡ Γk,ki,i (r), Γk,li (r) ≡ Γk,li,i (r), and Γki,j(r) ≡ Γk,ki,j (r). Similarly, denoting the pair of south beams by s1 and s2, the effective two-way radiation patterns for the corresponding power models take the following forms, Ws1(r) = |Rs1(rˆ)|2 k2 |Ts(rˆ)|2 Γuu(r), (7.12) Ws2(r) = |Rs2(rˆ)|2 k2 |Ts(rˆ)|2 Γud(r), (7.13) while the effective pattern for the cross-correlation model can be expressed as Ws1,s2(r) = Rs1(rˆ)R∗s2(rˆ) k2 |Ts(rˆ)|2 Γuu,d(r). (7.14) Above, Ts(rˆ) is the antenna-array factor of the transmitting antenna that is up-polarized, while Rs1(rˆ) and Rs2(rˆ) are the antenna-array factors of the receiving antennas that are up- and down-polarized, respectively. In the set of equations that model the power and correlation measurements, the vol- umetric incoherent scatter RCS σv(kˆ) and the two-way propagator matrix Π¯(r) are the only functions that depend on ionospheric parameters. The RCS σv(kˆ) depends mainly on electron density Ne, temperature ratio Te/Ti, and magnetic aspect angle α, while the matrix Π¯(r) only depends on the electron density and the magnetic field. Using the mag- netic field values provided by the IGRF model [e.g., Olsen et al., 2000], we have applied the power and cross-correlation models presented above to invert Ne and Te/Ti profiles from the data collected in the 3Ba and 3Bb experiments. For this purpose, we consider the ionosphere to be horizontally stratified. The electron 110 density profile Ne(z) is represented as a piecewise linear function Ne(z) = L∑ l=1 N le4 ( z − zl ∆z ) (7.15) for z ≥ z1, where 4(t) is the standard triangular function 4(t) =  1− |t|, |t| < 1 0, otherwise. (7.16) The density values at the nodes zl are N le = Ne(zl). The nodes are spaced every ∆z = 15 km (which is the length of the baud of the Barker code used in transmission). We have considered z1 = 210 km and zL = 600 km, thus L = 27. The density profile below z1 is also modeled by linear functions characterized by a single parameter N0e at a node z0 = 150 km. It is assumed that below z = 85 km the electron density vanishes. In addition, the temperature ratio profile Te(z)/Ti(z) is modeled as a layer with a Gaussian shape Te(z)/Ti(z) = Tm exp ( −(z − zT ) 2 2w2T ) + 1, (7.17) where the constraint Te/Ti ≥ 1 is implemented by default. This profile is characterized by the following parameters: Tm + 1 is the peak of the temperature ratio, zT is the altitude where Te/Ti peaks, and wT is the width of the Gaussian profile. Modeling the temperature ratio profile as a single layer can be considered a valid approximation for the range of altitudes probed in these experiments. Typical measurements of equatorial ionospheric temperatures (with other radar techniques at Jicamarca) show that, at daytime hours, there is a layer in the bottom-side of the F-region where Te/Ti ≥ 1. In the numerical implementation of the data model presented above, we have computed samples of the power and cross-correlation data every ∆r = 15 km at radar ranges rn = n∆r where n is an integer. For this purpose, equations (7.1) and (7.2) have been recast 111 as discrete convolutions taking the following forms: Pi[n] = |gi|2 |χ[n]|2 ∗ ( 1 4pir2n ˆ dΩWi(rn)σv(k) ) +Ni (7.18) and Ci,j [n] = gi g ∗ j |χ[n]|2 ∗ ( 1 4pir2n ˆ dΩWi,j(rn)σv(k) ) +Ni,j , (7.19) where χ[n] = 1 M c[−n] ∗ c[n] (7.20) is the discrete auto-correlation of the code c[n] of length M . In our 3Ba and 3Bb experi- ments, the radar pulses were modulated by a 3-baud Barker sequence + +− with a baud length equal to 15 km; thus, c[n] = [1, 1,−1] and M = 3. In this case, it can be verified that χ[n] is given by χ[n] = −1 3 δ[n+ 2] + δ[n]− 1 3 δ[n− 2]. (7.21) In order to numerically evaluate the solid angle integrals in (7.18) and (7.19), we reformu- late the integrals in terms of θx and θy, which are direction cosines of a coordinate system in which the x- and y-axes are aligned with the east and north directions (see Figures 7.2 and 7.3). Note that the differential solid angle dΩ becomes dθxdθy/ √ 1− θ2x − θ2y . Using this coordinate system, we can take advantage of the magnetic aspect angle dependence of the functions Wi(rn), Wi,j(rn), and σv(k). As we know, these functions vary quite rapidly in the north-south direction, but very slowly in the east-west direction. Therefore, values of these functions are computed in a grid that is finely sampled in the θy-direction with a separation between samples of 0.002 dc (∼ 0.11◦), and coarsely sampled in the θx-direction with a separation of 0.113 dc (∼ 6.5◦). A range of approximately ±20◦ with respect to the antenna on-axis direction is covered in both θx- and θy-directions. The grid is composed of 349 nodes in θy and 7 nodes in θx; thus, in order to calculate the integrals, a total of 2443 samples is computed for every function. A finite-element-like integration method based on rectangular elements is used to evaluate the integrals. Note that the same operation has 112 to be repeated for every range gate. A trust-region approach for minimizing non-linear least-squares problems subject to simple inequality conditions [e.g., Coleman and Li , 1996] was used for the inversion of the power and cross-correlation data. In comparison with the standard Levenberg-Marquardt algorithm, the trust-region method is more robust for large-scale problems. Additionally, the algorithm allows us to impose lower and upper bounds on the solution of the problem. The following list of parameters are fitted using this minimization algorithm. • Electron densities N le at the nodes zl for 0 ≤ l ≤ L. • Temperature ratio parameters Tm, zT , and wT . • Calibration factors gi for the six antenna channels. • Noise power levels Ni for the six antenna channels. • Noise cross-correlation values Ni,j for the three pairs of antenna channels. Note that gi and Ni,j are complex numbers. The parameters listed above are estimated by minimizing the following cost function E2 = E2w1,w2 + E2e1,e2 + E2s1,s2, (7.22) where E2i,j = ∑ n δmTi,j [n] M¯ −1 i,j [n] δmi,j [n] (7.23) is the error cost function of the data collected with the i-th and j-th antenna channels. In this expression, δmi,j [n] is the vector of the differences between the measured and modeled data at the radar range rn that is defined as δmi,j [n] =  δPi[n] δPj [n] δRi,j [n] δQi,j [n]  =  pi[n]− Pi[n] pj [n]− Pj [n] ri,j [n]−Ri,j [n] qi,j [n]−Qi,j [n]  , (7.24) 113 where pi[n] and pj [n] denote the power measurements, and ri,j [n] and qi,j [n] correspond to the real and imaginary parts of the cross-correlation data. In addition, Pi[n], Pj [n], Ri,j [n] = Re{Ci,j [n]}, and Qi,j [n] = Im{Ci,j [n]} are the modeled power (7.18) and cross- correlation (7.19). In the cost function (7.23), M¯i,j [n] is the covariance matrix of the vector of data measurement errors δmi,j [n], i.e., M¯i,j [n] = 〈 δmTi,j [n] δmi,j [n] 〉 . Assuming that δmi,j [n] is a zero-mean Gaussian random vector, M¯i,j [n] can be written in the following form [e.g., Feng et al., 2004] M¯i,j = 1 I  P 2i R 2 i,j +Q 2 i,j PiRi,j PiQi,j R2i,j +Q 2 i,j P 2 j PjRi,j PjQi,j PiRi,j PjRi,j 1 2 ( R2i,j −Q2i,j + PiPj ) Ri,jQi,j PiQi,j PjQi,j Ri,jQi,j 1 2 ( Q2i,j −R2i,j + PiPj )  , (7.25) where I is the number of samples used to compute the power and cross-correlation mea- surements. In this expression, indices n were omitted in order to simplify the notation. The assumption of a Gaussian vector δmi,j [n] can be justified because every element of δmi,j [n] results from the average of a large number of nearly uncorrelated samples. Note that the radar data is collected every 5 km; however, in the cost function (7.23), we have only considered samples that are spaced every 15 km. The form of the cost function (7.23) is valid if the samples at every range gate are not correlated. This is our case, as samples spaced at least a distance equal to the resolution of the pulse shape (which is 15 km) can be considered to be uncorrelated. However, in order to include the oversampled data, a more complicated cost function has to be considered because the data corresponding to neighboring gates spaced every 5 km are correlated. The implementation of the inversion algorithm including the oversampled data will be the subject of future research. The sum in (7.23) is taken over the set of data between 210 km and 600 km; thus 14 ≤ n ≤ 40, which implies that 27 range gates were fitted. In order to estimate the parameters listed above, we apply the following two-step strategy. 114 First step • Initial values for the electron density, noise, and calibration parameters are deter- mined. In this step, it is assumed that Te/Ti = 1 at all heights. • The following constraints are specified: electron densities N le and noise levels Ni must be positive, and gw1, ge1, and gs1 must be real and positive. The rest of the parameters are not constrained. • The trust-region method is used to minimize the cost function (7.22) subject to the constraints specified above. In this first step, we only fitted for N le, Ni, Ni,j , and gi. Second step • The estimated parameters obtained in the first step are used as an initial guess. In addition, we define initial values for Tm, zT , and wT . • The same constraints specified in the first step are considered in this one. Addition- ally, the Te/Ti parameters are restricted to be positive. • The trust-region method is used again to minimize the cost function (7.22) but including in the minimization the Te/Ti parameters. In this second step, we solve for all the parameters listed above, i.e., N le, Tm, zT , wT , Ni, Ni,j , and gi. We apply the two-step strategy because we notice that the inversion of N le is more robust than the inversion of the Te/Ti parameters. When fitting directly for the densities and the Te/Ti parameters all together in a single step, the solutions for Te/Ti tend to be erratic, which is probably related to the fact that the power and cross-correlation data is somewhat less dependent on Te/Ti than on the densities. Estimating the densities in a first step and subsequently the densities plus the Te/Ti parameters provides better results. Note that we have not considered any explicit regularization term in the cost function (7.22). On the other hand, modeling the Te/Ti profile as a Gaussian function (characterized with a few parameters) regularizes, in some sense, the problem because a continuity condition on 115 the Te/Ti profile is imposed. The mathematical relationship between our approach and standard regularization procedures will be the subject of future studies. In Figure 7.7, sample plots of the results obtained from the inversion of the 3Ba and 3Bb data are displayed. In these plots, the estimated electron density and Te/Ti profiles are presented in the left column. In the other columns, the power and cross-correlation measurements (dotted points) are compared with the modeled profiles resultant from the inversion (continuous lines). We can appreciate that the comparison is good, because measured and modeled profiles match closely; however, there are some discrepancies (par- ticularly in the cross-correlation data of the south-beams of the 3Ba experiment) that might be related to a poor characterization of the bottom part of the density profile. Fur- ther analysis and discussion regarding the characteristics of the results are presented in the following section. 7.4 Results and Comparisons In Figure 7.8, we present plots of the electron density and Te/Ti estimates obtained from the inversion of incoherent scatter power and cross-correlation data collected in the 3Ba experiment of June 19, 2008 (left column), and also in the 3Bb experiment of Jan- uary 9, 2009 (right column). To obtain these results, the data model and fitting technique described in the previous sections were applied. Both the Ne and Te/Ti profiles are plotted in linear scale. First, note that the estimated densities and temperature ratios are within the ranges of values expected for an equatorial ionosphere. Electron densities of the order of 1011 m−3 are typical of solar-minimum conditions. A layer of Te/Ti > 1 located at the bottom of the F-region is characteristic of the daytime ionosphere. The Te/Ti layer has peak values greater than two and relaxes toward thermal equilibrium (Te = Ti) as time reaches the night hours. However, we can also observe that the Te/Ti values fluc- tuate from one profile to the next one. The origin of these fluctuations is the subject of current analysis. Apparently, the magnitudes of the fluctuations are of the order of the statistical uncertainties expected for these estimates. Note that the radar data collected in the 3Ba and 3Bb experiments can be considered (to some extent) to be less dependent 116 0 1 2 3 4 x 1011 250 300 350 400 450 500 550 N e [m−3] H ei gh t [k m] 1 1.5 2 250 300 350 400 450 500 550 T e /Ti H ei gh t [k m] 400 600 800 250 300 350 400 450 500 550 Power (East−West) R an ge [k m] W1−Up W2−Up E1−Dn E2−Dn 150 200 250 300 250 300 350 400 450 500 550 Power (South) R an ge [k m] S−Up S−Dn 0 100 200 250 300 350 400 450 500 550 CCF (Real) R an ge [k m] W1 W2* E1 E2* SU SD* −5 0 5 10 15 250 300 350 400 450 500 550 CCF (Imag) R an ge [k m] DEWD 3Ba − Date: 19−Jun−2008 13:00:00 0 2 4 6 x 1011 250 300 350 400 450 500 550 600 N e [m−3] H ei gh t [k m] 1 1.5 2 250 300 350 400 450 500 550 600 T e /Ti H ei gh t [k m] 400 600 800 1000 250 300 350 400 450 500 550 600 Power (East−West) R an ge [k m] W1−Up W2−Up E1−Dn E2−Dn 400 600 800 250 300 350 400 450 500 550 600 Power (South) R an ge [k m] S−Up S−Dn −200 −100 0 100 200 250 300 350 400 450 500 550 600 CCF (Real) R an ge [k m] W1 W2* E1 E2* SU SD* 0 50 100 150 250 300 350 400 450 500 550 600 CCF (Imag) R an ge [k m] DEWD 3Bb − Date: 09−Jan−2009 13:20:00 Figure 7.7 Sample plots of the results obtained from the inversion of the 3Ba and 3Bb data (top and bottom panels, respectively). The estimated Ne and Te/Ti profiles are plotted in the left column. In the other columns, the power and cross-correlation measurements (dotted points) are compared with the modeled profiles (continuous lines). 117 Ti m e [H ou rs] Height [km] D EW D 3 Ba − E le ct ro n D en si ty (N e) − D ate : 1 9− Ju n− 20 08 6 8 10 12 14 16 18 20 22 24 25 0 30 0 35 0 40 0 45 0 50 0 55 0 60 0 Ne [m −3 ] 00. 5 11. 5 22. 5 33. 5 44. 5 55. 5 x 10 11 Ti m e [H ou rs] Height [km] D EW D 3 Bb − E le ct ro n De ns ity (N e) − D ate : 0 9− Ja n− 20 09 6 8 10 12 14 16 18 20 22 25 0 30 0 35 0 40 0 45 0 50 0 55 0 60 0 Ne [m −3 ] 01234567x 10 11 Ti m e [H ou rs] Height [km] D EW D 3 Ba − T em pe ra tu re ra tio (T e/T i) − D ate : 1 9− Ju n− 20 08 6 8 10 12 14 16 18 20 22 24 25 0 30 0 35 0 40 0 45 0 50 0 55 0 60 0 Te/Ti 11. 1 1. 2 1. 3 1. 4 1. 5 1. 6 1. 7 1. 8 1. 9 2 Ti m e [H ou rs] Height [km] D EW D 3 Bb − T em pe ra tu re ra tio (T e/T i) − D ate : 0 9− Ja n− 20 09 6 8 10 12 14 16 18 20 22 25 0 30 0 35 0 40 0 45 0 50 0 55 0 60 0 Te/Ti 11. 1 1. 2 1. 3 1. 4 1. 5 1. 6 1. 7 1. 8 1. 9 2 F ig ur e 7. 8 R an ge vs . ti m e pl ot s of th e el ec tr on de ns it y an d T e /T i es ti m at es ob ta in ed fr om th e in ve rs io n of in co he re nt sc at te r si gn al po w er an d cr os s- co rr el at io n da ta co lle ct ed in th e 3B a ex pe ri m en t of Ju ne 19 , 20 08 (l ef t co lu m n) , an d al so in th e 3B b ex pe ri m en t of Ja nu ar y 9, 20 09 (r ig ht co lu m n) . 118 on Te/Ti than on Ne. As we know, if Te/Ti > 1, the incoherent scatter RCS increases as the magnetic aspect angle α decreases and reaches its maximum value at α = 0◦ (i.e., at the angle of perpendicular-to-B propagation). However, the power increments mea- sured with perpendicular-to-B beams are somewhat smaller than the theoretical values for α = 0◦ because signal contributions with smaller RCS are averaged by the antenna beams. This reduces the sensitivity of the radar data to changes in Te/Ti. Additionally, note that the fluctuations of the Te/Ti values are larger in the 3Ba case than in the 3Bb case. This might be related to the fact that the 3Ba two-way radiation patterns are broader in the north-south direction than the 3Bb patterns, because a wider range of signal con- tributions with smaller RCS is collected by the 3Ba beams (see antenna description in Section 7.1). In addition, note that the incoherent scatter RCS is directly proportional to Ne (see Section 4.2). Moreover, the altitudinal variations of the measured signals due to magneto-ionic propagation effects are dependent on the electron density values but not on Te/Ti. Although it is likely that the Te/Ti fluctuations are due to the statistical errors of the estimated parameters, further studies need to be carried out in order to fully identify the origin of these fluctuations. In the density plots, we can observe patches of enhanced densities resembling “bead” shapes. Whether these beads were artifacts of the inversion procedure needed to be an- alyzed. First, we examined whether these events are correlated with fluctuations on the transmitted power. After analyzing records of the line voltage applied to the transmitters used in these experiments, we did not find any correlation between the fluctuations of the supplied voltage and the patches of enhanced electron density. We also noted that the voltage fluctuations follow closely the variations of the radar calibration parameters pre- sented in Figure 7.9. As there are no systematic variations that may be the cause for the density increments, we consider the origin of the density beads to be geophysical. Whether these variations are due to solar activity or due to perturbations of the neutral atmosphere will be the subjects of future investigations. It is also important to describe the characteristics of the radar calibration constants resulting from the inversion of the radar data. Note in Figure 7.9 that the calibration 119 −30 −25 −20 M ag ni tu de [d B] DEWD 3Ba − Calibration parameters − Date: 19−Jun−2008 W1−Up W2−Up E1−Dn E2−Dn S−Up S−Dn 6 8 10 12 14 16 18 20 22 24 −10 −5 0 5 Time [Hours] Ph as e [de g] W1−Up W2−Up E1−Dn E2−Dn S−Up S−Dn −24 −23 −22 −21 −20 −19 M ag ni tu de [d B] DEWD 3Bb − Calibration parameters − Date: 09−Jan−2009 W1−Up W2−Up E1−Dn E2−Dn S−Up S−Dn 6 8 10 12 14 16 18 20 22 −30 −20 −10 0 Time [Hours] Ph as e [de g] W1−Up W2−Up E1−Dn E2−Dn S−Up S−Dn Figure 7.9 Magnitudes and phases of the radar calibration constants obtained as part of the inversion of the 3Ba and 3Bb radar data. constants of the west and east antenna channels exhibit approximately the same variations in time. This is expected if we assume that the only system parameter that may fluctuate during the operation of the radar is the amount of transmitted power. However, we can observe that the features of the calibration parameters of the south channels do not follow necessarily the time variations of the calibration parameters of the other channels. The fact that two different transmitters were used in these experiments (one for the west and east beams and the other for the south beams) might be the reason for the distinct behaviors of the calibration parameters. However, we can also see that the behaviors of the calibration constants of the south beams are somewhat different from each other (particularly in the case of the 3Ba data). These differences become more evident at time intervals in which the electron temperature is greater than the ion temperature (i.e., Te/Ti > 1). The cause for these unexpected differences might be related to a poor description of the electron density profile at altitudes below 200 km (a segment of the profile that is characterized by a single parameter, as explained in Section 7.3). The shape of the lower part of the electron density profile determines the amount of Faraday rotation that the south-beam fields experience up to the first altitude of the data samples used in our inversions. If a simplified description of the bottom part of the density profile is causing an inaccurate estimation of these parameters, the use of ionosonde data to help on the inversion of the lower part of the density profile might be needed to properly model the radar data and improve the quality of the inversions. Further testing and analysis is needed in order to identify the origin of the differences between the calibration parameters. 120 6 8 10 12 14 16 18 20 22 24 0 2 4 6 8 10 12 14 16 18 Time [Hours] ε2 / N m e a su re m e n ts DEWD 3Ba − Goodness of Fit − Date: 19−Jun−2008 Te/Ti = 1 (Fixed) Te/Ti ≥ 1 6 8 10 12 14 16 18 20 22 0 2 4 6 8 10 12 Time [Hours] ε2 / N m e a su re m e n ts DEWD 3Bb − Goodness of Fit − Date: 09−Jan−2009 Te/Ti = 1 (Fixed) Te/Ti ≥ 1 Figure 7.10 Goodness-of-fit values resulting from the inversion of the 3Ba and 3Bb radar data. The blue curves correspond to the values obtained by fitting only for the electron density Ne and assuming that Te/Ti = 1 at all altitudes. The green curves correspond to the values obtained by fitting for both parameters, Ne and Te/Ti. In Figure 7.10, goodness-of-fit values obtained from the inversion of the 3Ba and 3Bb radar data are presented and can be used to analyze the quality of the inversion results. The goodness-of-fit is defined as E2/Nmeasurements. In the plots, the blue curves correspond to the values obtained by fitting only for the electron density Ne and assuming that Te/Ti = 1 at all altitudes. The green curves correspond to the values obtained by fitting for both parameters, Ne and Te/Ti. The goodness-of-fit values are around five in the case of the 3Ba results, while they are around three in the case of the 3Bb results. As expected, note that the goodness-of-fit values are smaller when both ionospheric parameters, Ne and Te/Ti, are fitted, implying that Te/Ti correction is needed. However, the goodness-of-fit values obtained are, in some sense, somewhat large, because a “good” fitting result corresponds to a goodness-of-fit value close to one. This indicates that the model has not totally captured the features of the data; therefore, there is still space for improvement. In addition, in order to validate our estimated densities, we have compared our results with ionosonde measurements conducted simultaneously at Jicamarca. In Figure 7.11, we are comparing the peak values of the plasma frequencies (foF2 = 12pi √ Nee2/meo) measured with the 3Ba and 3Bb radar techniques (blue points) and with the ionosonde measurements (red points). In addition, the altitudes corresponding to the peak values of the plasma frequencies are compared in Figure 7.12. Although the ionosonde data is noisy in the 3Ba case, we can observe that the agreement is very good, particularly in the 121 6 8 10 12 14 16 18 20 22 24 0 1 2 3 4 5 6 7 8 Time [Hours] fo F2 [M Hz ] DEWD 3Ba − F2 plasma frequency − Date: 19−Jun−2008 ISR Digisonde 6 8 10 12 14 16 18 20 22 0 1 2 3 4 5 6 7 8 9 10 11 Time [Hours] fo F2 [M Hz ] DEWD 3Bb − F2 plasma frequency − Date: 09−Jan−2009 ISR Digisonde Figure 7.11 Comparison of the peak values of plasma frequencies (foF2) measured with the 3Ba and 3Bb radar techniques (blue lines) and with an ionosonde system also located at Jicamarca (red dots). 6 8 10 12 14 16 18 20 22 24 0 50 100 150 200 250 300 350 400 Time [Hours] H ei gh t [k m] DEWD 3Ba − Height of foF2 − Date: 19−Jun−2008 ISR Digisonde 6 8 10 12 14 16 18 20 22 0 50 100 150 200 250 300 350 400 Time [Hours] H ei gh t [k m] DEWD 3Bb − Height of foF2 − Date: 09−Jan−2009 ISR Digisonde Figure 7.12 Comparison of the altitudes corresponding to the peak values of plasma frequency presented in Figure 7.11. 3Bb case. Note that the oscillations on the plasma frequency curves are correlated with the patches of enhanced densities described before. The fact that the same behavior can be observed in both the radar and the ionosonde measurements reinforces our conclusion that the origin of these patches is geophysical. 122 CHAPTER 8 CONCLUSIONS AND FUTURE WORK In this dissertation, we have studied in detail Coulomb collisions and magneto-ionic propagation effects on incoherent scatter radar measurements carried out at Jicamarca. Our studies show that both processes modify the shapes of the ISR spectra measured with antenna beams pointed perpendicular to B at 50 MHz; thus, their effects should not be neglected in the modeling of radar measurements. To perform this analysis, we developed a new incoherent scatter spectrum model that takes into account both types of effects. The model is valid for all magnetic aspect angles, including the direction perpendicular to the ambient magnetic field. In the first stage of our work, we focused on modeling the effects of Coulomb collisions. The procedure to model the IS spectrum is based on the simulation of the trajectories of charged particles moving in a magnetized plasma with suppressed collective interac- tions. The effects of Coulomb collisions are considered using the Fokker-Planck model of Rosenbluth et al. [1957]. The results of our studies show that, in contrast with the case of the ions, the statistics of the electron displacements cannot be fully approximated by simplified Coulomb collision models. Thus, in order to apply our collisional spectral model in routine computations of the IS spectrum, a numerical library of single-electron ACF’s had to be built. The library, however, was only developed for the case of oxygen plasmas. This limits the applicability of our model to the region around and below the peak of the F-region. In order to study the effects of Coulomb collisions on radar observations of the topside ionosphere, we have to extend our model to the case of hydrogen and helium plas- mas (because H+ and He+ are the dominant ion species at the top of the F-region). This 123 future work will require the simulation of a large number of particle trajectories for differ- ent plasma settings, a very computationally demanding task that will require the use of a cluster of computers, such as the Turing Cluster maintained by the Computational Science and Engineering (CSE) Program at the University of Illinois. Basic questions regarding the nature of Coulomb collisions will be addressed in these future studies. For instance, whether the Coulomb collision process in ionospheric H+ and He+ plasmas is Gaussian or not, it is an important question to answer. In the case of Gaussian particle trajec- tories, a Brownian-motion model characterized by constant collision coefficients could be used to simplify the description of the process. The collisional model resultant from this investigation will provide the Jicamarca Radio Observatory with unique capabilities. In the second stage of our studies, we have modeled the beam-weighted incoherent scatter radar spectrum taking into account magneto-ionic propagation effects. The model developed in Chapter 5 is effectively an extension of the Appleton-Hartree solution for a homogeneous magneto-plasma to the case of an inhomogeneous ionosphere. Simulation studies based on our model show that magneto-ionic propagation effectively modifies the shapes of the antenna beams used in our radar experiments. The model was developed based on the assumption that refraction effects can be neglected. Although refraction is weak at 50 MHz, given the long distances traveled by the radar signals, some features of our measurements might, to some extent, be the result of refraction effects, particularly in the case of radar observations above the F-region peak. Note that refraction may cause the bending of the propagating fields; therefore, the angular shifts of the antenna beams that were observed in our radar experiments might be caused not only by magneto-ionic polarization effects, but also by refraction effects. In order to analyze these possible effects, we are planning to develop a more sophisticated propagation model based on the WKB (Wentzel, Kramers, and Brillouin) approximation of Budden [1961] for the propagation of radiowaves in an anisotropic inhomogeneous ionosphere. In order to simplify the formulation of our spectrum model, we considered that the shape of the radar ambiguity function in frequency domain is much wider than the inco- herent scatter spectrum. Thus, the frequency dependence of the ambiguity function was 124 neglected by replacing this function with the ACF of the pulse waveform (see Chapter 5). In the case of perpendicular-to-B radar measurements, this assumption is well justified because the IS spectrum is very narrow; however, in the case of off-perpendicular radar measurements, the IS spectrum is wider, and, therefore, its shape might be filtered by the radar ambiguity function. Some preliminary calculations have shown that the power measured with off-perpendicular radar beams around the F-region peak may be attenu- ated by five percent due to ambiguity-function filtering effects. This attenuation, however, depends on the shape and length of the radar pulse, as well as on the width of the IS spectrum, which (at these aspect angles) is mainly controlled by the ionospheric temper- atures. Further studies need to be conducted in order to quantify more carefully these attenuation effects and verify if they will have an impact on the estimation of ionospheric plasma parameters. In Chapter 4, we studied the dependence of the collisional spectra on electron and ion temperature variations and noticed that the bandwidth of the perpendicular-to-B spec- trum is inversely proportional to Te/Ti; however, at larger aspect angles, the bandwidth dependence is different, because the spectral width increases with increasing Te/Ti. As the IS spectrum measured by the radar results from the combination of signal contribu- tions arriving from different directions, the dependence of the measured ISR spectrum might be more complicated, particularly because typical spectrum measurements with perpendicular-to-B beams are aliased. The analysis of the sensitivity of the beam-weighted IS spectrum model to electron and ion temperatures will be the subject of future studies. In Chapter 7, we applied our incoherent scatter data model to the estimation of electron densities and Te/Ti profiles from radar measurements carried out in multi-beam radar ex- periments at Jicamarca. In these experiments, different sections of the Jicamarca antenna array are phased to point west, east, or south of the local zenith. While the west and east beams are aimed perpendicular to the magnetic field, the south beams are pointed in a direction ∼ 4◦ to the south of the loci of perpendicularity. In addition to electron density and Te/Ti estimates, radar calibration parameters were also obtained from the inversion procedure. The quality of our inversion results was analyzed based on the goodness- 125 of-fit criterion. Moreover, comparisons between our density estimates and independent ionosonde measurements carried out simultaneously at Jicamarca showed excellent agree- ment. Although the quality of our inversion results can be considered good, a more careful analysis of the statistical errors of our estimates is still needed. There are certain features in our estimated results that require further investigation. For instance, in the inverted density maps, we can observe patches of enhanced electron density that resemble “bead” shapes, whose origin might be due to solar activity or perturbations of the neutral atmo- sphere. In addition, we can also see that the estimated Te/Ti values fluctuate from profile to profile. Although it is likely that the Te/Ti fluctuations are within the range of the statistical errors for these parameters, further analysis is needed. The ultimate application of our incoherent scatter model is the simultaneous estimation of drifts, densities, and temperatures of the equatorial ionosphere (an objective not yet achieved by the radar techniques available at Jicamarca). For this purpose, a composite experiment that interleaves multi-pulse and pulse-to-pulse radar measurements has been recently proposed for its application at Jicamarca. This experiment will be carried out with antenna beams pointed perpendicular to B. As mentioned before, the pulse-to- pulse incoherent scatter spectrum measured at Jicamarca is composed of narrow and wide bandwidth spectral components. The plasma drifts are determined from the Doppler shift of the narrow spectrum component. However, the wide spectrum component is typically aliased in pulse-to-pulse radar observations. As the bandwidth of the wide component is directly proportional to plasma temperatures, we will measure this spectrum component using a multi-pulse radar technique in an attempt to invert temperature profiles from these measurements. This will be the fulfillment of objectives set forth more than a decade ago, when spectral Te estimations were first tried [Bhattacharyya, 1998] and the inadequacy of the incoherent scatter theory close to perpendicular-to-B was first realized. 126 APPENDIX A PARTICLE DYNAMICS DESCRIPTION OF “BGK COLLISIONS” AS A POISSON PROCESS The Gordeyev integral for plasma particles colliding with neutrals is obtained using a particle dynamics formalism in which the collisions are modeled as a discrete Poisson process. The result leads to an electron density fluctuation spectrum model for partially ionized plasmas that is identical with the spectral model obtained from BGK plasma kinetic equations. This isomorphism between the Poisson process and the BGK operator is analogous to a similar relation between the Brownian-motion process and the Fokker- Planck operator with constant coefficients. We take advantage of this analogy to derive a collisional IS spectrum model that takes into account collisions with both neutrals and charged species. A.1 Introduction The purpose of the following discussion is to show that the incoherent scatter spec- trum model derived from the Boltzmann kinetic equation with the BGK collision opera- tor [Dougherty and Farley , 1963] can also be obtained by using a particle dynamics ap- proach where particle collisions are modeled as a Poisson process. This is analogous to the case of the Fokker-Planck operator (with constant coefficients) and the Brownian-motion (Ornstein-Uhlenbeck) collision process (described by the Langevin equation), which lead to identical spectral models when utilized in the Boltzmann equation and in the particle dynamics formalism, respectively [e.g., Chandrasekhar , 1943; Gillespie, 1996b; Holod et al., 2005]. 127 In general, incoherent scatter spectral theories pertinent to various types of ionospheric plasmas (magnetized, collisional, etc.) in thermal equilibrium can be expressed in terms of appropriately derived Gordeyev integrals utilized within the general framework described in Chapter 2 [e.g., Kudeki and Milla, 2006, 2010]. In this framework, Gordeyev integral Js(ω) ≡ ˆ ∞ 0 dτ e−jωτ 〈ejk·∆rs〉 (A.1) is defined as the one-sided Fourier transform of the characteristic function 〈ejk·∆rs〉 of random particle displacement vector ∆rs for species s over time intervals τ (in a plasma with suppressed collective interactions). Different types of plasmas are distinguished by different types of ∆rs statistics. Such statistics depend on the physical processes govern- ing individual particle motions and determine the species conductivities σs(k, ω) and the corresponding spectra of thermally impressed density fluctuations 〈|nts(k, ω)|2〉, functions that “collectively drive” the observed electron density spectrum 〈|ne(k, ω)|2〉. The char- acteristic function 〈ejk·∆rs〉 is also termed the single-particle signal correlation (or ACF), since signal returns from a single particle exposed to a radar pulse would be proportional to ejk·rs , with rs = rs(t) denoting the particle trajectory. In this appendix, we show that BGK-based incoherent scatter spectral models can also be derived from a particle dynamics approach and that they constitute no exception to the general framework presented in Chapter 2. In particular, the BGK-based Gordeyev integral is obtained, in Section A.2, by assuming a Poisson collision process, and it is shown that its inclusion in the general framework of Chapter 2 yields the usual BGK-model results [e.g., Dougherty , 1963; Dougherty and Farley , 1963] for the electron density spectrum. The appendix is concluded in Section A.3 with further discussions and implications of our results; for instance, an extension of our approach to include collisions with both neutrals and charged particles is proposed here. 128 A.2 Derivation of the BGK Gordeyev Integral following a Particle Dynamics Approach Clemmow and Dougherty [1969] state that the BGK model can be used to simulate the effect of binary collisions, e.g., collisions between charged and neutral particles, and also that these collisions could be imagined to be a Poisson process. However, this explanation was given more to provide a possible interpretation of the BGK model than to establish a direct relationship with the Poisson collision process. The mathematical proof that these two collision models are directly linked is provided next. Assume that in a time interval τ , a particle of type s (electron or ion) collides n times with neutral particles constituting a medium. If collision events are independent from each other and particles move in straight-line orbits in between collisions, we can then express the particle displacement over interval τ as the random vector ∆rs = n∑ l=0 vl(tl+1 − tl), (A.2) where t0 ≡ 0, tn+1 ≡ τ, and tl is the time of l-th collision such that 0 < t1 < ... < tn < τ . Let the number of collision events n invoked above be a Poisson random variable with a collision frequency ν and pmf p(n) = e−ντ (ντ)n n! for n ≥ 0. (A.3) Given that there are n collision events in an interval τ , the conditional pdf of collision times tl (for 1 ≤ l ≤ n) is given by [e.g., Hajek , 2009] f(t1, .., tn|n) =  n! τn if 0 < t1 < ... < tn < τ, 0 else. (A.4) This is in effect a uniform distribution over all ordered sets of collision times 0 < t1 < ... < 129 tn < τ , which can also be written as f(t1, .., tn|n) = n! τn n∏ l=0 u(tl+1 − tl) (A.5) in terms of unit-step function u(t) and satisfies ´ dtn · · · ´ dt1f(t1, .., tn|n) = 1 as all pdf’s do. Additionally, assume that the particle velocities vl between collision events, for l ∈ [0, n], constitute a set of independent and identically distributed random variables. Taking the distribution of vector velocities vl as Maxwellian, we have f(vl) = 1 (2piC2)3/2 e− 1 2C2 v2l , (A.6) where C ≡ √ KT m is the thermal speed of the particles. Let us next compute the single-particle ACF 〈ejk·∆rs〉, the characteristic function of displacements ∆rs, required in the general framework outlined in Chapter 2, where the expected value will be computed over random variables v0, · · · ,vn, t1, · · · , tn, and n using the probability distributions defined above. We note that 〈ejk·∆rs〉 = 〈 ejk· ∑n l=0 vl(tl+1−tl) 〉 v,t,n = 〈 n∏ l=0 ejk·vl(tl+1−tl) 〉 v,t,n , (A.7) where the subscripts on the right indicate successive expected value operations to be per- formed. Starting with the expectations over independent Gaussian random variables vl, we have 〈ejk·∆rs〉 = 〈 n∏ l=0 e− 1 2 k2C2(tl+1−tl)2 〉 t,n . (A.8) Expectations over t1, · · · , tn next yield〈 n! τn ˆ dtn · · · ˆ dt1 n∏ l=0 e− 1 2 k2C2(tl+1−tl)2u(tl+1 − tl) 〉 n , (A.9) 130 where integration limits run from −∞ to ∞. Finally, using p(n), we find that the ACF is ∞∑ n=0 νne−ντ ˆ dtn · · · ˆ dt1 n∏ l=0 e− 1 2 k2C2(tl+1−tl)2u(tl+1 − tl). (A.10) Now we define g(t) ≡ e−νt− 12k2C2t2 u(t), (A.11) to rearrange (A.10) as 〈ejk·∆rs〉 = ∞∑ n=0 νn ˆ dtn · · · ˆ dt1 n∏ l=0 g(tl+1 − tl) = ∞∑ n=0 νn ˆ dtng(τ − tn) · · · ˆ dt1g(t2 − t1)g(t1). (A.12) Clearly, the integral chain above is n successive convolutions of n + 1 realizations of g(t) evaluated at t = τ . Thus, the ACF of particles of type s with a Maxwellian velocity distribution undergoing a Poisson collision process reduces to 〈ejk·∆rs〉 = ∞∑ n=0 νn g(τ) ∗ · · · ∗ g(τ)︸ ︷︷ ︸ n+1 . (A.13) Since 〈ejk·∆rs〉 is the characteristic function of ∆rs, i.e., the Fourier transform of the pdf of ∆rs, our particle displacement model with Poisson collisions is uniquely described by (A.13). The corresponding Gordeyev integral (A.1), i.e., a one-sided Fourier transform of the single-particle ACF (A.13), is then (using the properties of the convolution and the Fourier transform) Js(ω) = G(ω) ∞∑ n=0 νnGn(ω), (A.14) where G(ω) ≡ ∞ˆ 0 dτe−jωτe−ντ− 1 2 k2C2τ2 (A.15) 131 is the Fourier transform of g(t). The convergence of the series in (A.14) is guaranteed as ∞∑ n=0 νnGn(ω) = 1 1− νG(ω) (A.16) since |νG(ω)| < 1 for any ν ≥ 0. Hence, Gordeyev integral (A.14) reduces to Js(ω) = G(ω) 1− νG(ω) , (A.17) a well-known result previously derived by Dougherty [1963] using the BGK collision oper- ator and Boltzmann kinetic equation for unmagnetized plasmas. A.3 Discussion We can argue, based on our result in Section A.2, that from a particle dynamics per- spective, the BGK collision process is a Poisson process. As a consequence, any of the expected quantities that can be obtained using the BGK kinetic equation can also be de- rived using our stochastic model for particle dynamics, for example, species conductivities σs(k, ω) and the spectra of thermally impressed density fluctuations 〈|nts(k, ω)|2〉. As we have shown above, the general framework for incoherent scatter spectrum models presented in Chapter 2 can be developed independent of plasma kinetic equations on the basis of only the following fundamental relations: the fluctuation-dissipation or Nyquist theorem [e.g., Callen and Welton, 1951], the Kramers-Kronig relations [e.g., Clemmow and Dougherty , 1969; Chew , 1990], and the Maxwell equations. Therefore, plasmas can also be studied based on these principles; whether we use kinetic equations or the particle dynamics approach in solving plasma problems concerning particles in thermal equilibrium (including the phenomenon of Landau damping) is a matter of choice. Let us now calculate, as an example, the covariance matrix of particle displacements in a Poisson collision process. The covariance matrix is defined as 〈∆rs∆rTs 〉 = 〈 n∑ l=0 n∑ l′=0 vlv T l′ (tl+1 − tl)(tl′+1 − tl′) 〉 v,t,n . (A.18) 132 Since vl are independent random variables, we have 〈∆rs∆rTs 〉 = 〈 n∑ l=0 vlv T l (tl+1 − tl)2 〉 v,t,n . (A.19) Furthermore, provided that the components of vl are Gaussian and independent, the co- variance matrix of vl is given by 〈vlvTl 〉 = C2 I, (A.20) where I is the identity matrix. Thus, 〈∆rs∆rTs 〉 = C2 I 〈 n∑ l=0 (tl+1 − tl)2 〉 t,n . (A.21) 133 After some math, it can be verified that1 〈 n∑ l=0 (tl+1 − tl)2 〉 t = 2τ2 n+ 2 . (A.22) Taking the expected value of (A.22) with respect to n and substituting in (A.21) gives 〈∆rs∆rTs 〉 = C2 I ∞∑ n=0 e−ντ (ντ)n n! 2τ2 n+ 2 = 2C2 ν2 I ∞∑ n=0 e−ντ (n+ 1)(ντ)n+2 (n+ 2)! , (A.23) which simplifies as 〈∆rs∆rTs 〉 = 2C2 ν2 I ( ντ − 1 + e−ντ) . (A.24) This is a very interesting result because the same mathematical expression can be 1Let sn(τ) ≡ 〈∑n l=0(tl+1 − tl)2 〉 t for n ≥ 0. Taking expectations and expanding terms, we find that sn(τ) = n! τn τˆ 0 dtn · · · t2ˆ 0 dt1 n∑ l=0 (tl+1 − tl)2 = n! τn τˆ 0 dtn · · · t2ˆ 0 dt1(τ − tn)2 + n! τn τˆ 0 dtn · · · t2ˆ 0 dt1 n−1∑ l=0 (tl+1 − tl)2 = n τn τˆ 0 dtn(τ − tn)2tn−1n + n τn τˆ 0 dtnt n−1 n (n− 1)! tn−1n tnˆ 0 dtn−1 · · · t2ˆ 0 dt1 n−1∑ l=0 (tl+1 − tl)2 ︸ ︷︷ ︸ sn−1(tn) = 2τ2 (n+ 1) (n+ 2) + n τn τˆ 0 dt tn−1sn−1(t). Using this recursive equation, we can compute s0(τ) = τ 2, s1(τ) = τ2 3 + 1 τ τˆ 0 dt τ2 = 2 3 τ2, s2(τ) = τ2 6 + 2 τ2 τˆ 0 dt 2 3 τ3 = 1 2 τ2, s3(τ) = τ2 10 + 3 τ3 τˆ 0 dt 1 2 τ4 = 2 5 τ2, and so forth. Then, by induction, we obtain sn(τ) = 2τ2 n+ 2 . 134 derived in the context of a Brownian-motion collision model [e.g., Chandrasekhar , 1943], an approach that it is often used to describe Coulomb collisions [e.g., Zagorodny and Holod , 2000]. In the Brownian-motion formalism, the effects of collisions on particle motion are considered to be caused by the action of a friction force and random diffusion forces. A parameter analogous to ν is also defined, but it is regarded as a friction coefficient. Although the expressions for 〈∆rs∆rTs 〉 are the same for the Poisson and the Brownian- motion models, the corresponding expressions for the single-particle ACF’s are not equal and lead to different incoherent scatter spectral shapes [e.g., Hagfors and Brockelman, 1971]. The differences, however, are only noticeable for intermediate values of ν since both spectra converge to the same asymptotic expressions in the collisionless/frictionless (ν → 0) as well as high collision/friction (ν → ∞) limits. The difference at intermediate ν can be attributed to ∆rs(τ) being a Gaussian random variable at each τ in case of a Brownian-motion process, but only so in τ → 0 and ∞ limits for a Poisson process. Note that when ∆rs(τ) is strictly Gaussian, and only then, the ACF 〈ejk·∆rs〉 can be shown to reduce to e− 1 2 k2〈∆r2s〉, where 〈∆r2s〉 is a diagonal element of 〈∆rs∆rTs 〉. A generalization of the results presented in Section A.2 can be easily performed. Notice that the assumption that between collisions the particles move in straight-line orbits was not a necessary condition and it was only considered in order to recover the classical results of the BGK collision model. In between neutral collision events, we could have considered the particles moving in helical orbits due the action of an external magnetic force or even move randomly because of Coulomb interactions with other charged particles constituting a plasma. Either of these assumptions would have led to different definitions of the function G(ω), but it can be shown that the form of the Gordeyev integral Js(ω), i.e., equation (A.14), would have remained the same. In general, it is found that G(ω) = ∞ˆ 0 dτe−jωτe−ντ 〈ejk·∆ris〉, (A.25) where 〈ejk·∆ris〉 is the single-particle ACF of the process that takes place between Poisson collision events. For instance, let us consider an unmagnetized plasma in which both 135 neutral and Coulomb collisions are relevant. Modeling the neutral collisions as a Poisson process of frequency ν, and Coulomb collisions as a Brownian-motion process with friction coefficient β, function G(ω) for a particle undergoing both types of collisions takes the following form G(ω) = ∞ˆ 0 dτe−jωτe−ντe− k2C2 β2 (βτ−1+e−βτ) . (A.26) This expression together with (A.14) provides us with a model for ionospheric incoherent scatter spectrum measurements detected from regions in which both neutral and Coulomb collisions are expected to be important, e.g., the equatorial 150-km region. More general extensions of the Poisson collision model can be pursued; for instance, a velocity dependent collision frequency ν(v) could be considered into the theory. Further generalizations of this type will be the subject of future studies. 136 APPENDIX B FRICTION AND DIFFUSION COEFFICIENTS FOR “TEST PARTICLES” IN NON-ISOTHERMAL PLASMAS In this appendix, it is shown that the steady-state solution of the Fokker-Planck equa- tion with the standard Spitzer coefficients for a plasma do not relax toward a Maxwellian distribution in the case of unequal electron and ion temperatures. However, if typical F-region plasma configurations are considered, deviations of the steady-state distributions from exact Maxwellian are very small. For instance, we have verified that the second order moments of the distributions (i.e., effective temperature values) have negligible dependence on the difference between the values of Te and Ti considered for the plasma. In addition, a modified version of the Spitzer friction coefficient is formulated in order to obtain steady- state Maxwellian distributions even in the case of unequal electron and ion temperatures (Te 6= Ti). B.1 Introduction Some of the problems in plasma physics can be reduced to the statistical analysis of the dynamics of one of the particles constituting the medium — a “test” particle problem [e.g., Rostoker and Rosenbluth, 1960]. For this analysis, a probability density function fs of a particle of type s is defined such that fs(r,v, t)drdv is the probability of finding the particle at a position r with velocity v at a time t. The time evolution of the probability 137 function fs is governed by the Boltzmann kinetic equation ∂fs ∂t + v · ∇rfs + F m · ∇vfs = ( δfs δt ) c , (B.1) where ( δfs δt ) c accounts for the changes in the distribution due to particle collisions. Note that fs must be properly normalized such that ´ fs(r,v, t)drdv = 1. If at a time to, position ro and velocity vo are set, the density function fs becomes a conditional pdf fs(r,v, t|ro,vo, to) (also referred as a transition probability function). This constitutes an initial value problem, since equation (B.1) is solved subject to the condition fs(r,v, to|ro,vo, to) = δ(r− ro)δ(v − vo). (B.2) The solution of this problem is of principal interest in plasma physics, and it can be used, for instance, to calculate the spectrum of electron density fluctuations [e.g., Rosenbluth and Rostoker , 1962; Woodman, 1967]. If Coulomb interactions are the dominant collision process, the right-hand side of (B.1) can take the form of a Fokker-Planck collision operator ( δfs δt ) c = −∇v · [〈 ∆v ∆t 〉 fs ] + 1 2 ∇v∇v : [〈 ∆v∆vT ∆t 〉 fs ] , (B.3) where ∇v is the gradient operator in velocity space, ∇v = xˆ (∂/∂vx) + yˆ (∂/∂vy) + zˆ (∂/∂vz) . 1 Above, 〈 ∆v ∆t 〉 and 〈 ∆v∆vT ∆t 〉 are the dynamical friction vector and velocity diffusion tensor, functions of particle velocity that can be formulated in terms of the Rosenbluth potentials [Rosenbluth et al., 1957]. A typical exercise in plasma books is the evaluation of 〈 ∆v ∆t 〉 and 〈 ∆v∆vT ∆t 〉 for the case of a multiple-component plasma in thermal equilibrium [e.g., Montgomery and Tidman, 1964; Gurnett and Bhattacharjee, 2005]. The resultant friction and diffusion coefficients are known as the Spitzer coefficients [Spitzer , 1The differential operations in equation (B.3) are defined as ∇v ·A(v) = ∑ i ∂ ∂vi Ai(v) and ∇v∇v : B¯(v) = ∑ i ∑ j ∂ ∂vi ∂ ∂vj Bi,j(v), where A(v) is a vector function and B¯(v) is a tensor. 138 1962]. Explicit expressions for these functions are provided in Section 2.5. In the literature, the Fokker-Planck collision model with the Spitzer coefficients is used to study the effects of Coulomb collisions on particle distributions in non-isothermal plasmas, i.e., in plasmas in which electron and ion temperatures are not equal [e.g., Milla and Kudeki , 2010]. Although this can be considered a valid extension of the model, we found that there is a discrepancy between the assumptions made to derive the coefficients of the Fokker-Planck equation and the solutions obtained from it. As we discuss in Section B.2, steady-state solutions of the Fokker-Planck equation (B.3) with the standard Spitzer coefficients do not relax towards Maxwellian distributions if electron and ion temperatures are not equal. However, to derive the Spitzer coefficients, it is assumed that the particle distributions for the different species in a non-isothermal plasma are Maxwellian. Despite the fact that deviations of the steady-state solutions of (B.3) from Maxwellian distributions are very small, this issue brings up the question of the validity of the model. In Section B.3, a modified version of the Spitzer friction coefficient is derived. This modified expression makes the steady-state solutions of the Fokker-Planck equation be exactly Maxwellian (even if electron and ion temperatures are different). Whether such a correction can be justified based on more fundamental physical principles than the ones exposed here will be the subject of a future investigation. A short discussion of our results presented in Section B.4 concludes this appendix. B.2 Steady-State Solution of the Fokker-Planck Equation with Spitzer Collision Coefficients In the case of a homogeneous stationary plasma, the friction vector and diffusion tensor of the Fokker-Planck equation (B.3) can be expressed in the forms 〈 ∆v ∆t 〉 = −β(v)v (B.4) and 〈 ∆v∆vT ∆t 〉 = D⊥(v) 2 I + ( D‖(v)− D⊥(v) 2 ) vvT v2 , (B.5) 139 where friction coefficient β(v) and diffusion coefficients D‖(v) and D⊥(v) are only depen- dent on particle speed v. Note that particle diffusion is characterized by two coefficients, D‖(v) and D⊥(v), one for the diffusion in the direction parallel to the particle trajectory and the other for the diffusion in the perpendicular plane. Plugging expressions (B.4) and (B.5) into the Fokker-Planck equation (B.3), we get ( δfs δt ) c = ∇v · [βv fs] + 1 2 ∇v∇v : [( D⊥ 2 I + ( D‖ − D⊥ 2 ) vvT v2 ) fs ] . (B.6) The second term in the right-hand side of this equation can be expanded using the following tensor identities ∇v∇v : [ gI ] = ∇2vg (B.7) and ∇v∇v : [ g vvT v2 ] = 2 v2 g + 4 v2 v · ∇vg + vv T v2 : ∇v∇vg, (B.8) where g(v) is a scalar function. As a result, (B.6) becomes ( δfs δt ) c = ∇v · [βv fs] +∇2v [ D⊥ 4 fs ] + 1 v2 [( D‖ − D⊥ 2 ) fs ] + 2v v2 · ∇v [( D‖ − D⊥ 2 ) fs ] + 1 2 vvT v2 : ∇v∇v [( D‖ − D⊥ 2 ) fs ] . (B.9) In steady state, the density function fs is time independent; therefore, ( δfs δt ) c = 0. (B.10) In addition, since β(v), D‖(v), and D⊥(v) are isotropic functions of particle speed v, the steady-state solution of (B.9) should also be an isotropic function, i.e., fs = fs(v). Then, velocity derivatives in (B.9) can be written as ∇vg = v v ∂g ∂v , ∇v · g = 1 v2 ∂ ∂v ( v2g · vˆ) , ∇2vg = 1v2 ∂∂v ( v2 ∂g ∂v ) , (B.11) 140 and ∇v∇vg = ( I− vv T v2 ) 1 v ∂g ∂v + vvT v2 ∂2g ∂v2 , (B.12) leading to ( δfs δt ) c = 1 v2 ∂ ∂v ( v3βfs ) + 1 v2 ∂ ∂v ( v2 ∂ ∂v ( D⊥ 4 fs )) + 1 v2 (( D‖ − D⊥ 2 ) fs ) + 2 v ∂ ∂v (( D‖ − D⊥ 2 ) fs ) + 1 2 ∂2 ∂v2 (( D‖ − D⊥ 2 ) fs ) = 0, (B.13) which, after simplification, becomes ( δfs δt ) c = 1 2v2 ∂ ∂v ( ∂ ∂v ( v2D‖fs ) + ( 2v3β − vD⊥ ) fs ) = 0. (B.14) Thus, the problem of computing the steady-state distribution fs reduces to solving ∂ ∂v ( v2D‖fs ) + ( 2v3β − vD⊥ ) fs = C, (B.15) where C is a constant of particle speed v. Among the different possible values that C can take, we choose C = 0 because it can be shown that, in this case, the solution of (B.15) does not have a singularity at the origin (i.e., at v = 0). Then, equation (B.15) becomes ∂ ∂v ( v2D‖fs ) + ( 2v3β − vD⊥ ) fs = 0, (B.16) or, equivalently, ∂fs ∂v + 1 v2D‖ ( ∂ ∂v ( v2D‖ ) + 2v3β − vD⊥ ) fs = 0. (B.17) 141 The solution of this equation is simply2 fs(v) = 1 Ws exp − vˆ 0 p(v)dv  , (B.18) where p(v) = 1 v2D‖ ( ∂ ∂v ( v2D‖ ) + 2v3β − vD⊥ ) (B.19) and Ws is a scalar such that ´ fs(v)dv 3 = 1. For the Spitzer coefficients defined in Section 2.5, namely β(v) = ∑ s′ (1 + ms ms′ )Ns′Γs/s′ 1 C2s′ ψ(zs′) v , (B.20) D||(v) = ∑ s′ 2Ns′Γs/s′ ψ(zs′) v , (B.21) D⊥(v) = ∑ s′ 2Ns′Γs/s′ φ(zs′)− ψ(zs′) v , (B.22) we can show, after some algebra, that p(v) = v C2s (∑ s′ Ts Ts′ Ns′q 2 s′ ln Λss′ψ(zs′)∑ s′ Ns′q 2 s′ ln Λss′ψ(zs′) ) , (B.23) which can also be written in the following form p(v) = v C2s 1 + ∑s′ 6=sNs′e2s′ ln Λss′ ( Ts Ts′ − 1 ) ψ(zs′)∑ s′ Ns′e 2 s′ ln Λss′ψ(zs′)  . (B.24) Then, placing this expression in (B.18), we can readily find that the steady-state density 2A first order differential equation of the form ∂f(x) ∂x + p(x)f(x) = q(x) has the following solution f(x) = e−g(x) xˆ 0 q(x)eg(x)dx+ f(0)e−g(x), where g(x) = ´ x 0 p(x)dx. 142 function fs is equal to fs(v) = 1 Ws exp − 1 C2s v2 2 + vˆ 0 v ∑ s′ 6=sNs′e 2 s′ ln Λss′ ( Ts Ts′ − 1 ) ψ(zs′)∑ s′ Ns′e 2 s′ ln Λss′ψ(zs′) dv  . (B.25) Note that, if the different species constituting a plasma are in thermal equilibrium (i.e., if Ts = Ts′), the density function (B.25) simplifies to fs(v) = ( 1 2piC2s )3/2 exp ( − v 2 2C2s ) , (B.26) thus, the steady-state distributions for all plasma species are exactly Maxwellian (which is in agreement with the assumption made to derive the Spitzer coefficients). On the other hand, however, if the plasma is not isothermal, the steady-state distributions deviate from the Maxwellian form as described by (B.25). For instance, let consider a single-ion plasma where Te > Ti. In this case, the electron and ion steady-state distributions are given by fe(v) = 1 We exp ( − 1 C2e ( v2 2 + ( Te Ti − 1 ) he(v) )) (B.27) and fi(v) = 1 Wi exp ( − 1 C2i ( v2 2 + ( Ti Te − 1 ) hi(v) )) , (B.28) respectively, where functions he and hi are defined as he(v) ≡ vˆ 0 v ln Λeiψ(zi) ln Λeeψ(ze) + ln Λeiψ(zi) dv (B.29) hi(v) ≡ vˆ 0 v ln Λieψ(ze) ln Λieψ(ze) + ln Λiiψ(zi) dv. (B.30) An analysis of these functions shows that the deviations of fe and fi from Maxwellian distributions are very small. To quantify the deviations, we have used the steady-state density functions (B.27) and (B.28) to calculate effective electron and ion temperatures and, then, compared these results to the values of Te and Ti considered to specify the Spitzer coefficients. The effective particle temperatures have been computed using the 143 Te/Ti T e f f e − T e T e [% ] 1 2 3 4 5 −5 −4 −3 −2 −1 0 ×10−3 Te = 500 K Te = 1000 K Te = 1500 K Te = 2000 K Te = 2500 K Te = 3000 K Te = 3500 K Figure B.1 Relative difference between the effective electron temperature T effe computed using (B.31) and the value of Te considered in the calculation of the Spitzer coefficients for an oxygen plasma. The difference is displayed as a function of Te/Ti. Each curve corresponds to different constant values of Te. following expression T effs ≡ ms 3κ ˆ v2fs(v)dv 3 = 4pims 3κ ∞ˆ 0 v4fs(v)dv, (B.31) wherems is the mass of particles of type s (either electrons or ions) and κ is the Boltzmann constant. In Figure B.1, the relative difference between the effective electron temperature T effe computed using (B.31) and the temperature Te considered in the calculation of the Spitzer coefficients for an oxygen plasma is displayed as a function of the ratio Te/Ti. Note that the different curves in Figure B.1, corresponding to different constant values of Te, overlap each other. As we can see in this plot, T effe is lower than Te for any Te > Ti and the magnitude of the relative difference between these values increases with Te/Ti; however, the difference is very small (of the order of 10−3 percent). In the case of the ion temperatures (see Figure B.2), T effi is greater than Ti for any Te > Ti. Although the differences between T effi and Ti are larger than the ones found in the case of the electron temperatures, such differences are still small (of the order of one percent). Note that the 144 Te/Ti T e f f i − T i T i [% ] 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1 Ti = 500 K Ti = 1000 K Ti = 1500 K Ti = 2000 K Ti = 2500 K Ti = 3000 K Ti = 3500 K Figure B.2 Same as Figure B.1 but for the case of the effective ion temperature T effi . ion temperature difference reaches its maximum value around Te/Ti = 3. Finally, for larger values of Te/Ti, the effective temperature T effi becomes closer to Ti. B.3 Modified Spitzer Friction Coefficient In the previous section, we found that the steady-state solution of the Fokker-Planck equation (B.6) is given by fs(v) = 1 Ws exp − vˆ 0 p(v)dv  , (B.32) where p(v) has been defined in (B.19). If we want the distribution fs to be Maxwellian, function p(v) is constrained to be equal to v/C2s , where Cs is the thermal speed of particle species s. Then, placing this condition in (B.19), we have that v C2s = 1 v2D‖(v) ( ∂ ∂v ( v2D‖(v) ) + 2v3β − vD⊥(v) ) . (B.33) 145 Solving for the friction coefficient β(v), we obtain β(v) = 1 2v3 ( v3 C2s D‖(v) + vD⊥(v)− ∂ ∂v ( v2D‖(v) )) . (B.34) Using Spitzer’s definitions for the diffusion coefficients D‖(v) and D⊥(v) in (B.34), we can find, after some manipulation, the following expression for the friction coefficient β(v) = ∑ s′ Ns′Γs/s′ ( 1 + Ts′ Ts ms ms′ ) 1 C2s′ ψ(zs′) v , (B.35) which can also be written in the following form: β(v) = ∑ s′ Ns′Γs/s′ ( 1 + ms ms′ ) 1 C2s′ ψ(zs′) v + ∑ s′ Ns′Γs/s′ ( 1− Ts Ts′ ) 1 C2s ψ(zs′) v . (B.36) Above, the first sum is equal to the standard Spitzer friction coefficient (defined in Section B.2), while the second term is the “correction” needed to make the distribution exactly Maxwellian. B.4 Discussion Concerns regarding the applicability of the Fokker-Planck equation with the Spitzer coefficients to study the effects of Coulomb collisions on the shape of the incoherent scatter spectrum were raised during our studies. First, we were concerned about the compatibility of this collision model with the general framework of incoherent scatter spectrum theories presented in Chapter 2. As we have discussed earlier in this dissertation, one of the funda- mentals principles on which our general framework is based is the fluctuation-dissipation or generalized Nyquist theorem. The theorem relates the spectrum of thermally driven density fluctuations and the conductivity function of each plasma species. However, the simple form of the theorem given in Chapter 2 (see expression (2.10)) is only valid if the velocity distribution of the particle species is Maxwellian. Although a more general formu- lation of the theorem is available [Sitenko, 1999], its application would further complicate the spectrum model. As we have shown in Section B.2, the steady-state solution of the 146 Fokker-Planck equation with the Spitzer coefficients is exactly Maxwellian only in the case of strict thermal equilibrium (i.e., when Te = Ti). However, for typical F-region plasma configurations, it was verified that, if the electron and ion temperatures are not equal, the deviation of the steady-state distribution from a Maxwellian form is very small and can be neglected for any practical purpose. Therefore, the fluctuation-dissipation theorem, as stated in equation (2.10), can be considered a good approximation for the Fokker-Planck Spitzer collision model used in our studies. A second concern that drew our attention was whether ionization and recombination processes should be considered in our collision model. As we have already mentioned, Spitzer coefficients come from a more general formulation of the Fokker-Planck equation in terms of the so-called Rosenbluth potentials [Rosenbluth et al., 1957]. In deriving these potentials, it is assumed that the plasma is fully ionized; thus, particles in the medium are only allowed to experience Coulomb collisions. Since ionization and recombination are disregarded, Rosenbluth’s model predicts the relaxation of the plasma to strict ther- mal equilibrium. Then, the reader may question the use of Spitzer coefficients to study non-isothermal plasmas, because this state is the result of a balance between ionization, recombination, and transport phenomena, which are processes that have been neglected in our model. Since production and recombination are relatively slow processes, their effects on incoherent scatter radar signals are typically ignored. This can be considered a valid approach for most radar viewing directions except, maybe, for observations perpendicular to the geomagnetic field, in which case the correlation time of the incoherent scatter signal becomes of the order of one second. This issue will be investigated in future studies. In an attempt to develop a more appropriate collision model for non-isothermal plas- mas, we derived a modified expression for the Spitzer friction coefficient based on the assumption that the particle distributions in a stationary plasma are Maxwellian even in the case of unequal electron and ion temperatures. The derivation of the modified coef- ficient is given in Section B.3. As we have already mentioned, the first term in (B.36) is the standard expression for the Spitzer friction coefficient, while the second term is the “correction” needed to make the distribution exactly Maxwellian. If we associate the cor- 147 rection term with ionization and recombination processes, we can interpret this term as the heating (cooling) rate required by the electrons (ions) to sustain the Te > Ti station- ary state. Whether the modified Spitzer friction coefficient can be derived from a physical model for ionization and recombination processes will be the subject of future research. 148 REFERENCES Aponte, N., M. P. Sulzer, and S. A. González (2001), Correction of the Jicamarca electron- ion temperature ratio problem: Verifying the effect of electron Coulomb collision on the incoherent scatter spectrum, Journal of Geophysical Research, 106 (A11), 24785–24793, doi:10.1029/2001JA000103. Bhatnagar, P. L., E. P. Gross, and M. Krook (1954), A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Physical Review, 94 (3), 511–525, doi:10.1103/PhysRev.94.511. Bhattacharyya, S. (1998), An inverse theory approach to incoherent scatter drift and tem- perature estimation at the Jicamarca Radio Observatory, Ph.D. dissertation, University of Illinois at Urbana-Champaign, Urbana, Illinois, 1998. Budden, K. G. (1961), Radio Waves in the Ionosphere, Cambridge University Press, Cam- bridge, United Kingdom. Callen, H. B., and T. A. Welton (1951), Irreversibility and generalized noise, Physical Review, 83 (1), 34–40. Callen, J. D. (2006), Fundamentals of Plasma Physics, Chapter 2 – Coulomb Collisions, University of Wisconsin, Madison, http://homepages.cae.wisc.edu/~callen/book. html. Chandrasekhar, S. (1943), Stochastic problems in physics and astronomy, Reviews of Mod- ern Physics, 15 (1), 1–89. Chew, W. C. (1990), Waves and Fields in Inhomogeneous Media, Van Nostrand Reinhold, New York, NY. Clemmow, P. C., and J. P. Dougherty (1969), Electrodynamics of Particles and Plasmas, Addison-Wesley Publishing Co., Inc., Redwood City, California. Coleman, T. F., and Y. Li (1996), An interior trust region approach for nonlinear minimization subject to bounds, SIAM Journal on Optimization, 6 (2), 418–445, doi: 10.1137/0806023. Dougherty, J. P. (1963), The conductivity of a partially ionized gas in alternating electric fields, Journal of Fluid Mechanics, 16 (1), 126–137. Dougherty, J. P. (1964), Model Fokker-Planck equation for a plasma an its solution, The Physics of Fluids, 7 (11), 1788–1799. 149 Dougherty, J. P., and D. T. Farley (1960), A theory of incoherent scattering of radio waves by a plasma, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 259 (1296), 79–99. Dougherty, J. P., and D. T. Farley (1963), A theory of incoherent scattering of radio waves by a plasma 3. Scattering in a partly ionized gas, Journal of Geophysical Research, 68, 5473–5486. Farley, D. T. (1964), The effect of Coulomb collisions on incoherent scattering of ra- dio waves by a plasma, Journal of Geophysical Research, 69 (1), 197–200, doi:10.1029/ JZ069i001p00197. Farley, D. T. (1966), A theory of incoherent scattering of radio waves by a plasma 4. The effect of unequal ion and electron temperatures, Journal of Geophysical Research, 71 (17), 4091–4098, doi:10.1029/JZ071i017p04091. Farley, D. T. (1969), Faraday rotation measurements using incoherent scatter, Radio Sci- ence, 4 (2), 143–152, doi:10.1029/RS004i002p00143. Farley, D. T., J. P. Dougherty, and D. W. Barron (1961), A theory of incoherent scattering of radio waves by a plasma ii. Scattering in a magnetic field, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 263 (1313), 238–258. Fejer, J. A. (1960), Scattering of radio waves by an ionized gas in thermal equilibrium, Canadian Journal of Physics, 38, 1114–1133. Fejer, J. A. (1961), Scattering of radio waves by an ionized gas in thermal equilibrium in the presence of a uniform magnetic field, Canadian Journal of Physics, 39, 716–740. Feng, Z. (2002), Equatorial F-region plasma density estimation with incoherent scatter radar using a transverse-mode differential-phase method, Ph.D. dissertation, University of Illinois at Urbana-Champaign, Urbana, Illinois, 2002. Feng, Z., E. Kudeki, R. F. Woodman, and J. L. Chau (2003), Transverse-beam incoherent scatter radar measurements of F -region plasma densities at Jicamarca, Radio Science, 38 (4), 1071, doi:10.1029/2002RS002722. Feng, Z., E. Kudeki, R. F. Woodman, J. L. Chau, and M. A. Milla (2004), F region plasma density estimation at Jicamarca using the complex cross-correlation of orthogonal po- larized backscatter fields, Radio Science, 39 (RS3015), 1–8, doi:10.1029/2003RS002963. Friedrich, M., K. M. Torkar, G. A. Lehmacher, C. L. Croskey, J. D. Mitchell, E. Kudeki, and M. A. Milla (2006), Rocket and incoherent scatter radar common-volume elec- tron measurements of the equatorial lower ionosphere, Geophysical Research Letters, 33, L08807, doi:10.1029/2005GL024622. Gillespie, D. T. (1996a), The mathematics of Brownian motion and Johnson noise, Amer- ican Journal of Physics, 64 (3), 225–240. Gillespie, D. T. (1996b), The multivariate Langevin and Fokker–Planck equations, Amer- ican Journal of Physics, 64 (10), 1246–1257. 150 Gurnett, D. A., and A. Bhattacharjee (2005), Introduction to Plasma Physics with Space and Laboratory Applications, Cambridge University Press, Cambridge, United Kingdom. Hagfors, T. (1961), Density fluctuations in a plasma in a magnetic field, with applications to the ionosphere, Journal of Geophysical Research, 66 (6), 1699–1712. Hagfors, T., and R. A. Brockelman (1971), A theory of collision dominated electron den- sity fluctuations in a plasma with applications to incoherent scattering, The Physics of Fluids, 14 (6), 1143–1151, doi:10.1063/1.1693578. Hajek, B. (2009), An Exploration of Random Processes for Engineers, Department of Electrical and Computer Engineering at the University of Illinois at Urbana-Champaign, Urbana, Illinois. Holod, I., A. Zagorodny, and J. Weiland (2005), Anisotropic diffusion across an external magnetic field and large-scale fluctuations in magnetized plasmas, Physical Review E, 71 (4), 046401, doi:10.1103/PhysRevE.71.046401. Ichimaru, S., and M. N. Rosenbluth (1970), Relaxation processes in plasmas with magnetic field. Temperature relaxations, The Physics of Fluids, 13 (11), 2778–2789, doi:10.1063/ 1.1692864. Kubo, R. (1966), The fluctuation-dissipation theorem, Reports on Progress in Physics, 29 (1), 255–284. Kudeki, E., and M. A. Milla (2006), Incoherent scatter spectrum theory for modes propagating perpendicular to the geomagnetic field, Journal of Geophysical Research, 111 (A06306), 1–3, doi:10.1029/2005JA011546. Kudeki, E., and M. A. Milla (2010), Incoherent scatter spectral theories: 1. A general framework and results for small magnetic aspect angles, IEEE Transactions on Geo- science and Remote Sensing, in press. Kudeki, E., S. Bhattacharyya, and R. F. Woodman (1999), A new approach in incoher- ent scatter F region E × B drift measurements at Jicamarca, Journal of Geophysical Research, 104 (A12), 28145–28162, doi:10.1029/1998JA900110. Kudeki, E., R. F. Woodman, and Z. Feng (2003), Incoherent scatter radar plasma density measurements at Jicamarca using a transverse-mode differential-phase method, Geo- physical Research Letters, 30 (5), 1255, doi:10.1029/2002GL015496. Kudeki, E., M. A. Milla, M. Friedrich, G. Lehmacher, and D. Sponseller (2006), ALTAIR incoherent scatter observations of the equatorial daytime ionosphere, Geophysical Re- search Letters, 33 (L08108), 1–5, doi:10.1029/2005GL025180. Lehmacher, G. A., C. L. Croskey, J. D. Mitchell, M. Friedrich, F.-J. Lubken, M. Rapp, E. Kudeki, and D. C. Fritts (2006), Intense turbulence observed above mesospheric temperature inversion at equatorial latitude, Geophysical Research Letters, 33 (L08808), 1–5, doi:10.1029/2005GL024345. Levanon, N. (1988), Radar Principles, Wiley-Interscience, New York, NY. 151 Li, Y. L., C. H. Liu, and S. J. Franke (1991), Adaptive evaluation of the Sommerfeld-type integral using the Chirp z-Transform, IEEE Transactions on Antennas and Propagation, 39 (12), 1788–1791, doi:10.1109/8.121603. Marsaglia, G., and W. W. Tsang (2000), The Ziggurat method for generating random variables, Journal of Statistical Software, 5 (8), 1–7. Milla, M. A., and E. Kudeki (2006), F -region electron density and Te/Ti measurements using incoherent scatter power data collected at ALTAIR, Annales Geophysicae, 24 (5), 1333–1342. Milla, M. A., and E. Kudeki (2009), Particle dynamics description of “BGK collisions” as a Poisson process, Journal of Geophysical Research, 114 (7), A07302, doi:10.1029/ 2009JA014200. Milla, M. A., and E. Kudeki (2010), Incoherent scatter spectral theories: 2. Modeling the spectrum for modes propagating perpendicular to B, IEEE Transactions on Geoscience and Remote Sensing, in press. Montgomery, D. C., and D. A. Tidman (1964), Plasma Kinetic Theory, 1st ed., McGraw- Hill Book Company, New York, NY. Ochs, G. R. (1965), The large 50 mc/s dipole array at Jicamarca radar observatory, NBS Report 8772, National Bureau of Standards, 1965. Olsen, N., T. J. Sabaka, and L. Tøffner-Clausen (2000), Determination of the IGRF 2000 model, Earth Planets Space, 52 (12), 1175–1182. Rosenbluth, M. N., and N. Rostoker (1962), Scattering of electromagnetic waves by a nonequilibrium plasma, The Physics of Fluids, 5 (7), 776–788, doi:10.1063/1.1724446. Rosenbluth, M. N., W. M. MacDonald, and D. L. Judd (1957), Fokker-Planck equation for an inverse-square force, Physical Review, 107 (1), 1–6, doi:10.1103/PhysRev.107.1. Rostoker, N., and M. N. Rosenbluth (1960), Test particles in a completely ionized plasma, The Physics of Fluids, 3 (1), 1–14, doi:10.1063/1.1705998. Salpeter, E. E. (1960), Electron density fluctuations in a plasma, Physical Review, 120 (5), 1528–1535. Sitenko, A. G. (1999), Extension of the fluctuation-dissipation theorem to non-equilibrium systems, Physics Letters A, 252 (6), 336–339. Spitzer, L., Jr. (1962), Physics of Fully Ionized Gases, 2nd ed., Wiley-Interscience, New York, NY. Sulzer, M. P., and S. A. González (1999), The effect of electron Coulomb collisions on the incoherent scatter spectrum in the F region at Jicamarca, Journal of Geophysical Research, 104 (A10), 22535–22551, doi:10.1029/1999JA900288. Uhlenbeck, G. E., and L. S. Ornstein (1930), On the theory of the Brownian motion, Physical Review, 36 (5), 823–841. 152 Woodman, R. F. (1967), Incoherent scattering of electromagnetic waves by a plasma, Ph.D. dissertation, Harvard University, Cambridge, Massachusetts, 1967. Woodman, R. F. (2004), On a proper electron collision frequency for a Fokker-Planck col- lision model with Jicamarca applications, Journal of Atmospheric and Solar-Terrestrial Physics, 66 (17), 1521–1541, doi:10.1016/j.jastp.2004.07.001. Woodman, R. F., and T. Hagfors (1969), Methods for the measurement of vertical iono- spheric motions near the magnetic Equator by incoherent scattering, Journal of Geo- physical Research, 74 (5), 1205–1212, doi:10.1029/JA074i005p01205. Yeh, K. C., H. Y. Chao, and K. H. Lin (1999a), A study of the generalized Faraday effect in several media, Radio Science, 34 (1), 139–153, doi:10.1029/98RS02442. Yeh, K. C., H. Y. Chao, and K. H. Lin (1999b), Polarization transformation of a wave field propagating in an anisotropic medium, IEEE Antennas and Propagation Magazine, 41 (5), 19–33, doi:10.1109/74.801511. Zagorodny, A., and I. Holod (2000), Large-scale fluctuations and particle diffusion across external magnetic field in turbulent plasmas, Condensed Matter Physics, 3 (2), 295–306. 153 AUTHOR’S BIOGRAPHY Marco Antonio Milla was born to Graciela Bravo and Felipe Milla on April 16, 1975, in Lima, Perú. He finished his bachelor studies in electrical engineering at the Pontificia Universidad Católica del Perú in 1997. He joined the Jicamarca Radio Observatory as a signal-processing engineer in 1998, and worked there for about five years in different research and engineering projects. In 2004, he started his graduate studies joining the Remote Sensing and Space Sciences group of the Department of Electrical and Computer Engineering at the University of Illinois at Urbana-Champaign. He received his M.S. degree in 2006 and is currently pursuing the Ph.D. degree. His doctoral research involves the development of incoherent scatter radar techniques for the estimation of ionospheric physical parameters. In particular, he has studied Coulomb collisions and magneto-ionic propagation effects on the incoherent scatter spec- trum measured with antenna beams pointed perpendicular to the Earth’s magnetic field. 154